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_DecompositionSystemVarReductImp.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 <typeinfo>
43 
44 #include "ConstrainedOptPack_DecompositionSystemVarReductImp.hpp"
45 #include "ConstrainedOptPack_MatrixIdentConcatStd.hpp"
46 #include "ConstrainedOptPack_MatrixVarReductImplicit.hpp"
47 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
48 #include "AbstractLinAlgPack_MatrixOpSubView.hpp"
49 #include "Teuchos_AbstractFactoryStd.hpp"
50 #include "Teuchos_dyn_cast.hpp"
51 #include "Teuchos_Assert.hpp"
52 
53 namespace ConstrainedOptPack {
54 
56  const VectorSpace::space_ptr_t &space_x
57  ,const VectorSpace::space_ptr_t &space_c
58  ,const basis_sys_ptr_t &basis_sys
59  ,const basis_sys_tester_ptr_t &basis_sys_tester
60  ,EExplicitImplicit D_imp
61  ,EExplicitImplicit Uz_imp
62  )
63  :DecompositionSystemVarReduct(D_imp,Uz_imp)
64  ,basis_sys_tester_(basis_sys_tester)
65  ,D_imp_used_(MAT_IMP_AUTO)
66 {
67  this->initialize(space_x,space_c,basis_sys);
68 }
69 
71  const VectorSpace::space_ptr_t &space_x
72  ,const VectorSpace::space_ptr_t &space_c
73  ,const basis_sys_ptr_t &basis_sys
74  )
75 {
76  namespace rcp = MemMngPack;
77  size_type num_vars = 0;
78 #ifdef TEUCHOS_DEBUG
80  basis_sys.get() != NULL && (space_x.get() == NULL || space_c.get() == NULL)
81  ,std::invalid_argument
82  ,"DecompositionSystemVarReductImp::initialize(...) : Error, "
83  "if basis_sys is set, then space_x and space_c must also be set!" );
84 #endif
85  if(basis_sys.get()) {
86  num_vars = basis_sys->var_dep().size() + basis_sys->var_indep().size();
87 #ifdef TEUCHOS_DEBUG
88  const size_type
89  space_x_dim = space_x->dim(),
90  space_c_dim = space_c->dim(),
91  num_equ = basis_sys->equ_decomp().size() + basis_sys->equ_undecomp().size();
93  num_vars != 0 && (space_x_dim != num_vars || space_c_dim != num_equ)
94  , std::invalid_argument
95  ,"DecompositionSystemVarReductImp::initialize(...) : Error, "
96  "the dimensions of space_x, space_c and basis_sys do not match up!" );
97 #endif
98  }
99  space_x_ = space_x;
100  space_c_ = space_c;
101  basis_sys_ = basis_sys;
102  if(num_vars) {
103  space_range_ = space_x_->sub_space(basis_sys->var_dep())->clone();
104  space_null_ = space_x_->sub_space(basis_sys->var_indep())->clone();
105  }
106  else {
107  space_range_ = Teuchos::null;
108  space_null_ = Teuchos::null;
109  }
110 }
111 
112 // Basis manipulation
113 
115  std::ostream *out
116  ,EOutputLevel olevel
117  ,ERunTests test_what
118  ,MatrixOp *Z
119  ,MatrixOp *Y
120  ,MatrixOpNonsing *R
121  ,MatrixOp *Uz
122  ,MatrixOp *Uy
124  ,Teuchos::RCP<MatrixOp> *D_ptr
125  )
126 {
127 
128  // ToDo: Implement undecomposed general equalities and general inequalities
129 
130  namespace rcp = MemMngPack;
131  using Teuchos::dyn_cast;
132 
133  if( out && olevel >= PRINT_BASIC_INFO )
134  *out << "\n****************************************************************"
135  << "\n*** DecompositionSystemVarReductImp::get_basis_matrices(...) ***"
136  << "\n****************************************************************\n";
137 
138  // ToDo: Validate input!
139 
140  //
141  // Get references to concrete matrix objects
142  //
143 
145  *Z_vr = Z ? &dyn_cast<MatrixIdentConcatStd>(*Z) : NULL;
146 
147  //
148  // Make all matrices but Z uninitialized to determine if we can
149  // reuse the Z.D matrix object (explicit or implicit).
150  // Also, return a reference to the basis matrix C from the
151  // R object.
152  //
153 
154  (*C_ptr) = uninitialize_matrices(out,olevel,Y,R,Uy);
155 
156  //
157  // Determine if we should be using an explicit or implicit D = -inv(C)*N object
158  // (if we are allowed to choose).
159  //
160  update_D_imp_used(&D_imp_used_);
161 
162  //
163  // Determine if we need to allocate a new matrix object for Z.D.
164  // Also, get a reference to the explicit D matrix (if one exits)
165  // and remove any reference to the basis matrix C by Z.D.
166  //
167 
168  bool new_D_mat_object = true; // Valgrind complains if this is not initialized.
169  if( Z_vr ) {
170  if( Z_vr->D_ptr().get() == NULL ) {
171  if( out && olevel >= PRINT_BASIC_INFO )
172  *out << "\nMust allocate a new matrix object for D = -inv(C)*N since "
173  << "one has not been allocated yet ...\n";
174  new_D_mat_object = true;
175  }
176  else {
178  *D_vr = dynamic_cast<MatrixVarReductImplicit*>(
179  const_cast<MatrixOp*>(Z_vr->D_ptr().get()) );
180  // We may have to reallocate a new object if someone is sharing it or
181  // if we are switching from implicit to explicit or visa-versa.
182  if( Z_vr->D_ptr().total_count() > 1 ) {
183  if( out && olevel >= PRINT_BASIC_INFO )
184  *out << "\nMust allocate a new matrix object for D = -inv(C)*N since someone "
185  << "else is using the current one ...\n";
186  new_D_mat_object = true;
187  }
188  else if( D_vr != NULL ) {
189  if( D_imp_used_ == MAT_IMP_EXPLICIT ) {
190  if( out && olevel >= PRINT_BASIC_INFO )
191  *out << "\nMust allocate a new matrix object for D = -inv(C)*N since we "
192  << "are switching from implicit to explicit ...\n";
193  new_D_mat_object = true;
194  }
195  }
196  else if( D_imp_used_ == MAT_IMP_IMPLICIT ) {
197  if( out && olevel >= PRINT_BASIC_INFO )
198  *out << "\nMust allocate a new matrix object for D = -inv(C)*N since we "
199  << "are switching from explicit to implicit ...\n";
200  new_D_mat_object = true;
201  }
202  // Remove the reference to the basis matrix C
203  if(D_vr)
204  D_vr->set_uninitialized();
205  }
206  }
207  else {
208  if( out && olevel >= PRINT_BASIC_INFO )
209  *out << "\nMust allocate a new matrix object for D = -inv(C)*N since "
210  << " Z == NULL was passed in ...\n";
211  new_D_mat_object = true;
212  }
213 
214  //
215  // Get the matrix object of D and allocate a new matrix object if needed.
216  //
217 
218  Teuchos::RCP<MatrixOp> _D_ptr;
219  if( new_D_mat_object ) {
220  // Create a new matrix object!
221  alloc_new_D_matrix( out, olevel, &_D_ptr );
222  }
223  else {
224  // Use current matrix object!
225  _D_ptr = Teuchos::rcp_const_cast<MatrixOp>(Z_vr->D_ptr());
226  }
227 
228  // Set cached implicit D matrix or external explicit D matrix
229  if(D_imp_used_ == MAT_IMP_IMPLICIT) {
230  D_ptr_ = _D_ptr; // Need to remember this for when update_decomp() is called!
231  if(D_ptr)
232  (*D_ptr) = Teuchos::null; // No explicit D matrix
233  }
234  else {
235  D_ptr_ = Teuchos::null; // This matrix will be passed back in set_basis_matrices(...) before update_decomp(...) is called.
236  if(D_ptr)
237  (*D_ptr) = _D_ptr; // Externalize explicit D matrix.
238  }
239 
240  //
241  // Determine if we need to allocate a new basis matrix object C.
242  //
243  // At this point, if a matrix object for C already exits and noone
244  // outside has a reference to this basis matrix object then the only
245  // references to it should be in C_ptr
246  //
247 
248  bool new_C_mat_object; // compiler should warn if used before initialized!
249  if( (*C_ptr).get() == NULL ) {
250  if( out && olevel >= PRINT_BASIC_INFO )
251  *out << "\nMust allocate a new basis matrix object for C since "
252  << "one has not been allocated yet ...\n";
253  new_C_mat_object = true;
254  }
255  else {
256  if( (*C_ptr).total_count() > 1 ) {
257  if( out && olevel >= PRINT_BASIC_INFO )
258  *out << "\nMust allocate a new basis matrix object C since someone "
259  << "else is using the current one ...\n";
260  new_C_mat_object = true;
261  }
262  else {
263  new_C_mat_object = false; // The current C matrix is owned only by us!
264  }
265  }
266 
267  //
268  // Get the basis matrix object C and allocate a new object for if we have to.
269  //
270 
271  if( new_C_mat_object) {
272  (*C_ptr) = basis_sys_->factory_C()->create();
273  if( out && olevel >= PRINT_BASIC_INFO )
274  *out << "\nAllocated a new basis matrix object C "
275  << "of type \'" << typeName(*(*C_ptr)) << "\' ...\n";
276  }
277 
278 
279  if( out && olevel >= PRINT_BASIC_INFO )
280  *out << "\nEnd DecompositionSystemVarReductImp::get_basis_matrices(...)\n";
281 
282 }
283 
285  std::ostream *out
286  ,EOutputLevel olevel
287  ,ERunTests test_what
288  ,const Teuchos::RCP<MatrixOpNonsing> &C_ptr
289  ,const Teuchos::RCP<MatrixOp> &D_ptr
290  ,MatrixOp *Uz
291  ,const basis_sys_ptr_t &basis_sys
292  )
293 {
294  C_ptr_ = C_ptr;
295  if( D_ptr.get() )
296  D_ptr_ = D_ptr;
297  if(basis_sys.get()) {
298  basis_sys_ = basis_sys;
299  space_range_ = space_x_->sub_space(basis_sys->var_dep())->clone();
300  space_null_ = space_x_->sub_space(basis_sys->var_indep())->clone();
301  }
302 }
303 
304 // Overridden from DecompositionSystem
305 
307 {
308  if(space_x_.get()) {
309  return space_x_->dim();
310  }
311  return 0; // Not fully initialized!
312 }
313 
315 {
316  if(space_c_.get()) {
317  return space_c_->dim();
318  }
319  return 0; // Not fully initialized!
320 }
321 
323 {
324  if(basis_sys_.get()) {
325  return basis_sys_->equ_decomp().size();
326  }
327  return 0; // Not fully initialized!
328 }
329 
330 // range and null spaces
331 
332 const VectorSpace::space_ptr_t
334 {
335  return space_range_;
336 }
337 
338 const VectorSpace::space_ptr_t
340 {
341  return space_null_;
342 }
343 
344 // matrix factories
345 
348 {
349  namespace rcp = MemMngPack;
350  return Teuchos::rcp(
352  );
353 }
354 
357 {
359 }
360 
362  std::ostream *out
363  ,EOutputLevel olevel
364  ,ERunTests test_what
365  ,const MatrixOp &Gc
366  ,MatrixOp *Z
367  ,MatrixOp *Y
368  ,MatrixOpNonsing *R
369  ,MatrixOp *Uz
370  ,MatrixOp *Uy
371  ,EMatRelations mat_rel
372  ) const
373 {
374  namespace rcp = MemMngPack;
375  using Teuchos::dyn_cast;
376 
377  if( out && olevel >= PRINT_BASIC_INFO ) {
378  *out << "\n***********************************************************"
379  << "\n*** DecompositionSystemVarReductImp::update_decomp(...) ***"
380  << "\n************************************************************\n";
381 
382  if(mat_rel != MATRICES_INDEP_IMPS)
383  *out << "\nWarning!!! mat_rel != MATRICES_INDEP_IMPS; The decompsition matrix "
384  << "objects may not be independent of each other!\n";
385  }
386 
387 #ifdef TEUCHOS_DEBUG
388  // Validate setup
390  basis_sys_.get() == NULL, std::logic_error
391  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
392  "a BasisSystem object has not been set yet!" );
393 #endif
394 
395  const size_type
396  n = this->n(),
397  m = this->m(),
398  r = this->r();
399  const Range1D
400  var_indep(r+1,n),
401  equ_decomp = this->equ_decomp();
402 
403 #ifdef TEUCHOS_DEBUG
404  // Validate input
406  Z==NULL&&Y==NULL&&R==NULL&&Uz==NULL&&Uy==NULL
407  , std::invalid_argument
408  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
409  "at least one of Z, Y, R, Uz and Uycan not be NULL!" );
411  m == r && Uz != NULL, std::invalid_argument
412  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
413  "Uz must be NULL if m==r is NULL!" );
415  m == r && Uy != NULL, std::invalid_argument
416  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
417  "Uy must be NULL if m==r is NULL!" );
418 #endif
419 
420  //
421  // Get references to concrete matrix objects
422  //
423 
425  *Z_vr = Z ? &dyn_cast<MatrixIdentConcatStd>(*Z) : NULL;
426  TEUCHOS_TEST_FOR_EXCEPT( !( Uz == NULL ) ); // ToDo: Implement undecomposed general equalities
427 
428  //
429  // Get smart pointers to unreferenced C and D matrix objects.
430  //
431 
434 
435  if( C_ptr_.get() ) {
436  //
437  // The function get_basis_matrices() was called by external client
438  // so we don't need to call it here or update the decomposition matrices.
439  //
440  C_ptr = C_ptr_; // This was set by set_basis_matrices()
441  }
442  else {
443  //
444  // Make all matrices uninitialized and get unreferenced smart
445  // pointers to C and D (explicit only!).
446  //
448  out,olevel,test_what,Z,Y,R,Uz,Uy,&C_ptr,&D_ptr);
449  }
450 
451  // Get the D matrix created by get_basis_matrices() and set by
452  // get_basis_matrices() if implicit or set by set_basis_matrices()
453  // if explicit. This matrix may not be allocated yet in which
454  // case we need to allocate it.
455  if( D_ptr.get() == NULL ) {
456  // D_ptr was not set in get_basis_matrix() in code above but
457  // it may be cashed (if explicit) in D_ptr.
458  if( D_ptr_.get() != NULL )
459  D_ptr = D_ptr_;
460  else
461  alloc_new_D_matrix( out, olevel, &D_ptr );
462  }
463 
464  if( C_ptr_.get() == NULL ) {
465 
466  //
467  // The basis matrices were not updated by an external client
468  // so we must do it ourselves here using basis_sys.
469  //
470 
471  if( out && olevel >= PRINT_BASIC_INFO )
472  *out << "\nUpdating the basis matrix C and other matices using the BasisSystem object ...\n";
473 
474  TEUCHOS_TEST_FOR_EXCEPT( !( D_ptr.get() ) ); // local programming error only!
475  TEUCHOS_TEST_FOR_EXCEPT( !( C_ptr.get() ) ); // local programming error only!
476 
477  basis_sys_->update_basis(
478  Gc // Gc
479  ,C_ptr.get() // C
480  ,D_imp_used_ == MAT_IMP_EXPLICIT ? D_ptr.get() : NULL // D?
481  ,NULL // GcUP == Uz
482  ,(mat_rel == MATRICES_INDEP_IMPS
483  ? BasisSystem::MATRICES_INDEP_IMPS
484  : BasisSystem::MATRICES_ALLOW_DEP_IMPS )
485  ,out && olevel >= PRINT_BASIC_INFO ? out : NULL
486  );
487  }
488 
489  //
490  // Create the matrix object: N = Gc(var_indep,cond_decomp)'
491  //
493  N_ptr = Teuchos::null;
494  if( D_imp_used_ == MAT_IMP_IMPLICIT ) {
496  GcDd_ptr = Gc.sub_view(var_indep,equ_decomp);
498  GcDd_ptr.get() == NULL, std::logic_error
499  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
500  "The matrix class used for the gradient of constraints matrix Gc of type \'"
501  << typeName(Gc) << "\' must return return.get() != NULL from "
502  "Gc.sub_view(var_indep,equ_decomp)!" );
503  if(mat_rel == MATRICES_INDEP_IMPS) {
504  GcDd_ptr = GcDd_ptr->clone();
506  GcDd_ptr.get() == NULL, std::logic_error
507  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
508  "The matrix class used for the gradient of constraints matrix Gc.sub_view(var_indep,equ_decomp) "
509  "of type \'" << typeName(*GcDd_ptr) << "\' must return return.get() != NULL from \n"
510  "Gc.sub_view(var_indep,equ_decomp)->clone() since mat_rel == MATRICES_INDEP_IMPS!" );
511  }
512  N_ptr = Teuchos::rcp(
513  new MatrixOpSubView(
514  Teuchos::rcp_const_cast<MatrixOp>(GcDd_ptr) // M_full
515  ,Range1D() // rng_rows
516  ,Range1D() // rng_cols
517  ,BLAS_Cpp::trans // M_trans
518  )
519  );
520  }
521 
522  //
523  // Test the basis matrix C and the other matrices.
524  //
525 
526  if( test_what == RUN_TESTS ) {
527  if( out && olevel >= PRINT_BASIC_INFO )
528  *out << "\nTesting the basis matrix C and other matices updated using the BasisSystem object ...\n";
530  basis_sys_tester_.get() == NULL, std::logic_error
531  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
532  "test_what == RUN_TESTS but this->basis_sys_tester().get() == NULL!" );
533  if( basis_sys_tester_->print_tests() == BasisSystemTester::PRINT_NOT_SELECTED ) {
534  BasisSystemTester::EPrintTestLevel
535  print_tests;
536  switch(olevel) {
537  case PRINT_NONE:
538  print_tests = BasisSystemTester::PRINT_NONE;
539  break;
540  case PRINT_BASIC_INFO:
541  print_tests = BasisSystemTester::PRINT_BASIC;
542  break;
543  case PRINT_MORE_INFO:
544  print_tests = BasisSystemTester::PRINT_MORE;
545  break;
546  case PRINT_VECTORS:
547  case PRINT_EVERY_THING:
548  print_tests = BasisSystemTester::PRINT_ALL;
549  break;
550  }
551  basis_sys_tester_->print_tests(print_tests);
552  basis_sys_tester_->dump_all( olevel == PRINT_EVERY_THING );
553  }
554  const bool passed = basis_sys_tester_->test_basis_system(
555  *basis_sys_ // basis_sys
556  ,&Gc // Gc
557  ,C_ptr.get() // C
558  ,N_ptr.get() // N
559  ,D_imp_used_ == MAT_IMP_EXPLICIT ? D_ptr.get() : NULL // D?
560  ,NULL // GcUP == Uz
561  ,out // out
562  );
564  !passed, TestFailed
565  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
566  "Test of the basis system failed!" );
567  }
568 
569  //
570  // Initialize the implicit D = -inv(C)*N matrix object.
571  //
572 
573  TEUCHOS_TEST_FOR_EXCEPT( !( D_ptr.get() ) ); // local programming error only?
574  if( D_imp_used_ == MAT_IMP_IMPLICIT ) {
575  if( !C_ptr.has_ownership() && mat_rel == MATRICES_INDEP_IMPS ) {
576  C_ptr = C_ptr->clone_mwons();
578  C_ptr.get() == NULL, std::logic_error
579  ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
580  "The matrix class used for the basis matrix C (from the BasisSystem object) "
581  "must return return.get() != NULL from the clone_mwons() method since "
582  "mat_rel == MATRICES_INDEP_IMPS!" );
583  }
585  C_ptr // C
586  ,N_ptr // N
587  ,Teuchos::null // D_direct
588  );
589  // Above, if the basis matrix object C is not owned by the smart
590  // reference counted pointer object C_ptr, then we need to make the
591  // reference that the implicit D = -inv(C)*N matrix object has to
592  // C independent and the clone() method does that very nicely.
593  }
594 
595  //
596  // Call on the subclass to construct Y, R and Uy given C and D.
597  //
598 
599  initialize_matrices(out,olevel,C_ptr,D_ptr,Y,R,Uy,mat_rel);
600 
601  //
602  // Reconstruct Z, Uz and Uy.
603  //
604 
605  if( Z_vr ) {
606  Z_vr->initialize(
607  space_x_ // space_cols
608  ,space_x_->sub_space(var_indep)->clone() // space_rows
609  ,MatrixIdentConcatStd::TOP // top_or_bottom
610  ,1.0 // alpha
611  ,D_ptr // D_ptr
612  ,BLAS_Cpp::no_trans // D_trans
613  );
614  }
615 
616  TEUCHOS_TEST_FOR_EXCEPT( !( Uz == NULL ) ); // ToDo: Implement for undecomposed general equalities
617 
618  // Clear cache for basis matrices.
619  C_ptr_ = Teuchos::null;
620  D_ptr_ = Teuchos::null;
621 
622  if( out && olevel >= PRINT_BASIC_INFO )
623  *out << "\nEnd DecompositionSystemVarReductImp::update_decomp(...)\n";
624 }
625 
627  std::ostream& out, const std::string& L
628  ) const
629 {
630  out
631  << L << "*** Variable reduction decomposition (class DecompositionSytemVarReductImp)\n"
632  << L << "C = Gc(var_dep,equ_decomp)' (using basis_sys)\n"
633  << L << "if C is nearly singular then throw SingularDecomposition exception\n"
634  << L << "if D_imp == MAT_IMP_IMPICIT then\n"
635  << L << " D = -inv(C)*N represented implicitly (class MatrixVarReductImplicit)\n"
636  << L << "else\n"
637  << L << " D = -inv(C)*N computed explicity (using basis_sys)\n"
638  << L << "end\n"
639  << L << "Z = [ D; I ] (class MatrixIdentConcatStd)\n"
640  << L << "Uz = Gc(var_indep,equ_undecomp)' - Gc(var_dep,equ_undecomp)'*D\n"
641  << L << "begin update Y, R and Uy\n"
642  ;
643  print_update_matrices( out, L + " " );
644  out
645  << L << "end update of Y, R and Uy\n"
646  ;
647 }
648 
649 // Overridden from DecompositionSystemVarReduct
650 
652 {
653  return basis_sys_.get() ? basis_sys_->var_indep() : Range1D::Invalid;
654 }
655 
657 {
658  return basis_sys_.get() ? basis_sys_->var_dep() : Range1D::Invalid;
659 }
660 
661 // protected
662 
664 {
665  *D_imp_used = ( D_imp() == MAT_IMP_AUTO
666  ? MAT_IMP_IMPLICIT // Without better info, use implicit by default!
667  : D_imp() );
668 }
669 
670 // private
671 
672 void DecompositionSystemVarReductImp::alloc_new_D_matrix(
673  std::ostream *out
674  ,EOutputLevel olevel
675  ,Teuchos::RCP<MatrixOp> *D_ptr
676  ) const
677 {
678  if(D_imp_used_ == MAT_IMP_IMPLICIT) {
679  (*D_ptr) = Teuchos::rcp(new MatrixVarReductImplicit());
680  if( out && olevel >= PRINT_BASIC_INFO )
681  *out << "\nAllocated a new implicit matrix object for D = -inv(C)*N "
682  << "of type \'MatrixVarReductImplicit\' ...\n";
683  }
684  else {
685  (*D_ptr) = basis_sys_->factory_D()->create();
686  if( out && olevel >= PRINT_BASIC_INFO )
687  *out << "\nAllocated a new explicit matrix object for D = -inv(C)*N "
688  << "of type \'" << typeName(*(*D_ptr)) << "\' ...\n";
689  }
690 }
691 
692 } // end namespace ConstrainedOptPack
DecompositionSystemVarReductImp(const VectorSpace::space_ptr_t &space_x, const VectorSpace::space_ptr_t &space_c, const basis_sys_ptr_t &basis_sys, const basis_sys_tester_ptr_t &basis_sys_tester, EExplicitImplicit D_imp, EExplicitImplicit Uz_imp)
Construct a variable reduction decomposition.
virtual void update_D_imp_used(EExplicitImplicit *D_imp_used) const
Update D_imp_used.
virtual mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t uninitialize_matrices(std::ostream *out, EOutputLevel olevel, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uy) const =0
Overridden by subclasses to uninitialized Y, R and Uy then return C if referenced.
const VectorSpace::space_ptr_t space_range() const
Returns this->space_x()->sub_space(var_dep)
void print_update_decomp(std::ostream &out, const std::string &leading_str) const
EOutputLevel
Enumeration for the amount of output to create from update_decomp().
bool has_ownership() const
virtual void initialize(const mat_nonsing_ptr_t &C, const mat_ptr_t &N, const mat_ptr_t &D_direct)
Initialize this matrix object.
ERunTests
Enumeration for if to run internal tests or not.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void initialize_matrices(std::ostream *out, EOutputLevel olevel, const mat_nonsing_fcty_ptr_t::element_type::obj_ptr_t &C_ptr, const mat_fcty_ptr_t::element_type::obj_ptr_t &D_ptr, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uy, EMatRelations mat_rel) const =0
Overridden by subclasses to initialize Y, R and Uy given C and D.
T_To & dyn_cast(T_From &from)
const basis_sys_ptr_t & basis_sys() const
T * get() const
const VectorSpace::space_ptr_t & space_c() const
void get_basis_matrices(std::ostream *out, EOutputLevel olevel, ERunTests test_what, MatrixOp *Z, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uz, MatrixOp *Uy, Teuchos::RCP< MatrixOpNonsing > *C_ptr, Teuchos::RCP< MatrixOp > *D_ptr)
Called by client to uninitialize decomposition matrices in prepairation for selecting a different bas...
void initialize(const VectorSpace::space_ptr_t &space_x, const VectorSpace::space_ptr_t &space_c, const basis_sys_ptr_t &basis_sys)
Initialize.
void update_decomp(std::ostream *out, EOutputLevel olevel, ERunTests test_what, const MatrixOp &Gc, MatrixOp *Z, MatrixOp *Y, MatrixOpNonsing *R, MatrixOp *Uz, MatrixOp *Uy, EMatRelations mat_rel) const
Preconditions:
virtual Range1D equ_decomp() const
Returns the range of the decomposed equalities.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void print_update_matrices(std::ostream &out, const std::string &leading_str) const =0
Print the sub-algorithm by which the matrices Y, R, Uy and Uy are updated.
Concrete implementation class for a matrix vertically concatonated with an identity matrix...
size_t size_type
Specialization of DecompositionSystem for variable reduction decompositions.
Teuchos::RCP< const Teuchos::AbstractFactory< MatrixOp > > mat_fcty_ptr_t
Implements D = - inv(C) * N for a variable reduction projection.
virtual void initialize(const VectorSpace::space_ptr_t &space_cols, const VectorSpace::space_ptr_t &space_rows, ETopBottom top_or_bottom, value_type alpha, const D_ptr_t &D_ptr, BLAS_Cpp::Transp D_trans)
Setup with a matrix object.
size_type r() const
Returns this->basis_sys()->equ_decomp().size().
virtual const D_ptr_t & D_ptr() const
Return the smart reference counted point to the D matrix.
const VectorSpace::space_ptr_t & space_x() const
const VectorSpace::space_ptr_t space_null() const
Returns this->space_x()->sub_space(var_indep)
Specialization node implementation subclass of DecompositionSystem for variable reduction decompositi...
int total_count() const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
virtual void set_uninitialized()
Set the matrix to uninitialized.
std::string typeName(const T &t)
void set_basis_matrices(std::ostream *out, EOutputLevel olevel, ERunTests test_what, const Teuchos::RCP< MatrixOpNonsing > &C_ptr, const Teuchos::RCP< MatrixOp > &D_ptr, MatrixOp *Uz, const basis_sys_ptr_t &basis_sys=Teuchos::null)
Set updated basis matrices along with a possibly updated basis system object.