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