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"
53 namespace ConstrainedOptPack {
56 const VectorSpace::space_ptr_t &space_x
57 ,
const VectorSpace::space_ptr_t &space_c
59 ,
const basis_sys_tester_ptr_t &basis_sys_tester
64 ,basis_sys_tester_(basis_sys_tester)
65 ,D_imp_used_(MAT_IMP_AUTO)
71 const VectorSpace::space_ptr_t &space_x
72 ,
const VectorSpace::space_ptr_t &space_c
76 namespace rcp = MemMngPack;
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!" );
86 num_vars = basis_sys->var_dep().size() + basis_sys->var_indep().size();
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!" );
103 space_range_ = space_x_->sub_space(basis_sys->var_dep())->clone();
104 space_null_ = space_x_->sub_space(basis_sys->var_indep())->clone();
107 space_range_ = Teuchos::null;
108 space_null_ = Teuchos::null;
130 namespace rcp = MemMngPack;
133 if( out && olevel >= PRINT_BASIC_INFO )
134 *out <<
"\n****************************************************************"
135 <<
"\n*** DecompositionSystemVarReductImp::get_basis_matrices(...) ***"
136 <<
"\n****************************************************************\n";
168 bool new_D_mat_object =
true;
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;
179 const_cast<MatrixOp*
>(Z_vr->
D_ptr().
get()) );
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;
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;
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;
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;
219 if( new_D_mat_object ) {
221 alloc_new_D_matrix( out, olevel, &_D_ptr );
225 _D_ptr = Teuchos::rcp_const_cast<MatrixOp>(Z_vr->
D_ptr());
229 if(D_imp_used_ == MAT_IMP_IMPLICIT) {
232 (*D_ptr) = Teuchos::null;
235 D_ptr_ = Teuchos::null;
248 bool new_C_mat_object;
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;
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;
263 new_C_mat_object =
false;
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";
279 if( out && olevel >= PRINT_BASIC_INFO )
280 *out <<
"\nEnd DecompositionSystemVarReductImp::get_basis_matrices(...)\n";
297 if(basis_sys.
get()) {
299 space_range_ = space_x_->sub_space(basis_sys->var_dep())->clone();
300 space_null_ = space_x_->sub_space(basis_sys->var_indep())->clone();
309 return space_x_->dim();
317 return space_c_->dim();
324 if(basis_sys_.get()) {
325 return basis_sys_->equ_decomp().size();
332 const VectorSpace::space_ptr_t
338 const VectorSpace::space_ptr_t
349 namespace rcp = MemMngPack;
374 namespace rcp = MemMngPack;
377 if( out && olevel >= PRINT_BASIC_INFO ) {
378 *out <<
"\n***********************************************************"
379 <<
"\n*** DecompositionSystemVarReductImp::update_decomp(...) ***"
380 <<
"\n************************************************************\n";
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";
390 basis_sys_.get() == NULL, std::logic_error
391 ,
"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
392 "a BasisSystem object has not been set yet!" );
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!" );
448 out,olevel,test_what,Z,Y,R,Uz,Uy,&C_ptr,&D_ptr);
455 if( D_ptr.
get() == NULL ) {
458 if( D_ptr_.get() != NULL )
461 alloc_new_D_matrix( out, olevel, &D_ptr );
464 if( C_ptr_.get() == NULL ) {
471 if( out && olevel >= PRINT_BASIC_INFO )
472 *out <<
"\nUpdating the basis matrix C and other matices using the BasisSystem object ...\n";
477 basis_sys_->update_basis(
480 ,D_imp_used_ == MAT_IMP_EXPLICIT ? D_ptr.
get() : NULL
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
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!" );
514 Teuchos::rcp_const_cast<MatrixOp>(GcDd_ptr)
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
538 print_tests = BasisSystemTester::PRINT_NONE;
540 case PRINT_BASIC_INFO:
541 print_tests = BasisSystemTester::PRINT_BASIC;
543 case PRINT_MORE_INFO:
544 print_tests = BasisSystemTester::PRINT_MORE;
547 case PRINT_EVERY_THING:
548 print_tests = BasisSystemTester::PRINT_ALL;
551 basis_sys_tester_->print_tests(print_tests);
552 basis_sys_tester_->dump_all( olevel == PRINT_EVERY_THING );
554 const bool passed = basis_sys_tester_->test_basis_system(
559 ,D_imp_used_ == MAT_IMP_EXPLICIT ? D_ptr.
get() : NULL
565 ,
"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
566 "Test of the basis system failed!" );
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!" );
609 ,MatrixIdentConcatStd::TOP
619 C_ptr_ = Teuchos::null;
620 D_ptr_ = Teuchos::null;
622 if( out && olevel >= PRINT_BASIC_INFO )
623 *out <<
"\nEnd DecompositionSystemVarReductImp::update_decomp(...)\n";
627 std::ostream& out,
const std::string& L
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"
637 << L <<
" D = -inv(C)*N computed explicity (using basis_sys)\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"
645 << L <<
"end update of Y, R and Uy\n"
653 return basis_sys_.get() ? basis_sys_->var_indep() : Range1D::Invalid;
658 return basis_sys_.get() ? basis_sys_->var_dep() : Range1D::Invalid;
665 *D_imp_used = ( D_imp() == MAT_IMP_AUTO
672 void DecompositionSystemVarReductImp::alloc_new_D_matrix(
678 if(D_imp_used_ == MAT_IMP_IMPLICIT) {
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";
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";
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
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.
const mat_fcty_ptr_t factory_Uz() const
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...
Specialization of DecompositionSystem for variable reduction decompositions.
Teuchos::RCP< const Teuchos::AbstractFactory< MatrixOp > > mat_fcty_ptr_t
const mat_fcty_ptr_t factory_Z() const
Implements D = - inv(C) * N for a variable reduction projection.
Range1D var_indep() const
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...
size_type n() const
Returns this->space_x()->dim().
size_type m() const
Returns this->space_c()->dim().
#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.