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