ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_MatrixHessianRelaxed.cpp
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include <assert.h>
45 
46 #include "ConstrainedOptPack_MatrixHessianRelaxed.hpp"
47 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymOp.hpp"
48 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
49 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
50 #include "AbstractLinAlgPack_SpVectorOp.hpp"
51 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
52 
53 namespace LinAlgOpPack {
56 }
57 
58 namespace ConstrainedOptPack {
59 
61  :
62  n_(0)
63  ,H_(NULL)
64  ,bigM_(0.0)
65 {}
66 
68  const MatrixSymOp &H
69  , value_type bigM
70  )
71 {
72  n_ = H.rows();
73  H_ = &H;
74  bigM_ = bigM;
75 }
76 
77 // Overridden from Matrix
78 
80 {
81  return n_ ? n_ + 1 : 0;
82 }
83 
84 // Overridden from MatrixOp
85 
87  DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
88  , const DVectorSlice& x, value_type b ) const
89 {
90  using BLAS_Cpp::no_trans;
91  using BLAS_Cpp::trans;
93  //
94  // y = b*y + a * M * x
95  //
96  // = b*y + a * [ H 0 ] * [ x1 ]
97  // [ 0 bigM ] [ x2 ]
98  //
99  // =>
100  //
101  // y1 = b*y1 + a*H*x1
102  //
103  // y2 = b*y2 + bigM * x2
104  //
105  LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),rows(),cols(),M_trans,x.size());
106 
107  DVectorSlice
108  y1 = (*y)(1,n_);
109  value_type
110  &y2 = (*y)(n_+1);
111  const DVectorSlice
112  x1 = x(1,n_);
113  const value_type
114  x2 = x(n_+1);
115 
116  // y1 = b*y1 + a*H*x1
117  Vp_StMtV( &y1, a, *H_, no_trans, x1, b );
118 
119  // y2 = b*y2 + bigM * x2
120  if( b == 0.0 )
121  y2 = bigM_ * x2;
122  else
123  y2 = b*y2 + bigM_ * x2;
124 
125 }
126 
128  DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
129  , const SpVectorSlice& x, value_type b ) const
130 {
131  using BLAS_Cpp::no_trans;
132  using BLAS_Cpp::trans;
134  //
135  // y = b*y + a * M * x
136  //
137  // = b*y + a * [ H 0 ] * [ x1 ]
138  // [ 0 bigM ] [ x2 ]
139  //
140  // =>
141  //
142  // y1 = b*y1 + a*H*x1
143  //
144  // y2 = b*y2 + bigM * x2
145  //
146  LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),rows(),cols(),M_trans,x.size());
147 
148  DVectorSlice
149  y1 = (*y)(1,n_);
150  value_type
151  &y2 = (*y)(n_+1);
152  const SpVectorSlice
153  x1 = x(1,n_);
154  const SpVectorSlice::element_type
155  *x2_ele = x.lookup_element(n_+1);
156  const value_type
157  x2 = x2_ele ? x2_ele->value() : 0.0;
158 
159  // y1 = b*y1 + a*H*x1
160  Vp_StMtV( &y1, a, *H_, no_trans, x1, b );
161 
162  // y2 = b*y2 + bigM * x2
163  if( b == 0.0 )
164  y2 = bigM_ * x2;
165  else
166  y2 = b*y2 + bigM_ * x2;
167 
168 }
169 
171  DVectorSlice* y, value_type a
172  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
173  , BLAS_Cpp::Transp M_trans
174  , const DVectorSlice& x, value_type b ) const
175 {
176  using BLAS_Cpp::no_trans;
177  using BLAS_Cpp::trans;
178  namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
179  //
180  // y = b*y + a * op(P) * M * x
181  //
182  // = b*y + a * [ op(P1) op(P2) ] * [ H 0 ] * [ x1 ]
183  // [ 0 bigM ] [ x2 ]
184  //
185  // =>
186  //
187  // y = b*y + a*op(P1)*H*x1 + a*op(P2)*bigM*x2
188  //
189  LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans
190  , BLAS_Cpp::rows( rows(), cols(), M_trans) );
191  LinAlgOpPack::Vp_MtV_assert_sizes( BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
192  ,rows(),cols(),M_trans,x.size());
193 
194  // For this to work (as shown below) we need to have P sorted by
195  // row if op(P) = P' or sorted by column if op(P) = P. If
196  // P is not sorted properly, we will just use the default
197  // implementation of this operation.
198  if( ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == no_trans )
199  || ( P.ordered_by() == GPMSIP::BY_COL && P_trans == trans ) )
200  {
201  // Call the default implementation
202  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b);
203  return;
204  }
205 
206  if( P.is_identity() )
207  TEUCHOS_TEST_FOR_EXCEPT( !( BLAS_Cpp::rows( P.rows(), P.cols(), P_trans ) == n_ ) );
208 
209  const GenPermMatrixSlice
210  P1 = ( P.is_identity()
211  ? GenPermMatrixSlice( n_, n_, GenPermMatrixSlice::IDENTITY_MATRIX )
212  : P.create_submatrix(Range1D(1,n_),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
213  ),
214  P2 = ( P.is_identity()
215  ? GenPermMatrixSlice(
216  P_trans == no_trans ? n_ : 1
217  , P_trans == no_trans ? 1 : n_
218  , GenPermMatrixSlice::ZERO_MATRIX )
219  : P.create_submatrix(Range1D(n_+1,n_+1),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
220  );
221 
222  const DVectorSlice
223  x1 = x(1,n_);
224  const value_type
225  x2 = x(n_+1);
226  // y = b*y + a*op(P1)*H*x1
227  AbstractLinAlgPack::Vp_StPtMtV( y, a, P1, P_trans, *H_, no_trans, x1, b );
228  // y += a*op(P2)*bigM*x2
229  if( P2.nz() ){
230  TEUCHOS_TEST_FOR_EXCEPT( !( P2.nz() == 1 ) );
231  const size_type
232  i = P_trans == no_trans ? P2.begin()->row_i() : P2.begin()->col_j();
233  (*y)(i) += a * bigM_ * x2;
234  }
235 }
236 
238  DVectorSlice* y, value_type a
239  , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
240  , BLAS_Cpp::Transp M_trans
241  , const SpVectorSlice& x, value_type b ) const
242 {
243  using BLAS_Cpp::no_trans;
244  using BLAS_Cpp::trans;
245  namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
246  //
247  // y = b*y + a * op(P) * M * x
248  //
249  // = b*y + a * [ op(P1) op(P2) ] * [ H 0 ] * [ x1 ]
250  // [ 0 bigM ] [ x2 ]
251  //
252  // =>
253  //
254  // y = b*y + a*op(P1)*H*x1 + a*op(P2)*bigM*x2
255  //
256  LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans
257  , BLAS_Cpp::rows( rows(), cols(), M_trans) );
258  LinAlgOpPack::Vp_MtV_assert_sizes( BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
259  ,rows(),cols(),M_trans,x.size());
260 
261  // For this to work (as shown below) we need to have P sorted by
262  // row if op(P) = P' or sorted by column if op(P) = P. If
263  // P is not sorted properly, we will just use the default
264  // implementation of this operation.
265  if( ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == no_trans )
266  || ( P.ordered_by() == GPMSIP::BY_COL && P_trans == trans ) )
267  {
268  // Call the default implementation
269  MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b);
270  return;
271  }
272 
273  if( P.is_identity() )
274  TEUCHOS_TEST_FOR_EXCEPT( !( BLAS_Cpp::rows( P.rows(), P.cols(), P_trans ) == n_ ) );
275 
276  const GenPermMatrixSlice
277  P1 = ( P.is_identity()
278  ? GenPermMatrixSlice( n_, n_, GenPermMatrixSlice::IDENTITY_MATRIX )
279  : P.create_submatrix(Range1D(1,n_),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
280  ),
281  P2 = ( P.is_identity()
282  ? GenPermMatrixSlice(
283  P_trans == no_trans ? n_ : 1
284  , P_trans == no_trans ? 1 : n_
285  , GenPermMatrixSlice::ZERO_MATRIX )
286  : P.create_submatrix(Range1D(n_+1,n_+1),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
287  );
288 
289  const SpVectorSlice
290  x1 = x(1,n_);
291  const SpVectorSlice::element_type
292  *x2_ele = x.lookup_element(n_+1);
293  const value_type
294  x2 = x2_ele ? x2_ele->value() : 0.0;
295  // y = b*y + a*op(P1)*H*x1
296  AbstractLinAlgPack::Vp_StPtMtV( y, a, P1, P_trans, *H_, no_trans, x1, b );
297  // y += a*op(P2)*bigM*x2
298  if( P2.nz() ){
299  TEUCHOS_TEST_FOR_EXCEPT( !( P2.nz() == 1 ) );
300  const size_type
301  i = P_trans == no_trans ? P2.begin()->row_i() : P2.begin()->col_j();
302  (*y)(i) += a * bigM_ * x2;
303  }
304 }
305 
307  const SpVectorSlice& x1, BLAS_Cpp::Transp M_trans
308  , const SpVectorSlice& x2 ) const
309 {
310  using BLAS_Cpp::no_trans;
311  //
312  // a = x1' * M * x2
313  //
314  // = [ x11' x12' ] * [ H 0 ] * [ x21 ]
315  // [ 0 bigM ] [ x22 ]
316  //
317  // =>
318  //
319  // a = x11'*H*x21 + x12'*bigM*x22
320  //
321  DenseLinAlgPack::Vp_MtV_assert_sizes(x1.size(),rows(),cols(),M_trans,x2.size());
322 
323  if( &x1 == &x2 ) {
324  // x1 and x2 are the same sparse vector
325  const SpVectorSlice
326  x11 = x1(1,n_);
327  const SpVectorSlice::element_type
328  *x12_ele = x1.lookup_element(n_+1);
329  const value_type
330  x12 = x12_ele ? x12_ele->value() : 0.0;
331  return AbstractLinAlgPack::transVtMtV( x11, *H_, no_trans, x11)
332  + x12 * bigM_ * x12;
333  }
334  else {
335  // x1 and x2 could be different sparse vectors
336  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this!
337  }
338  return 0.0; // Will never be executed!
339 }
340 
341 } // end namespace ConstrainedOptPack
342 
343 #endif // 0
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
MatrixHessianRelaxed()
Construct to uninitialized.
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
size_t size_type
void initialize(const MatrixSymOp &H, value_type bigM)
Initialize.
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
void Vp_StPtMtV(DVectorSlice *vs_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_rhs2_trans, const DVectorSlice &vs_rhs3, value_type beta) const
Transp
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
value_type transVtMtV(const SpVectorSlice &sv_rhs1, BLAS_Cpp::Transp trans_rhs2, const SpVectorSlice &sv_rhs3) const