ConstrainedOptPack_DecompositionSystemVarReductImp.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <typeinfo>
00030 
00031 #include "ConstrainedOptPack_DecompositionSystemVarReductImp.hpp"
00032 #include "ConstrainedOptPack_MatrixIdentConcatStd.hpp"
00033 #include "ConstrainedOptPack_MatrixVarReductImplicit.hpp"
00034 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00035 #include "AbstractLinAlgPack_MatrixOpSubView.hpp"
00036 #include "Teuchos_AbstractFactoryStd.hpp"
00037 #include "Teuchos_dyn_cast.hpp"
00038 #include "Teuchos_TestForException.hpp"
00039 
00040 namespace ConstrainedOptPack {
00041 
00042 DecompositionSystemVarReductImp::DecompositionSystemVarReductImp(
00043   const VectorSpace::space_ptr_t     &space_x
00044   ,const VectorSpace::space_ptr_t    &space_c
00045   ,const basis_sys_ptr_t             &basis_sys
00046   ,const basis_sys_tester_ptr_t      &basis_sys_tester
00047   ,EExplicitImplicit                 D_imp
00048   ,EExplicitImplicit                 Uz_imp
00049   )
00050   :DecompositionSystemVarReduct(D_imp,Uz_imp)
00051   ,basis_sys_tester_(basis_sys_tester)
00052   ,D_imp_used_(MAT_IMP_AUTO)
00053 {
00054   this->initialize(space_x,space_c,basis_sys);
00055 }
00056 
00057 void DecompositionSystemVarReductImp::initialize(
00058   const VectorSpace::space_ptr_t     &space_x
00059   ,const VectorSpace::space_ptr_t    &space_c
00060   ,const basis_sys_ptr_t             &basis_sys
00061   )
00062 {
00063   namespace rcp = MemMngPack;
00064   size_type num_vars = 0;
00065 #ifdef TEUCHOS_DEBUG
00066   TEST_FOR_EXCEPTION(
00067     basis_sys.get() != NULL && (space_x.get() == NULL || space_c.get() == NULL)
00068     ,std::invalid_argument
00069     ,"DecompositionSystemVarReductImp::initialize(...) : Error, "
00070     "if basis_sys is set, then space_x and space_c must also be set!" );
00071 #endif
00072   if(basis_sys.get()) {
00073     num_vars = basis_sys->var_dep().size() + basis_sys->var_indep().size();
00074 #ifdef TEUCHOS_DEBUG
00075     const size_type
00076       space_x_dim = space_x->dim(),
00077       space_c_dim = space_c->dim(),
00078       num_equ     = basis_sys->equ_decomp().size() + basis_sys->equ_undecomp().size();
00079     TEST_FOR_EXCEPTION(
00080       num_vars != 0 && (space_x_dim != num_vars || space_c_dim != num_equ)
00081       , std::invalid_argument
00082       ,"DecompositionSystemVarReductImp::initialize(...) : Error, "
00083       "the dimensions of space_x, space_c and basis_sys do not match up!" );
00084 #endif
00085   }
00086   space_x_    = space_x;
00087   space_c_    = space_c;
00088   basis_sys_  = basis_sys;
00089   if(num_vars) {
00090     space_range_ = space_x_->sub_space(basis_sys->var_dep())->clone();
00091     space_null_  = space_x_->sub_space(basis_sys->var_indep())->clone();
00092   }
00093   else {
00094     space_range_ = Teuchos::null;
00095     space_null_  = Teuchos::null;
00096   }
00097 }
00098 
00099 // Basis manipulation
00100 
00101 void DecompositionSystemVarReductImp::get_basis_matrices(
00102   std::ostream                                      *out
00103   ,EOutputLevel                                     olevel
00104   ,ERunTests                                        test_what
00105   ,MatrixOp                                         *Z
00106   ,MatrixOp                                         *Y
00107   ,MatrixOpNonsing                                  *R
00108   ,MatrixOp                                         *Uz
00109   ,MatrixOp                                         *Uy
00110   ,Teuchos::RCP<MatrixOpNonsing>       *C_ptr
00111   ,Teuchos::RCP<MatrixOp>              *D_ptr
00112   )
00113 {
00114 
00115   // ToDo: Implement undecomposed general equalities and general inequalities
00116 
00117   namespace rcp = MemMngPack;
00118   using Teuchos::dyn_cast;
00119 
00120   if( out && olevel >= PRINT_BASIC_INFO )
00121     *out << "\n****************************************************************"
00122        << "\n*** DecompositionSystemVarReductImp::get_basis_matrices(...) ***"
00123        << "\n****************************************************************\n";
00124 
00125   // ToDo: Validate input!
00126 
00127   //
00128   // Get references to concrete matrix objects
00129   //
00130 
00131   MatrixIdentConcatStd
00132     *Z_vr = Z ? &dyn_cast<MatrixIdentConcatStd>(*Z) : NULL;
00133 
00134   //
00135   // Make all matrices but Z uninitialized to determine if we can
00136   // reuse the Z.D matrix object (explicit or implicit).
00137   // Also, return a reference to the basis matrix C from the
00138   // R object.
00139   //
00140 
00141   (*C_ptr) = uninitialize_matrices(out,olevel,Y,R,Uy);
00142 
00143   //
00144   // Determine if we should be using an explicit or implicit D = -inv(C)*N object
00145   // (if we are allowed to choose).
00146   //
00147   update_D_imp_used(&D_imp_used_);
00148 
00149   //
00150   // Determine if we need to allocate a new matrix object for Z.D.
00151   // Also, get a reference to the explicit D matrix (if one exits)
00152   // and remove any reference to the basis matrix C by Z.D.
00153   //
00154 
00155   bool new_D_mat_object = true; // Valgrind complains if this is not initialized.
00156   if( Z_vr ) {
00157     if( Z_vr->D_ptr().get() == NULL ) {
00158       if( out && olevel >= PRINT_BASIC_INFO )
00159         *out << "\nMust allocate a new matrix object for D = -inv(C)*N since "
00160            << "one has not been allocated yet ...\n";
00161       new_D_mat_object = true;
00162     }
00163     else {
00164       MatrixVarReductImplicit
00165         *D_vr = dynamic_cast<MatrixVarReductImplicit*>(
00166           const_cast<MatrixOp*>(Z_vr->D_ptr().get()) );
00167       // We may have to reallocate a new object if someone is sharing it or
00168       // if we are switching from implicit to explicit or visa-versa.
00169       if( Z_vr->D_ptr().count() > 1 ) {
00170         if( out && olevel >= PRINT_BASIC_INFO )
00171           *out << "\nMust allocate a new matrix object for D = -inv(C)*N since someone "
00172              << "else is using the current one ...\n";
00173         new_D_mat_object = true;
00174       }
00175       else if( D_vr != NULL ) {
00176         if( D_imp_used_ == MAT_IMP_EXPLICIT ) {
00177           if( out && olevel >= PRINT_BASIC_INFO )
00178             *out << "\nMust allocate a new matrix object for D = -inv(C)*N since we "
00179                << "are switching from implicit to explicit ...\n";
00180           new_D_mat_object = true;
00181         }
00182       }
00183       else if( D_imp_used_ == MAT_IMP_IMPLICIT ) {
00184         if( out && olevel >= PRINT_BASIC_INFO )
00185           *out << "\nMust allocate a new matrix object for D = -inv(C)*N since we "
00186              << "are switching from explicit to implicit ...\n";
00187         new_D_mat_object = true;
00188       }
00189       // Remove the reference to the basis matrix C
00190       if(D_vr)
00191         D_vr->set_uninitialized();
00192     }
00193   }
00194   else {
00195     if( out && olevel >= PRINT_BASIC_INFO )
00196       *out << "\nMust allocate a new matrix object for D = -inv(C)*N since "
00197          << " Z == NULL was passed in ...\n";
00198     new_D_mat_object = true;
00199   }
00200 
00201   //
00202   // Get the matrix object of D and allocate a new matrix object if needed.
00203   //
00204 
00205   Teuchos::RCP<MatrixOp> _D_ptr;
00206   if( new_D_mat_object ) {
00207     // Create a new matrix object!
00208     alloc_new_D_matrix( out, olevel, &_D_ptr );
00209   }
00210   else {
00211     // Use current matrix object!
00212     _D_ptr = Teuchos::rcp_const_cast<MatrixOp>(Z_vr->D_ptr());
00213   }
00214 
00215   // Set cached implicit D matrix or external explicit D matrix
00216   if(D_imp_used_ == MAT_IMP_IMPLICIT) {
00217     D_ptr_ = _D_ptr;     // Need to remember this for when update_decomp() is called!
00218     if(D_ptr)
00219       (*D_ptr) = Teuchos::null;  // No explicit D matrix
00220   }
00221   else {
00222     D_ptr_ = Teuchos::null;  // This matrix will be passed back in set_basis_matrices(...) before update_decomp(...) is called.
00223     if(D_ptr)
00224       (*D_ptr) = _D_ptr;  // Externalize explicit D matrix.
00225   }
00226 
00227   //
00228   // Determine if we need to allocate a new basis matrix object C.
00229   //
00230   // At this point, if a matrix object for C already exits and noone
00231   // outside has a reference to this basis matrix object then the only
00232   // references to it should be in C_ptr
00233   //
00234 
00235   bool new_C_mat_object; // compiler should warn if used before initialized!
00236   if( (*C_ptr).get() == NULL ) {
00237     if( out && olevel >= PRINT_BASIC_INFO )
00238       *out << "\nMust allocate a new basis matrix object for C since "
00239          << "one has not been allocated yet ...\n";
00240     new_C_mat_object = true;
00241   }
00242   else {
00243     if( (*C_ptr).count() > 1 ) {
00244       if( out && olevel >= PRINT_BASIC_INFO )
00245         *out << "\nMust allocate a new basis matrix object C since someone "
00246            << "else is using the current one ...\n";
00247       new_C_mat_object = true;
00248     }
00249     else {
00250       new_C_mat_object = false; // The current C matrix is owned only by us!
00251     }
00252   }
00253   
00254   //
00255   // Get the basis matrix object C and allocate a new object for if we have to.
00256   //
00257 
00258   if( new_C_mat_object) {
00259     (*C_ptr) = basis_sys_->factory_C()->create();
00260     if( out && olevel >= PRINT_BASIC_INFO )
00261       *out << "\nAllocated a new basis matrix object C "
00262          << "of type \'" << typeName(*(*C_ptr)) << "\' ...\n";
00263   }
00264 
00265 
00266   if( out && olevel >= PRINT_BASIC_INFO )
00267     *out << "\nEnd DecompositionSystemVarReductImp::get_basis_matrices(...)\n";
00268 
00269 }
00270 
00271 void DecompositionSystemVarReductImp::set_basis_matrices(
00272   std::ostream                                           *out
00273   ,EOutputLevel                                          olevel
00274   ,ERunTests                                             test_what
00275   ,const Teuchos::RCP<MatrixOpNonsing>      &C_ptr
00276   ,const Teuchos::RCP<MatrixOp>             &D_ptr
00277   ,MatrixOp                                              *Uz
00278   ,const basis_sys_ptr_t                                 &basis_sys
00279   )
00280 {
00281   C_ptr_ = C_ptr;
00282   if( D_ptr.get() )
00283     D_ptr_ = D_ptr;
00284   if(basis_sys.get()) {
00285     basis_sys_   = basis_sys;
00286     space_range_ = space_x_->sub_space(basis_sys->var_dep())->clone();
00287     space_null_  = space_x_->sub_space(basis_sys->var_indep())->clone();
00288   }
00289 }
00290 
00291 // Overridden from DecompositionSystem
00292 
00293 size_type DecompositionSystemVarReductImp::n() const
00294 {
00295   if(space_x_.get()) {
00296     return space_x_->dim();
00297   }
00298   return 0; // Not fully initialized!
00299 }
00300 
00301 size_type DecompositionSystemVarReductImp::m() const
00302 {
00303   if(space_c_.get()) {
00304     return space_c_->dim();
00305   }
00306   return 0; // Not fully initialized!
00307 }
00308 
00309 size_type DecompositionSystemVarReductImp::r() const
00310 {
00311   if(basis_sys_.get()) {
00312     return basis_sys_->equ_decomp().size();
00313   }
00314   return 0; // Not fully initialized!
00315 }
00316 
00317 // range and null spaces
00318 
00319 const VectorSpace::space_ptr_t
00320 DecompositionSystemVarReductImp::space_range() const
00321 {
00322   return space_range_;
00323 }
00324 
00325 const VectorSpace::space_ptr_t
00326 DecompositionSystemVarReductImp::space_null() const
00327 {
00328   return space_null_;
00329 }
00330 
00331 // matrix factories
00332 
00333 const DecompositionSystem::mat_fcty_ptr_t
00334 DecompositionSystemVarReductImp::factory_Z() const
00335 {
00336   namespace rcp = MemMngPack;
00337   return Teuchos::rcp(
00338     new Teuchos::AbstractFactoryStd<MatrixOp,MatrixIdentConcatStd>()
00339     );
00340 }
00341 
00342 const DecompositionSystem::mat_fcty_ptr_t
00343 DecompositionSystemVarReductImp::factory_Uz() const
00344 {
00345   return Teuchos::rcp(  new Teuchos::AbstractFactoryStd<MatrixOp,MatrixOpSubView>() );
00346 }
00347 
00348 void DecompositionSystemVarReductImp::update_decomp(
00349   std::ostream          *out
00350   ,EOutputLevel         olevel
00351   ,ERunTests            test_what
00352   ,const MatrixOp       &Gc
00353   ,MatrixOp             *Z
00354   ,MatrixOp             *Y
00355   ,MatrixOpNonsing      *R
00356   ,MatrixOp             *Uz
00357   ,MatrixOp             *Uy
00358   ,EMatRelations        mat_rel
00359   ) const
00360 {
00361   namespace rcp = MemMngPack;
00362   using Teuchos::dyn_cast;
00363 
00364   if( out && olevel >= PRINT_BASIC_INFO ) {
00365     *out << "\n***********************************************************"
00366        << "\n*** DecompositionSystemVarReductImp::update_decomp(...) ***"
00367        << "\n************************************************************\n";
00368 
00369     if(mat_rel != MATRICES_INDEP_IMPS)
00370       *out << "\nWarning!!! mat_rel != MATRICES_INDEP_IMPS; The decompsition matrix "
00371          << "objects may not be independent of each other!\n"; 
00372   }
00373 
00374 #ifdef TEUCHOS_DEBUG
00375   // Validate setup
00376   TEST_FOR_EXCEPTION(
00377     basis_sys_.get() == NULL, std::logic_error
00378     ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00379     "a BasisSystem object has not been set yet!" );
00380 #endif
00381   
00382   const size_type
00383     n = this->n(),
00384     m = this->m(),
00385     r = this->r();
00386   const Range1D
00387     var_indep(r+1,n),
00388     equ_decomp   = this->equ_decomp();
00389 
00390 #ifdef TEUCHOS_DEBUG
00391   // Validate input
00392   TEST_FOR_EXCEPTION(
00393       Z==NULL&&Y==NULL&&R==NULL&&Uz==NULL&&Uy==NULL
00394     , std::invalid_argument
00395     ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00396     "at least one of Z, Y, R, Uz and Uycan not be NULL!" );
00397   TEST_FOR_EXCEPTION(
00398     m == r && Uz != NULL, std::invalid_argument
00399     ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00400     "Uz must be NULL if m==r is NULL!" );
00401   TEST_FOR_EXCEPTION(
00402     m == r && Uy != NULL, std::invalid_argument
00403     ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00404     "Uy must be NULL if m==r is NULL!" );
00405 #endif
00406 
00407   //
00408   // Get references to concrete matrix objects
00409   //
00410 
00411   MatrixIdentConcatStd
00412     *Z_vr = Z ? &dyn_cast<MatrixIdentConcatStd>(*Z) : NULL;
00413   TEST_FOR_EXCEPT( !( Uz == NULL ) ); // ToDo: Implement undecomposed general equalities
00414 
00415   //
00416   // Get smart pointers to unreferenced C and D matrix objects.
00417   //
00418 
00419   Teuchos::RCP<MatrixOpNonsing>    C_ptr;
00420   Teuchos::RCP<MatrixOp>               D_ptr;
00421 
00422   if( C_ptr_.get() ) {
00423     //
00424     // The function get_basis_matrices() was called by external client
00425     // so we don't need to call it here or update the decomposition matrices.
00426     //
00427     C_ptr = C_ptr_; // This was set by set_basis_matrices()
00428   }
00429   else {
00430     //
00431     // Make all matrices uninitialized and get unreferenced smart
00432     // pointers to C and D (explicit only!).
00433     //
00434     const_cast<DecompositionSystemVarReductImp*>(this)->get_basis_matrices(
00435       out,olevel,test_what,Z,Y,R,Uz,Uy,&C_ptr,&D_ptr);
00436   }
00437 
00438   // Get the D matrix created by get_basis_matrices() and set by
00439   // get_basis_matrices() if implicit or set by set_basis_matrices()
00440   // if explicit.  This matrix may not be allocated yet in which
00441   // case we need to allocate it.
00442   if( D_ptr.get() == NULL ) {
00443     // D_ptr was not set in get_basis_matrix() in code above but
00444     // it may be cashed (if explicit) in D_ptr.
00445     if( D_ptr_.get() != NULL )
00446       D_ptr = D_ptr_;
00447     else
00448       alloc_new_D_matrix( out, olevel, &D_ptr );
00449   }
00450 
00451   if( C_ptr_.get() == NULL ) {
00452 
00453     //
00454     // The basis matrices were not updated by an external client
00455     // so we must do it ourselves here using basis_sys.
00456     //
00457   
00458     if( out && olevel >= PRINT_BASIC_INFO )
00459       *out << "\nUpdating the basis matrix C and other matices using the BasisSystem object ...\n";
00460   
00461     TEST_FOR_EXCEPT( !(  D_ptr.get()  ) ); // local programming error only!
00462     TEST_FOR_EXCEPT( !(  C_ptr.get()  ) ); // local programming error only!
00463   
00464     basis_sys_->update_basis(
00465       Gc                                                       // Gc
00466       ,C_ptr.get()                                             // C
00467       ,D_imp_used_ == MAT_IMP_EXPLICIT ? D_ptr.get() : NULL    // D?
00468       ,NULL                                                    // GcUP == Uz
00469       ,(mat_rel == mat_rel == MATRICES_INDEP_IMPS
00470         ? BasisSystem::MATRICES_INDEP_IMPS
00471         : BasisSystem::MATRICES_ALLOW_DEP_IMPS )
00472       ,out && olevel >= PRINT_BASIC_INFO ? out : NULL
00473       );
00474   }
00475     
00476   //
00477   // Create the matrix object: N = Gc(var_indep,cond_decomp)' 
00478   //
00479   Teuchos::RCP<const MatrixOp>
00480     N_ptr = Teuchos::null;
00481   if( D_imp_used_ == MAT_IMP_IMPLICIT ) {
00482     Teuchos::RCP<const MatrixOp>
00483       GcDd_ptr = Gc.sub_view(var_indep,equ_decomp);
00484     TEST_FOR_EXCEPTION(
00485       GcDd_ptr.get() == NULL, std::logic_error
00486       ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00487       "The matrix class used for the gradient of constraints matrix Gc of type \'"
00488       << typeName(Gc) << "\' must return return.get() != NULL from "
00489       "Gc.sub_view(var_indep,equ_decomp)!" );
00490     if(mat_rel == MATRICES_INDEP_IMPS) {
00491       GcDd_ptr = GcDd_ptr->clone();
00492       TEST_FOR_EXCEPTION(
00493         GcDd_ptr.get() == NULL, std::logic_error
00494         ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00495         "The matrix class used for the gradient of constraints matrix Gc.sub_view(var_indep,equ_decomp) "
00496         "of type \'" << typeName(*GcDd_ptr) << "\' must return return.get() != NULL from \n"
00497         "Gc.sub_view(var_indep,equ_decomp)->clone() since mat_rel == MATRICES_INDEP_IMPS!" );
00498     }
00499     N_ptr = Teuchos::rcp(
00500       new MatrixOpSubView(
00501         Teuchos::rcp_const_cast<MatrixOp>(GcDd_ptr)  // M_full
00502         ,Range1D()                                   // rng_rows
00503         ,Range1D()                                   // rng_cols
00504         ,BLAS_Cpp::trans                             // M_trans
00505         )
00506       );
00507   }
00508 
00509   //
00510   // Test the basis matrix C and the other matrices.
00511   //
00512   
00513   if( test_what == RUN_TESTS ) {
00514     if( out && olevel >= PRINT_BASIC_INFO )
00515       *out << "\nTesting the basis matrix C and other matices updated using the BasisSystem object ...\n";
00516     TEST_FOR_EXCEPTION(
00517       basis_sys_tester_.get() == NULL, std::logic_error
00518       ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00519       "test_what == RUN_TESTS but this->basis_sys_tester().get() == NULL!" );
00520     if( basis_sys_tester_->print_tests() == BasisSystemTester::PRINT_NOT_SELECTED ) {
00521       BasisSystemTester::EPrintTestLevel
00522         print_tests;
00523       switch(olevel) {
00524         case PRINT_NONE:
00525           print_tests = BasisSystemTester::PRINT_NONE;
00526           break;
00527         case PRINT_BASIC_INFO:
00528           print_tests = BasisSystemTester::PRINT_BASIC;
00529           break;
00530         case PRINT_MORE_INFO:
00531           print_tests = BasisSystemTester::PRINT_MORE;
00532           break;
00533         case PRINT_VECTORS:
00534         case PRINT_EVERY_THING:
00535           print_tests = BasisSystemTester::PRINT_ALL;
00536           break;
00537       }
00538       basis_sys_tester_->print_tests(print_tests);
00539       basis_sys_tester_->dump_all( olevel == PRINT_EVERY_THING );
00540     }
00541     const bool passed = basis_sys_tester_->test_basis_system(
00542       *basis_sys_                                              // basis_sys
00543       ,&Gc                                                     // Gc
00544       ,C_ptr.get()                                             // C
00545       ,N_ptr.get()                                             // N
00546       ,D_imp_used_ == MAT_IMP_EXPLICIT ? D_ptr.get() : NULL    // D?
00547       ,NULL                                                    // GcUP == Uz
00548       ,out                                                     // out
00549       );
00550     TEST_FOR_EXCEPTION(
00551       !passed, TestFailed
00552       ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00553       "Test of the basis system failed!" );
00554   }
00555 
00556   //
00557   // Initialize the implicit D = -inv(C)*N matrix object.
00558   //
00559 
00560   TEST_FOR_EXCEPT( !( D_ptr.get() ) ); // local programming error only?
00561   if( D_imp_used_ == MAT_IMP_IMPLICIT ) {
00562     if( !C_ptr.has_ownership() && mat_rel == MATRICES_INDEP_IMPS ) {
00563       C_ptr = C_ptr->clone_mwons();
00564       TEST_FOR_EXCEPTION(
00565         C_ptr.get() == NULL, std::logic_error
00566         ,"DecompositionSystemVarReductImp::update_decomp(...) : Error, "
00567         "The matrix class used for the basis matrix C (from the BasisSystem object) "
00568         "must return return.get() != NULL from the clone_mwons() method since "
00569         "mat_rel == MATRICES_INDEP_IMPS!" );
00570     }
00571     dyn_cast<MatrixVarReductImplicit>(*D_ptr).initialize(
00572       C_ptr       // C
00573       ,N_ptr      // N
00574       ,Teuchos::null  // D_direct
00575       );
00576     // Above, if the basis matrix object C is not owned by the smart
00577     // reference counted pointer object C_ptr, then we need to make the
00578     // reference that the implicit D = -inv(C)*N matrix object has to
00579     // C independent and the clone() method does that very nicely.
00580   }
00581 
00582   //
00583   // Call on the subclass to construct Y, R and Uy given C and D.
00584   //
00585 
00586   initialize_matrices(out,olevel,C_ptr,D_ptr,Y,R,Uy,mat_rel);
00587 
00588   //
00589   // Reconstruct Z, Uz and Uy.
00590   //
00591 
00592   if( Z_vr ) {
00593     Z_vr->initialize(
00594       space_x_                                   // space_cols
00595       ,space_x_->sub_space(var_indep)->clone()   // space_rows
00596       ,MatrixIdentConcatStd::TOP                 // top_or_bottom
00597       ,1.0                                       // alpha
00598       ,D_ptr                                     // D_ptr
00599       ,BLAS_Cpp::no_trans                        // D_trans
00600       );
00601   }
00602 
00603   TEST_FOR_EXCEPT( !( Uz == NULL ) ); // ToDo: Implement for undecomposed general equalities
00604 
00605   // Clear cache for basis matrices.
00606   C_ptr_ = Teuchos::null;
00607   D_ptr_ = Teuchos::null;
00608 
00609   if( out && olevel >= PRINT_BASIC_INFO )
00610     *out << "\nEnd DecompositionSystemVarReductImp::update_decomp(...)\n";
00611 }
00612 
00613 void DecompositionSystemVarReductImp::print_update_decomp(
00614   std::ostream& out, const std::string& L
00615   ) const
00616 {
00617   out
00618     << L << "*** Variable reduction decomposition (class DecompositionSytemVarReductImp)\n"
00619     << L << "C = Gc(var_dep,equ_decomp)' (using basis_sys)\n"
00620     << L << "if C is nearly singular then throw SingularDecomposition exception\n"
00621     << L << "if D_imp == MAT_IMP_IMPICIT then\n"
00622     << L << "  D = -inv(C)*N represented implicitly (class MatrixVarReductImplicit)\n"
00623     << L << "else\n"
00624     << L << "  D = -inv(C)*N computed explicity (using basis_sys)\n"
00625     << L << "end\n"
00626     << L << "Z = [ D; I ] (class MatrixIdentConcatStd)\n"
00627     << L << "Uz = Gc(var_indep,equ_undecomp)' - Gc(var_dep,equ_undecomp)'*D\n"
00628     << L << "begin update Y, R and Uy\n"
00629     ;
00630   print_update_matrices( out, L + "  " );
00631   out
00632     << L << "end update of Y, R and Uy\n"
00633     ;
00634 }
00635 
00636 // Overridden from DecompositionSystemVarReduct
00637 
00638 Range1D DecompositionSystemVarReductImp::var_indep() const
00639 {
00640   return basis_sys_.get() ? basis_sys_->var_indep() : Range1D::Invalid;
00641 }
00642 
00643 Range1D DecompositionSystemVarReductImp::var_dep() const
00644 {
00645   return basis_sys_.get() ? basis_sys_->var_dep() : Range1D::Invalid;
00646 }
00647 
00648 // protected
00649 
00650 void DecompositionSystemVarReductImp::update_D_imp_used(EExplicitImplicit *D_imp_used) const
00651 {
00652   *D_imp_used = ( D_imp() == MAT_IMP_AUTO
00653            ? MAT_IMP_IMPLICIT     // Without better info, use implicit by default!
00654            : D_imp() );
00655 }
00656 
00657 // private
00658 
00659 void DecompositionSystemVarReductImp::alloc_new_D_matrix( 
00660   std::ostream                             *out
00661   ,EOutputLevel                            olevel
00662   ,Teuchos::RCP<MatrixOp> *D_ptr
00663   ) const
00664 {
00665   if(D_imp_used_ == MAT_IMP_IMPLICIT) {
00666     (*D_ptr) = Teuchos::rcp(new MatrixVarReductImplicit());
00667     if( out && olevel >= PRINT_BASIC_INFO )
00668       *out << "\nAllocated a new implicit matrix object for D = -inv(C)*N "
00669          << "of type \'MatrixVarReductImplicit\' ...\n";
00670   }
00671   else {
00672     (*D_ptr) = basis_sys_->factory_D()->create();
00673     if( out && olevel >= PRINT_BASIC_INFO )
00674       *out << "\nAllocated a new explicit matrix object for D = -inv(C)*N "
00675          << "of type \'" << typeName(*(*D_ptr)) << "\' ...\n";
00676   }
00677 }
00678 
00679 } // end namespace ConstrainedOptPack

Generated on Tue Oct 20 12:51:44 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7