MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_BasisSystemPermDirectSparse.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 // 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 "AbstractLinAlgPack_BasisSystemPermDirectSparse.hpp"
00043 #include "AbstractLinAlgPack_MatrixOp.hpp"
00044 #include "AbstractLinAlgPack_MatrixOpNonsingAggr.hpp"
00045 #include "AbstractLinAlgPack_MatrixPermAggr.hpp"
00046 #include "AbstractLinAlgPack_MatrixConvertToSparseEncap.hpp"
00047 #include "AbstractLinAlgPack_MatrixExtractSparseElements.hpp"
00048 #include "AbstractLinAlgPack_MultiVectorMutableDense.hpp"
00049 #include "AbstractLinAlgPack_PermutationSerial.hpp"
00050 #include "AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp"
00051 #include "Teuchos_AbstractFactoryStd.hpp"
00052 #include "Teuchos_Assert.hpp"
00053 #include "Teuchos_dyn_cast.hpp"
00054 
00055 namespace AbstractLinAlgPack {
00056 
00057 BasisSystemPermDirectSparse::BasisSystemPermDirectSparse(
00058   const direct_solver_ptr_t&   direct_solver
00059   )
00060   :BasisSystemPerm(
00061     Teuchos::rcp(
00062       new Teuchos::AbstractFactoryStd<MatrixSymOp,MatrixSymPosDefCholFactor>())         // D'*D
00063     ,Teuchos::rcp(
00064       new Teuchos::AbstractFactoryStd<MatrixSymOpNonsing,MatrixSymPosDefCholFactor>())  // S
00065     )
00066 {
00067   this->initialize(direct_solver);
00068 }
00069   
00070 void BasisSystemPermDirectSparse::initialize(
00071   const direct_solver_ptr_t&   direct_solver
00072   )
00073 {
00074   direct_solver_ = direct_solver;
00075   n_ = m_ = r_ = Gc_nz_ = 0;
00076   init_var_inv_perm_.resize(0);
00077   init_equ_inv_perm_.resize(0);
00078   init_var_rng_    = Range1D::Invalid;
00079     init_equ_rng_    = Range1D::Invalid;
00080   var_dep_         = Range1D::Invalid;
00081   var_indep_       = Range1D::Invalid;
00082     equ_decomp_      = Range1D::Invalid;
00083     equ_undecomp_    = Range1D::Invalid;
00084 }
00085 
00086 // Overridden from BasisSystem
00087 
00088 const BasisSystem::mat_nonsing_fcty_ptr_t
00089 BasisSystemPermDirectSparse::factory_C() const
00090 {
00091   return Teuchos::rcp(
00092     new Teuchos::AbstractFactoryStd<MatrixOpNonsing,MatrixOpNonsingAggr>()
00093     );
00094 }
00095 
00096 const BasisSystem::mat_fcty_ptr_t
00097 BasisSystemPermDirectSparse::factory_D() const
00098 {
00099   return Teuchos::rcp(
00100     new Teuchos::AbstractFactoryStd<MatrixOp,MultiVectorMutableDense>()
00101     );
00102 }
00103 
00104 const BasisSystem::mat_fcty_ptr_t
00105 BasisSystemPermDirectSparse::factory_GcUP() const
00106 {
00107   return Teuchos::rcp(
00108     new Teuchos::AbstractFactoryStd<MatrixOp,MultiVectorMutableDense>()
00109     );
00110 }
00111 
00112 Range1D BasisSystemPermDirectSparse::var_dep() const
00113 {
00114   return var_dep_;
00115 }
00116 
00117 Range1D BasisSystemPermDirectSparse::var_indep() const
00118 {
00119   return var_indep_;
00120 }
00121 
00122 Range1D BasisSystemPermDirectSparse::equ_decomp() const
00123 {
00124   return equ_decomp_;
00125 }
00126 
00127 Range1D BasisSystemPermDirectSparse::equ_undecomp() const
00128 {
00129   return equ_undecomp_;
00130 }
00131 
00132 void BasisSystemPermDirectSparse::update_basis(
00133   const MatrixOp          &Gc
00134   ,MatrixOpNonsing        *C
00135   ,MatrixOp               *D
00136   ,MatrixOp               *GcUP
00137   ,EMatRelations          mat_rel
00138   ,std::ostream           *out
00139   ) const
00140 {
00141   using Teuchos::dyn_cast;
00142   if(out)
00143     *out << "\nUsing a direct sparse solver to update basis ...\n";
00144   const size_type
00145     n  = Gc.rows(),
00146     m  = Gc.cols();
00147 #ifdef TEUCHOS_DEBUG
00148   const size_type Gc_rows = n, Gc_cols = m, Gc_nz = Gc.nz();
00149   TEUCHOS_TEST_FOR_EXCEPTION(
00150     Gc_rows != n_ || Gc_cols != m_ || Gc_nz != Gc_nz_, std::invalid_argument
00151     ,"BasisSystemPermDirectSparse::set_basis(...) : Error, "
00152     "This matrix object is not compatible with last call to set_basis() or select_basis()!" );
00153   TEUCHOS_TEST_FOR_EXCEPTION(
00154     C == NULL, std::invalid_argument
00155     ,"BasisSystemPermDirectSparse::set_basis(...) : Error!" );
00156 #endif
00157   // Get the aggregate matrix object for Gc
00158   const MatrixPermAggr  
00159     &Gc_pa = dyn_cast<const MatrixPermAggr>(Gc);
00160   // Get the basis matrix object from the aggregate or allocate one
00161   MatrixOpNonsingAggr
00162     &C_aggr = dyn_cast<MatrixOpNonsingAggr>(*C);
00163   Teuchos::RCP<DirectSparseSolver::BasisMatrix>
00164     C_bm = get_basis_matrix(C_aggr);
00165   // Setup the encapulated convert-to-sparse matrix object
00166   MatrixConvertToSparseEncap A_mctse;
00167   set_A_mctse( n, m, Gc_pa, &A_mctse );
00168   // Refactor this basis (it had better be full rank)!
00169   try {
00170     direct_solver_->factor(
00171       A_mctse
00172       ,C_bm.get()
00173       ,Teuchos::null // Same factorization structure as before
00174       ,out
00175       );
00176   }
00177   catch(const DirectSparseSolver::FactorizationFailure& excpt) {
00178     if(out)
00179       *out << "\nCurrent basis is singular : " << excpt.what() << std::endl
00180          << "Throwing SingularBasis exception to client ...\n";
00181     TEUCHOS_TEST_FOR_EXCEPTION(
00182       true, SingularBasis
00183       ,"BasisSystemPermDirectSparse::update_basis(...) : Error, the current basis "
00184       "is singular : " << excpt.what() );
00185   }
00186   // Update the aggregate basis matrix and compute the auxiliary projected matrices
00187   update_basis_and_auxiliary_matrices( Gc, C_bm, &C_aggr, D, GcUP );
00188 }
00189 
00190 // Overridded from BasisSystemPerm
00191 
00192 const AbstractLinAlgPack::BasisSystemPerm::perm_fcty_ptr_t
00193 BasisSystemPermDirectSparse::factory_P_var() const
00194 {
00195   TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement using PermutationSerial
00196   return Teuchos::null;
00197 }
00198 
00199 const AbstractLinAlgPack::BasisSystemPerm::perm_fcty_ptr_t
00200 BasisSystemPermDirectSparse::factory_P_equ() const
00201 {
00202   TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement using PermutationSerial
00203   return Teuchos::null;
00204 }
00205 
00206 const AbstractLinAlgPack::BasisSystemPerm::perm_fcty_ptr_t
00207 BasisSystemPermDirectSparse::factory_P_inequ() const
00208 {
00209   TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement using PermutationSerial
00210   return Teuchos::null;
00211 }
00212 
00213 void BasisSystemPermDirectSparse::set_basis(
00214   const Permutation          &P_var
00215   ,const Range1D             &var_dep
00216   ,const Permutation         *P_equ
00217   ,const Range1D             *equ_decomp
00218   ,const MatrixOp            &Gc
00219   ,MatrixOpNonsing           *C
00220   ,MatrixOp                  *D
00221   ,MatrixOp                  *GcUP
00222   ,EMatRelations             mat_rel
00223   ,std::ostream              *out
00224   )
00225 {
00226   using Teuchos::dyn_cast;
00227   if(out)
00228     *out << "\nUsing a direct sparse solver to set a new basis ...\n";
00229   const size_type
00230     n  = Gc.rows(),
00231     m  = Gc.cols();
00232 #ifdef TEUCHOS_DEBUG
00233   const size_type Gc_rows = n, Gc_cols = m, Gc_nz = Gc.nz();
00234   TEUCHOS_TEST_FOR_EXCEPTION(
00235     P_equ == NULL || equ_decomp == NULL, std::invalid_argument
00236     ,"BasisSystemPermDirectSparse::set_basis(...) : Error!" );
00237   TEUCHOS_TEST_FOR_EXCEPTION(
00238     C == NULL, std::invalid_argument
00239     ,"BasisSystemPermDirectSparse::set_basis(...) : Error!" );
00240 #endif
00241   // Get the aggreate matrix object for Gc
00242   const MatrixPermAggr  
00243     &Gc_pa = dyn_cast<const MatrixPermAggr>(Gc);
00244   // Get the basis matrix object from the aggregate or allocate one
00245   MatrixOpNonsingAggr
00246     &C_aggr = dyn_cast<MatrixOpNonsingAggr>(*C);
00247   Teuchos::RCP<DirectSparseSolver::BasisMatrix>
00248     C_bm = get_basis_matrix(C_aggr);
00249   // Get at the concreate permutation vectors
00250   const PermutationSerial
00251     &P_var_s = dyn_cast<const PermutationSerial>(P_var),
00252     &P_equ_s = dyn_cast<const PermutationSerial>(*P_equ);
00253   // Setup the encapulated convert-to-sparse matrix object
00254   init_var_inv_perm_  = *P_var_s.inv_perm();
00255   init_var_rng_       = var_dep;
00256   init_equ_inv_perm_  = *P_equ_s.inv_perm();
00257     init_equ_rng_       = *equ_decomp;
00258   MatrixConvertToSparseEncap A_mctse;
00259   set_A_mctse( n, m, Gc_pa, &A_mctse );
00260   // Analyze and factor this basis (it had better be full rank)!
00261   IVector row_perm_ds, col_perm_ds; // Must store these even though we don't want them!
00262   size_type rank = 0;
00263   direct_solver_->analyze_and_factor(
00264     A_mctse
00265     ,&row_perm_ds
00266     ,&col_perm_ds
00267     ,&rank
00268     ,C_bm.get()
00269     ,out
00270     );
00271   if( rank < var_dep.size() ) {
00272     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Throw an exception with a good error message!
00273   }
00274   // Update the rest of the basis stuff
00275   do_some_basis_stuff(Gc,var_dep,*equ_decomp,C_bm,&C_aggr,D,GcUP);
00276 }
00277 
00278 void BasisSystemPermDirectSparse::select_basis(
00279   const Vector               *nu
00280   ,MatrixOp                  *Gc
00281   ,Permutation               *P_var
00282   ,Range1D                   *var_dep
00283   ,Permutation               *P_equ
00284   ,Range1D                   *equ_decomp
00285   ,MatrixOpNonsing           *C
00286   ,MatrixOp                  *D
00287   ,MatrixOp                  *GcUP
00288   ,EMatRelations             mat_rel
00289   ,std::ostream              *out
00290   )
00291 {
00292   using Teuchos::dyn_cast;
00293   if(out)
00294     *out << "\nUsing a direct sparse solver to select a new basis ...\n";
00295 #ifdef TEUCHOS_DEBUG
00296   // Validate input
00297   const char msg_err_head[] = "BasisSystemPermDirectSparse::set_basis(...) : Error!";
00298   TEUCHOS_TEST_FOR_EXCEPTION(
00299     Gc == NULL, std::invalid_argument
00300     ,msg_err_head << " Must have equality constriants in this current implementation! " );
00301 #endif
00302   const size_type
00303     n  = Gc->rows(),
00304     m  = Gc->cols();
00305 #ifdef TEUCHOS_DEBUG
00306   // Validate input
00307   const size_type Gc_rows = Gc->rows(), Gc_cols = Gc->cols(), Gc_nz = Gc->nz();
00308   TEUCHOS_TEST_FOR_EXCEPTION(
00309     P_var == NULL || var_dep == NULL, std::invalid_argument, msg_err_head );
00310   TEUCHOS_TEST_FOR_EXCEPTION(
00311     P_equ == NULL || equ_decomp == NULL, std::invalid_argument, msg_err_head );
00312   TEUCHOS_TEST_FOR_EXCEPTION(
00313     C == NULL, std::invalid_argument, msg_err_head );
00314 #endif
00315   // Get the aggreate matrix object for Gc
00316   MatrixPermAggr  
00317     &Gc_pa = dyn_cast<MatrixPermAggr>(*Gc);
00318   // Get the basis matrix object from the aggregate or allocate one
00319   MatrixOpNonsingAggr
00320     &C_aggr = dyn_cast<MatrixOpNonsingAggr>(*C);
00321   Teuchos::RCP<DirectSparseSolver::BasisMatrix>
00322     C_bm = get_basis_matrix(C_aggr);
00323   // Setup the encapulated convert-to-sparse matrix object
00324   // ToDo: Use nu to exclude variables that are at a bound!
00325   init_var_rng_       = Range1D(1,n);
00326   init_var_inv_perm_.resize(0); // Not used since above is full range
00327   init_equ_rng_       = Range1D(1,m);
00328   init_equ_inv_perm_.resize(0); // Not used since above is full range
00329   MatrixConvertToSparseEncap A_mctse;
00330   set_A_mctse( n, m, Gc_pa, &A_mctse );
00331   // Analyze and factor this basis (it had better be full rank)!
00332   Teuchos::RCP<IVector>
00333     var_perm_ds = Teuchos::rcp(new IVector),
00334     equ_perm_ds = Teuchos::rcp(new IVector);
00335   size_type rank = 0;
00336   direct_solver_->analyze_and_factor(
00337     A_mctse
00338     ,equ_perm_ds.get()
00339     ,var_perm_ds.get()
00340     ,&rank
00341     ,C_bm.get()
00342     ,out
00343     );
00344   if( rank == 0 ) {
00345     TEUCHOS_TEST_FOR_EXCEPT( !(  rank == 0  ) ); // ToDo: Throw exception with good error message!
00346   }
00347   // Return the selected basis
00348   // ToDo: Use var_perm_ds and equ_perm_ds together with nu to
00349   // get the overall permuations for all of the variables.
00350   PermutationSerial
00351     &P_var_s = dyn_cast<PermutationSerial>(*P_var),
00352     &P_equ_s = dyn_cast<PermutationSerial>(*P_equ);
00353   // Create the overall permutations to set to the permutation matrices!
00354   *var_dep = Range1D(1,rank);
00355   P_var_s.initialize( var_perm_ds, Teuchos::null ); 
00356   *equ_decomp = Range1D(1,rank);
00357   P_equ_s.initialize( equ_perm_ds, Teuchos::null );
00358   // Setup Gc_aggr with Gc_perm
00359   const int          num_row_part = 2;
00360   const int          num_col_part = 2;
00361   const index_type   row_part[3]  = { 1, rank, n+1 };
00362   const index_type   col_part[3]  = { 1, rank, m+1 };
00363   Gc_pa.initialize(
00364     Gc_pa.mat_orig()
00365     ,Teuchos::rcp(new PermutationSerial(var_perm_ds,Teuchos::null)) // var_perm_ds reuse is okay!
00366     ,Teuchos::rcp(new PermutationSerial(equ_perm_ds,Teuchos::null)) // equ_perm_ds resue is okay!
00367     ,Gc_pa.mat_orig()->perm_view(
00368       P_var,row_part,num_row_part
00369       ,P_equ,col_part,num_col_part
00370       )
00371     );
00372   // Update the rest of the basis stuff
00373   do_some_basis_stuff(*Gc,*var_dep,*equ_decomp,C_bm,&C_aggr,D,GcUP);
00374 }
00375 
00376 // private
00377 
00378 Teuchos::RCP<DirectSparseSolver::BasisMatrix>
00379 BasisSystemPermDirectSparse::get_basis_matrix( MatrixOpNonsingAggr &C_aggr ) const
00380 {
00381   using Teuchos::dyn_cast;
00382   Teuchos::RCP<DirectSparseSolver::BasisMatrix> C_bm;
00383   if( C_aggr.mns().get() ) {
00384     C_bm = Teuchos::rcp_dynamic_cast<DirectSparseSolver::BasisMatrix>(
00385       Teuchos::rcp_const_cast<MatrixNonsing>(C_aggr.mns() ) );
00386     if(C_bm.get() == NULL)
00387       dyn_cast<const DirectSparseSolver::BasisMatrix>(*C_aggr.mns()); // Throws exception!
00388   }
00389   else {
00390     C_bm = direct_solver_->basis_matrix_factory()->create();
00391   }
00392   return C_bm;
00393 }
00394 
00395 void BasisSystemPermDirectSparse::set_A_mctse(
00396   size_type                    n
00397   ,size_type                   m
00398   ,const MatrixPermAggr        &Gc_pa
00399   ,MatrixConvertToSparseEncap  *A_mctse
00400   ) const
00401 {
00402   A_mctse->initialize(
00403     Teuchos::rcp_dynamic_cast<const MatrixExtractSparseElements>(Gc_pa.mat_orig())
00404     ,Teuchos::rcp( init_var_rng_.size()    < n ? &init_var_inv_perm_ : NULL, false )
00405     ,init_var_rng_
00406     ,Teuchos::rcp( init_equ_rng_.size() < m ? &init_equ_inv_perm_ : NULL, false )
00407     ,init_equ_rng_
00408     ,BLAS_Cpp::trans
00409     );
00410 }
00411 
00412 void BasisSystemPermDirectSparse::update_basis_and_auxiliary_matrices(
00413   const MatrixOp& Gc
00414   ,const Teuchos::RCP<DirectSparseSolver::BasisMatrix>& C_bm
00415   ,MatrixOpNonsingAggr *C_aggr
00416   ,MatrixOp* D, MatrixOp* GcUP
00417   ) const
00418 {
00419   using Teuchos::dyn_cast;
00420   // Initialize the aggregate basis matrix object.
00421   C_aggr->initialize(
00422     Gc.sub_view(var_dep_,equ_decomp_)
00423     ,BLAS_Cpp::trans
00424     ,C_bm
00425     ,BLAS_Cpp::no_trans
00426     );
00427   // Compute the auxiliary projected matrices
00428   // Get the concreate type of the direct sensitivity matrix (if one was passed in)
00429   if( D ) {
00430     MultiVectorMutableDense *D_mvd = &dyn_cast<MultiVectorMutableDense>(*D);
00431     TEUCHOS_TEST_FOR_EXCEPT( !(  D  ) ); // ToDo: Throw exception!
00432     // D = -inv(C) * N
00433     D_mvd->initialize(var_dep_.size(),var_indep_.size());
00434     AbstractLinAlgPack::M_StInvMtM(
00435       D_mvd, -1.0, *C_bm, BLAS_Cpp::no_trans
00436       ,*Gc.sub_view(var_indep_,equ_decomp_),BLAS_Cpp::trans // N = Gc(var_indep,equ_decomp)'
00437       );
00438   }
00439   if( GcUP ) {
00440     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00441   }
00442 }
00443 
00444 void BasisSystemPermDirectSparse::do_some_basis_stuff(
00445   const MatrixOp& Gc
00446   ,const Range1D& var_dep, const Range1D& equ_decomp
00447   ,const Teuchos::RCP<DirectSparseSolver::BasisMatrix>& C_bm
00448   ,MatrixOpNonsingAggr *C_aggr
00449   ,MatrixOp* D, MatrixOp* GcUP
00450   )
00451 {
00452   const size_type
00453     n  = Gc.rows(),
00454     m  = Gc.cols();
00455   // Set the ranges
00456   var_dep_      = var_dep;
00457   var_indep_   = ( var_dep.size() == n
00458             ? Range1D::Invalid
00459             : ( var_dep.lbound() == 1
00460               ? Range1D(var_dep.size()+1,n)
00461               : Range1D(1,n-var_dep.size()) ) );
00462   equ_decomp_   = equ_decomp;
00463   equ_undecomp_ = ( equ_decomp.size() == m
00464             ? Range1D::Invalid
00465             : ( equ_decomp.lbound() == 1
00466               ? Range1D(equ_decomp.size()+1,m)
00467               : Range1D(1,m-equ_decomp.size()) ) );
00468   // Set the basis system dimensions
00469   n_     = n;
00470   m_     = m;
00471   r_     = var_dep.size();
00472   Gc_nz_ = Gc.nz();
00473   // Update the aggregate basis matrix and compute the auxiliary projected matrices
00474   update_basis_and_auxiliary_matrices( Gc, C_bm, C_aggr, D, GcUP );
00475 }
00476   
00477 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines