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_MatrixSymHessianRelaxNonSing.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 
44 #include "ConstrainedOptPack_MatrixSymHessianRelaxNonSing.hpp"
45 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
46 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
47 #include "AbstractLinAlgPack_SpVectorClass.hpp"
48 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
49 #include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"
50 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
51 #include "DenseLinAlgPack_AssertOp.hpp"
52 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
53 #include "ProfileHackPack_profile_hack.hpp"
54 #include "Teuchos_Assert.hpp"
55 
56 namespace {
57 
58 //
59 template<class V>
60 void Vp_StPtMtV_imp(
61  DenseLinAlgPack::DVectorSlice* y, DenseLinAlgPack::value_type a
64  ,const V& x, DenseLinAlgPack::value_type b
65  )
66 {
67  using BLAS_Cpp::no_trans;
68  using BLAS_Cpp::trans;
69  using BLAS_Cpp::trans_not;
74  namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
75 
76 #ifdef TEUCHOS_DEBUG
77  TEUCHOS_TEST_FOR_EXCEPT(y==NULL);
78 #endif
79 
81  no = H.G().rows(), // number of original variables
82  nr = H.M().rows(), // number of relaxation variables
83  nd = no + nr, // total number of variables
84  ydim = y->dim(); // y->dim() == number of rows in op(P)
85 
86  DenseLinAlgPack::Vp_MtV_assert_sizes(
87  y->dim(),P.rows(),P.cols(),P_trans
88  ,BLAS_Cpp::rows( nd, nd, H_trans) );
89  DenseLinAlgPack::Vp_MtV_assert_sizes(
90  BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
91  ,nd, nd, H_trans, x.dim() );
92 
93  //
94  // y = b*y + a * op(P) * H * x
95  //
96  // y = b*y + a * [op(P1) op(P2) ] * [ G 0 ] * [ x1 ]
97  // [ 0 M ] [ x2 ]
98  //
99  // =>
100  //
101  // y = b*y + a*op(P1)*G*x1 + a*op(P2)*H*x2
102  //
103  // For this to work op(P) must be sorted by column.
104  //
105  if( ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == BLAS_Cpp::no_trans )
106  || ( P.ordered_by() == GPMSIP::BY_COL && P_trans == BLAS_Cpp::trans )
107  || ( P.ordered_by() == GPMSIP::UNORDERED ) )
108  {
109  // Call the default implementation
110  //H.MatrixOp::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b);
111  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
112  return;
113  }
114  const DenseLinAlgPack::Range1D
115  o_rng(1,no),
116  r_rng(no+1,no+nr);
118  P1 = ( P.is_identity()
119  ? GenPermMatrixSlice(
120  P_trans == no_trans ? ydim : no
121  ,P_trans == no_trans ? no : ydim
122  ,GenPermMatrixSlice::IDENTITY_MATRIX )
123  : P.create_submatrix(o_rng,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
124  ),
125  P2 = ( P.is_identity()
126  ? GenPermMatrixSlice(
127  P_trans == no_trans ? ydim : nr
128  ,P_trans == no_trans ? nr : ydim
129  ,GenPermMatrixSlice::ZERO_MATRIX )
130  : P.create_submatrix(r_rng,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
131  );
132  const V
133  x1 = x(o_rng),
134  x2 = x(r_rng);
135  // y = b*y
136  LinAlgOpPack::Vt_S(y,b);
137  // y += a*op(P1)*G*x1
138  if( P1.nz() )
139  LinAlgOpPack::Vp_StPtMtV( y, a, P1, P_trans, H.G(), H_trans, x1, b );
140  // y += a*op(P2)*M*x2
141  if( P2.nz() )
142  LinAlgOpPack::Vp_StPtMtV( y, a, P2, P_trans, H.M(), H_trans, x2, 1.0 );
143 }
144 
145 } // end namespace
146 
147 namespace ConstrainedOptPack {
148 
150  : vec_space_(Teuchos::null)
151 {}
152 
154  const G_ptr_t &G_ptr
155  ,const vec_mut_ptr_t &M_diag_ptr
156  ,const space_ptr_t &space
157  )
158  : vec_space_(Teuchos::null)
159 {
160  initialize(G_ptr,M_diag_ptr,space);
161 }
162 
164  const G_ptr_t &G_ptr
165  ,const vec_mut_ptr_t &M_diag_ptr
166  ,const space_ptr_t &space
167  )
168 {
169  namespace mmp = MemMngPack;
170 #ifdef TEUCHOS_DEBUG
171  const char err_msg_head[] = "MatrixSymHessianRelaxNonSing::initialize(...) : Error!";
172  TEUCHOS_TEST_FOR_EXCEPTION(G_ptr.get()==NULL, std::invalid_argument, err_msg_head);
173  TEUCHOS_TEST_FOR_EXCEPTION(M_diag_ptr.get()==NULL, std::invalid_argument, err_msg_head);
174  const size_type G_rows = G_ptr->rows(), M_diag_dim = M_diag_ptr->dim();
175  TEUCHOS_TEST_FOR_EXCEPTION(G_rows==0, std::invalid_argument, err_msg_head);
176  TEUCHOS_TEST_FOR_EXCEPTION(M_diag_dim==0, std::invalid_argument, err_msg_head);
177 #endif
178  if( space.get() ) {
179 #ifdef TEUCHOS_DEBUG
180  const size_type space_dim = space->dim();
181  TEUCHOS_TEST_FOR_EXCEPTION(space_dim != G_rows + M_diag_dim, std::invalid_argument, err_msg_head);
182 #endif
183  vec_space_ = space;
184  }
185  else {
186  VectorSpace::space_ptr_t spaces[]
187  = { Teuchos::rcp(&G_ptr->space_cols(),false), Teuchos::rcp(&M_diag_ptr->space(),false) };
188  vec_space_ = Teuchos::rcp(new VectorSpaceBlocked( spaces, 2 ) );
189  }
190  G_ptr_ = G_ptr;
191  M_.initialize(M_diag_ptr);
192 }
193 
194 // Overridden from MatrixOp
195 
197 {
198  assert_initialized();
199  return *vec_space_;
200 }
201 
203  MatrixOp* C, value_type a, BLAS_Cpp::Transp H_trans
204  ) const
205 {
206 #ifdef PROFILE_HACK_ENABLED
207  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Mp_StM(...)" );
208 #endif
209  assert_initialized();
210  return MatrixOp::Mp_StM(C,a,H_trans); // ToDo: Update below code!
211 /*
212  const size_type
213  nG = G_ptr_->rows(),
214  nM = M_.rows();
215  AbstractLinAlgPack::Mp_StM( &(*C)(1,nG,1,nG), a, *G_ptr_, H_trans);
216  AbstractLinAlgPack::Mp_StM( &(*C)(nG+1,nG+nM,nG+1,nG+nM), a, M_, H_trans);
217 */
218 }
219 
221  VectorMutable* y, value_type a, BLAS_Cpp::Transp H_trans
222  ,const Vector& x, value_type b
223  ) const
224 {
225 #ifdef PROFILE_HACK_ENABLED
226  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Vp_StMtV(...Vector...)" );
227 #endif
228  assert_initialized();
229  const size_type
230  nG = G_ptr_->rows(),
231  nM = M_.rows();
233  AbstractLinAlgPack::Vp_StMtV( y->sub_view(1,nG).get(), a, *G_ptr_, H_trans, *x.sub_view(1,nG) );
234  AbstractLinAlgPack::Vp_StMtV( y->sub_view(nG+1,nG+nM).get(), a, M_, H_trans, *x.sub_view(nG+1,nG+nM) );
235 }
236 
238  VectorMutable* y, value_type a, BLAS_Cpp::Transp H_trans
239  ,const SpVectorSlice& x, value_type b
240  ) const
241 {
242 #ifdef PROFILE_HACK_ENABLED
243  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Vp_StMtV(...SpVectorSlice...)" );
244 #endif
245  assert_initialized();
246  const size_type
247  nG = G_ptr_->rows(),
248  nM = M_.rows();
249  AbstractLinAlgPack::Vt_S(y,b); // Takes care of b == 0.0 and y uninitialized
250  AbstractLinAlgPack::Vp_StMtV( y->sub_view(1,nG).get(), a, *G_ptr_, H_trans, x(1,nG) );
251  AbstractLinAlgPack::Vp_StMtV( y->sub_view(nG+1,nG+nM).get(), a, M_, H_trans, x(nG+1,nG+nM) );
252 }
253 
255  VectorMutable* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
256  ,BLAS_Cpp::Transp H_trans, const Vector& x, value_type b
257  ) const
258 {
259 #ifdef PROFILE_HACK_ENABLED
260  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Vp_StPtMtV(...Vector...)" );
261 #endif
262  assert_initialized();
263  //MatrixOp::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // Uncomment for this default implementation
266  Vp_StPtMtV_imp(&y_d(),a,P,P_trans,*this,H_trans,x_d(),b);
267 }
268 
270  VectorMutable* y, value_type a, const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
271  ,BLAS_Cpp::Transp H_trans, const SpVectorSlice& x, value_type b
272  ) const
273 {
274 #ifdef PROFILE_HACK_ENABLED
275  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Vp_StPtMtV(...SpVectorSlice...)" );
276 #endif
277  assert_initialized();
278  //MatrixOp::Vp_StPtMtV(y,a,P,P_trans,H_trans,x,b); // Uncomment for this default implementation
280  Vp_StPtMtV_imp(&y_d(),a,P,P_trans,*this,H_trans,x,b);
281 }
282 
283 // Overridden form MatrixSymOp
284 
286  MatrixSymOp* S, value_type a
287  ,EMatRhsPlaceHolder dummy_place_holder
288  ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
289  ,value_type b
290  ) const
291 {
292  using BLAS_Cpp::no_trans;
293  using BLAS_Cpp::trans;
294  using BLAS_Cpp::trans_not;
295  namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
296 #ifdef PROFILE_HACK_ENABLED
297  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::Mp_StPtMtP(...)" );
298 #endif
299  assert_initialized();
300 
301  MatrixSymOp::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b); // ToDo: Override when needed!
302  return;
303 /* ToDo: Update below code!
304  const DenseLinAlgPack::size_type
305  no = G().rows(), // number of original variables
306  nr = M().rows(), // number of relaxation variables
307  nd = no + nr; // total number of variables
308 
309  DenseLinAlgPack::Mp_MtM_assert_sizes( S->rows(), S->cols(), no_trans
310  , P.rows(), P.cols(), trans_not(P_trans)
311  , P.rows(), P.cols(), P_trans );
312  DenseLinAlgPack::Vp_V_assert_sizes( BLAS_Cpp::rows( P.rows(), P.cols(), P_trans), nd );
313 
314  //
315  // S = b*S + a * op(P)' * H * op(P)
316  //
317  // S = b*S + a * [op(P1)' op(P2)' ] * [ G 0 ] * [ op(P1) ]
318  // [ 0 M ] [ op(P2) ]
319  //
320  // =>
321  //
322  // S = b*S
323  // S1 += op(P1)' * G * op(P1)
324  // S2 += op(P2)' * M * op(P2)
325  //
326  // For this to work op(P) must be sorted by row.
327  //
328  if( ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == BLAS_Cpp::trans )
329  || ( P.ordered_by() == GPMSIP::BY_COL && P_trans == BLAS_Cpp::no_trans )
330  || ( P.ordered_by() == GPMSIP::UNORDERED ) )
331  {
332  // Call the default implementation
333  MatrixSymOp::Mp_StPtMtP(S,a,dummy_place_holder,P,P_trans,b);
334  return;
335  }
336  const DenseLinAlgPack::Range1D
337  o_rng(1,no),
338  r_rng(no+1,no+nr);
339  const AbstractLinAlgPack::GenPermMatrixSlice
340  P1 = ( P.is_identity()
341  ? GenPermMatrixSlice(
342  P_trans == no_trans ? nd : no
343  ,P_trans == no_trans ? no : nd
344  ,GenPermMatrixSlice::IDENTITY_MATRIX )
345  : P.create_submatrix(o_rng,P_trans==no_trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
346  ),
347  P2 = ( P.is_identity()
348  ? GenPermMatrixSlice(
349  P_trans == no_trans ? nd : nr
350  ,P_trans == no_trans ? nr : nd
351  ,GenPermMatrixSlice::ZERO_MATRIX )
352  : P.create_submatrix(r_rng,P_trans==no_trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
353  );
354  // S = b*S
355  DenseLinAlgPack::Mt_S( &DMatrixSliceTriEle(S->gms(),S->uplo()),b); // Handles b == 0.0 properly!
356 
357  // S1 += a*op(P1)'*G*op(P1)
358  if( P1.nz() )
359  AbstractLinAlgPack::Mp_StPtMtP(
360  &DMatrixSliceSym( S->gms()(1,no,1,no), S->uplo() )
361  , a, dummy_place_holder, G(), P1, P_trans );
362  // S2 += a*op(P2)'*M*op(P2)
363  if( P2.nz() )
364  AbstractLinAlgPack::Mp_StPtMtP(
365  &DMatrixSliceSym( S->gms()(no+1,nd,no+1,nd), S->uplo() )
366  , a, dummy_place_holder, M(), P2, P_trans );
367 */
368 }
369 
370 // Overridden from MatrixOpNonsing
371 
373  VectorMutable* y, BLAS_Cpp::Transp H_trans, const Vector& x
374  ) const
375 {
376 #ifdef PROFILE_HACK_ENABLED
377  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::V_InvMtV(...Vector...)" );
378 #endif
379  assert_initialized();
380  const size_type
381  nG = G_ptr_->rows(),
382  nM = M_.rows();
383  AbstractLinAlgPack::V_InvMtV( y->sub_view(1,nG).get(), *G_ptr_, H_trans, *x.sub_view(1,nG) );
384  AbstractLinAlgPack::V_InvMtV( y->sub_view(nG+1,nG+nM).get(), M_, H_trans, *x.sub_view(nG+1,nG+nM) );
385 }
386 
388  VectorMutable* y, BLAS_Cpp::Transp H_trans, const SpVectorSlice& x
389  ) const
390 {
391 #ifdef PROFILE_HACK_ENABLED
392  ProfileHackPack::ProfileTiming profile_timing( "MatrixSymHessianRelaxNonSing::V_InvMtV(...SpVectorSlice...)" );
393 #endif
394  assert_initialized();
395  const size_type
396  nG = G_ptr_->rows(),
397  nM = M_.rows();
398  AbstractLinAlgPack::V_InvMtV( y->sub_view(1,nG).get(), *G_ptr_, H_trans, x(1,nG) );
399  AbstractLinAlgPack::V_InvMtV( y->sub_view(nG+1,nG+nM).get(), M_, H_trans, x(nG+1,nG+nM) );
400 }
401 
402 // private
403 
404 void MatrixSymHessianRelaxNonSing::assert_initialized() const
405 {
407  G_ptr_.get() == NULL, std::logic_error
408  ,"MatrixSymHessianRelaxNonSing::assert_initialized(): Error, Not initalized yet!" );
409 }
410 
411 } // end namespace ConstrainedOptPack
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
T * get() const
void initialize(const G_ptr_t &G_ptr, const vec_mut_ptr_t &M_diag_ptr, const space_ptr_t &space=Teuchos::null)
Initialize the Hessian and the relaxation terms.
void Mp_StPtMtP(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans)
void Mp_StPtMtP(MatrixSymOp *sym_lhs, value_type alpha, EMatRhsPlaceHolder dummy_place_holder, const GenPermMatrixSlice &gpms_rhs, BLAS_Cpp::Transp gpms_rhs_trans, value_type beta) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_t size_type
Matrix class for non-singular Hessian matrix augmented with a terms for "Big M" relaxation variables...
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
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)
Transp trans_not(Transp _trans)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta) const
Transp
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta) const
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
bool Mp_StM(MatrixOp *mwo_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)