AbstractLinAlgPack_BasisSystemComposite.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 <assert.h>
00030 
00031 #include "AbstractLinAlgPack_BasisSystemComposite.hpp"
00032 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00033 #include "AbstractLinAlgPack_MatrixComposite.hpp"
00034 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00035 #include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"
00036 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00037 #include "ReleaseResource_ref_count_ptr.hpp"
00038 #include "Teuchos_AbstractFactoryStd.hpp"
00039 #include "Teuchos_dyn_cast.hpp"
00040 #include "Teuchos_TestForException.hpp"
00041 
00042 namespace {
00043 
00044 // Allocates a MultiVectorMutable object given a vector space
00045 // object and the number of columns (num_vecs).
00046 
00047 class AllocatorMultiVectorMutable {
00048 public:
00049   AllocatorMultiVectorMutable(
00050     const AbstractLinAlgPack::VectorSpace::space_ptr_t&  vec_space
00051     ,AbstractLinAlgPack::size_type                       num_vecs
00052     )
00053     :vec_space_(vec_space)
00054     ,num_vecs_(num_vecs)
00055   {}
00056   typedef Teuchos::RCP<
00057     AbstractLinAlgPack::MultiVectorMutable>               ptr_t;
00058   ptr_t allocate() const
00059   {
00060     return vec_space_->create_members(num_vecs_);
00061   }
00062 private:
00063   AbstractLinAlgPack::VectorSpace::space_ptr_t  vec_space_;
00064   AbstractLinAlgPack::size_type                 num_vecs_;
00065 }; // end class AllocatorMultiVectorMutable
00066 
00067 } // end namespace
00068 
00069 namespace AbstractLinAlgPack {
00070 
00071 // Static member functions
00072 
00073 void BasisSystemComposite::initialize_space_x(
00074   const VectorSpace::space_ptr_t    &space_xD
00075   ,const VectorSpace::space_ptr_t   &space_xI
00076   ,Range1D                          *var_dep
00077   ,Range1D                          *var_indep
00078   ,VectorSpace::space_ptr_t         *space_x
00079   )
00080 {
00081   namespace mmp = MemMngPack;
00082 #ifdef TEUCHOS_DEBUG
00083   TEST_FOR_EXCEPTION(
00084     space_xD.get() == NULL, std::invalid_argument
00085     ,"BasisSystemComposite::initialize_space_x(...): Error!" );
00086   TEST_FOR_EXCEPTION(
00087       var_dep == NULL, std::invalid_argument
00088     ,"BasisSystemComposite::initialize_space_x(...): Error!" );
00089   TEST_FOR_EXCEPTION(
00090     space_xI.get() != NULL && var_indep == NULL, std::invalid_argument
00091     ,"BasisSystemComposite::initialize_space_x(...): Error!" );
00092   TEST_FOR_EXCEPTION(
00093     space_x == NULL, std::invalid_argument
00094     ,"BasisSystemComposite::initialize_space_x(...): Error!" );
00095 #endif
00096   *var_dep   = Range1D(1,space_xD->dim());
00097   if(space_xI.get()) *var_indep = Range1D(var_dep->ubound()+1,var_dep->ubound()+space_xI->dim());
00098   else               *var_indep = Range1D::Invalid;
00099   if(space_xI.get()) {
00100     const VectorSpace::space_ptr_t
00101       vec_spaces[2] = { space_xD, space_xI };
00102     *space_x   = Teuchos::rcp(new VectorSpaceBlocked(vec_spaces,2));
00103   }
00104   else {
00105     *space_x = space_xD;
00106   }
00107 }
00108 
00109 const BasisSystemComposite::fcty_Gc_ptr_t
00110 BasisSystemComposite::factory_Gc()
00111 {
00112   namespace mmp = MemMngPack;
00113   return Teuchos::rcp( new Teuchos::AbstractFactoryStd<MatrixOp,MatrixComposite>() );
00114 }
00115 
00116 void BasisSystemComposite::initialize_Gc(
00117   const VectorSpace::space_ptr_t    &space_x
00118   ,const Range1D                    &var_dep
00119   ,const Range1D                    &var_indep
00120   ,const VectorSpace::space_ptr_t   &space_c
00121   ,const C_ptr_t                    &C
00122   ,const N_ptr_t                    &N
00123   ,MatrixOp                         *Gc
00124   )
00125 {
00126   namespace mmp = MemMngPack;
00127   using Teuchos::dyn_cast;
00128 #ifdef TEUCHOS_DEBUG
00129   TEST_FOR_EXCEPTION(
00130     space_x.get() == NULL, std::invalid_argument
00131     ,"BasisSystemComposite::initialize_Gc(...): Error!" );
00132   TEST_FOR_EXCEPTION(
00133     space_c.get() == NULL, std::invalid_argument
00134     ,"BasisSystemComposite::initialize_Gc(...): Error!" );
00135 #endif
00136   const size_type
00137     n            = space_x->dim(),
00138     m            = space_c->dim(),
00139     var_dep_size = var_dep.size();
00140 #ifdef TEUCHOS_DEBUG
00141   TEST_FOR_EXCEPTION(
00142     C.get() == NULL, std::invalid_argument
00143     ,"BasisSystemComposite::initialize_Gc(...): Error!" );
00144   TEST_FOR_EXCEPTION(
00145     var_dep_size < n && N.get() == NULL, std::invalid_argument
00146     ,"BasisSystemComposite::initialize_Gc(...): Error!" );
00147   TEST_FOR_EXCEPTION(
00148     Gc == NULL, std::invalid_argument
00149     ,"BasisSystemComposite::initialize_Gc(...): Error!" );
00150 #endif
00151 
00152   MatrixComposite
00153     &Gc_comp = dyn_cast<MatrixComposite>(*Gc);
00154     
00155   //
00156   // Gc = [ C'; N' ]
00157   //
00158   
00159   Gc_comp.reinitialize(n,m);
00160   // Add the C matrix object
00161   typedef Teuchos::RCP<mmp::ReleaseResource_ref_count_ptr<MatrixOpNonsing> > C_rr_ptr_ptr_t;
00162   C_rr_ptr_ptr_t
00163     C_rr_ptr_ptr = Teuchos::rcp(new mmp::ReleaseResource_ref_count_ptr<MatrixOpNonsing>(C));
00164   Gc_comp.add_matrix(
00165     var_dep.lbound()-1, 0    // row_offset, col_offset
00166     ,1.0                     // alpha
00167     ,C_rr_ptr_ptr->ptr.get() // A
00168     ,C_rr_ptr_ptr            // A_release
00169     ,BLAS_Cpp::trans         // A_trans
00170     );
00171   if( n > m ) {
00172     // Add the N matrix object
00173     typedef Teuchos::RCP<mmp::ReleaseResource_ref_count_ptr<MatrixOp> > N_rr_ptr_ptr_t;
00174     N_rr_ptr_ptr_t
00175       N_rr_ptr_ptr = Teuchos::rcp(new mmp::ReleaseResource_ref_count_ptr<MatrixOp>(N));
00176     Gc_comp.add_matrix(
00177       var_indep.lbound()-1, 0  // row_offset, col_offset
00178       ,1.0                     // alpha
00179       ,N_rr_ptr_ptr->ptr.get() // A
00180       ,N_rr_ptr_ptr            // A_release
00181       ,BLAS_Cpp::trans         // A_trans
00182       );
00183   }
00184   // Finish construction
00185   Gc_comp.finish_construction( space_x, space_c );
00186 }
00187 
00188 void BasisSystemComposite::get_C_N(
00189   MatrixOp               *Gc
00190   ,MatrixOpNonsing       **C
00191   ,MatrixOp              **N
00192   )
00193 {
00194   using Teuchos::dyn_cast;
00195 #ifdef TEUCHOS_DEBUG
00196   TEST_FOR_EXCEPTION(
00197     Gc == NULL, std::invalid_argument
00198     ,"BasisSystemComposite::get_C_N(...): Error!" );
00199 #endif
00200   const size_type
00201     n = Gc->rows(),
00202     m = Gc->cols();
00203 #ifdef TEUCHOS_DEBUG
00204   TEST_FOR_EXCEPTION(
00205     C == NULL, std::invalid_argument
00206     ,"BasisSystemComposite::get_C_N(...): Error!" );
00207   TEST_FOR_EXCEPTION(
00208     n > m && N == NULL, std::invalid_argument
00209     ,"BasisSystemComposite::get_C_N(...): Error!" );
00210 #endif
00211   // Get reference to concrete Gc matrix subclass
00212   MatrixComposite
00213     &Gc_comp = dyn_cast<MatrixComposite>(*Gc);
00214   // Get referencs to the aggregate C and N matrtices
00215   MatrixComposite::matrix_list_t::const_iterator
00216     mat_itr = Gc_comp.matrices_begin(),
00217     mat_end = Gc_comp.matrices_end();
00218   if( mat_itr != mat_end ) {
00219     TEST_FOR_EXCEPT( !( mat_itr != mat_end ) );
00220     *C = &dyn_cast<MatrixOpNonsing>(
00221       const_cast<MatrixOp&>(*(mat_itr++)->A_) );
00222     if( n > m ) {
00223       TEST_FOR_EXCEPT( !( mat_itr != mat_end ) );
00224       *N = &const_cast<MatrixOp&>(*(mat_itr++)->A_);
00225     }
00226     TEST_FOR_EXCEPT( !( mat_itr == mat_end ) );
00227   }
00228   else {
00229     *C = NULL;
00230     *N = NULL;
00231   }
00232 }
00233 
00234 void BasisSystemComposite::get_C_N(
00235   const MatrixOp               &Gc
00236   ,const MatrixOpNonsing       **C
00237   ,const MatrixOp              **N
00238   )
00239 {
00240   using Teuchos::dyn_cast;
00241   const size_type
00242     n = Gc.rows(),
00243     m = Gc.cols();
00244 #ifdef TEUCHOS_DEBUG
00245   TEST_FOR_EXCEPTION(
00246     C == NULL, std::invalid_argument
00247     ,"BasisSystemComposite::get_C_N(...): Error!" );
00248   TEST_FOR_EXCEPTION(
00249     n > m && N == NULL, std::invalid_argument
00250     ,"BasisSystemComposite::get_C_N(...): Error!" );
00251 #endif
00252   // Get reference to concrete Gc matrix subclass
00253   const AbstractLinAlgPack::MatrixComposite
00254     &Gc_comp = dyn_cast<const AbstractLinAlgPack::MatrixComposite>(Gc);
00255   // Get referencs to the aggregate C and N matrtices
00256   MatrixComposite::matrix_list_t::const_iterator
00257     mat_itr = Gc_comp.matrices_begin(),
00258     mat_end = Gc_comp.matrices_end();
00259   if( mat_itr != mat_end ) {
00260     TEST_FOR_EXCEPT( !( mat_itr != mat_end ) );
00261     *C = &dyn_cast<const MatrixOpNonsing>(*(mat_itr++)->A_);
00262     if( n > m ) {
00263       TEST_FOR_EXCEPT( !( mat_itr != mat_end ) );
00264       *N = &dyn_cast<const MatrixOp>(*(mat_itr++)->A_);
00265     }
00266     TEST_FOR_EXCEPT( !( mat_itr == mat_end ) );
00267   }
00268   else {
00269     TEST_FOR_EXCEPTION(
00270       true, std::invalid_argument
00271       ,"BasisSystemComposite::get_C_N(...): Error, "
00272       "The Gc matrix object has not been initialized with C and N!" );
00273   }
00274 }
00275 
00276 // Constructors / initializers
00277 
00278 BasisSystemComposite::BasisSystemComposite()
00279   :BasisSystem(Teuchos::null,Teuchos::null)
00280 {}
00281 
00282 BasisSystemComposite::BasisSystemComposite(
00283   const VectorSpace::space_ptr_t       &space_x
00284   ,const VectorSpace::space_ptr_t      &space_c
00285   ,const mat_nonsing_fcty_ptr_t        &factory_C
00286   ,const mat_sym_fcty_ptr_t            &factory_transDtD
00287   ,const mat_sym_nonsing_fcty_ptr_t    &factory_S
00288   )
00289   :BasisSystem(Teuchos::null,Teuchos::null)
00290 {
00291   namespace mmp = MemMngPack;
00292   const size_type n = space_x->dim(), m = space_c->dim();
00293   this->initialize(
00294     space_x,Range1D(1,space_c->dim()),(n>m ? Range1D(space_c->dim()+1,space_x->dim()) : Range1D::Invalid)
00295     ,space_c,factory_C,factory_transDtD,factory_S
00296     );
00297 }
00298 
00299 BasisSystemComposite::BasisSystemComposite(
00300   const VectorSpace::space_ptr_t       &space_x
00301   ,const Range1D                       &var_dep
00302   ,const Range1D                       &var_indep
00303   ,const VectorSpace::space_ptr_t      &space_c
00304   ,const mat_nonsing_fcty_ptr_t        &factory_C
00305   ,const mat_sym_fcty_ptr_t            &factory_transDtD
00306   ,const mat_sym_nonsing_fcty_ptr_t    &factory_S
00307   ,const mat_fcty_ptr_t                &factory_D
00308   )
00309   :BasisSystem(Teuchos::null,Teuchos::null)
00310 {
00311   this->initialize(
00312     space_x,var_dep,var_indep,space_c,factory_C,factory_transDtD,factory_S
00313     ,factory_D
00314     );
00315 }
00316   
00317 void BasisSystemComposite::initialize(
00318   const VectorSpace::space_ptr_t       &space_x
00319   ,const Range1D                       &var_dep
00320   ,const Range1D                       &var_indep
00321   ,const VectorSpace::space_ptr_t      &space_c
00322   ,const mat_nonsing_fcty_ptr_t        &factory_C
00323   ,const mat_sym_fcty_ptr_t            &factory_transDtD
00324   ,const mat_sym_nonsing_fcty_ptr_t    &factory_S
00325   ,const mat_fcty_ptr_t                &factory_D
00326   )
00327 {
00328   namespace mmp = MemMngPack;
00329 #ifdef TEUCHOS_DEBUG
00330   TEST_FOR_EXCEPTION(
00331     space_x.get() == NULL, std::invalid_argument
00332     ,"BasisSystemComposite::initialize(...): Error!" );
00333   TEST_FOR_EXCEPTION(
00334     space_c.get() == NULL, std::invalid_argument
00335     ,"BasisSystemComposite::initialize(...): Error!" );
00336 #endif
00337   const size_type n = space_x->dim(), m = space_c->dim();
00338 #ifdef TEUCHOS_DEBUG
00339   TEST_FOR_EXCEPTION(
00340     var_dep.size() + var_indep.size() != space_x->dim(), std::invalid_argument
00341     ,"BasisSystemComposite::initialize(...): Error!" );
00342   TEST_FOR_EXCEPTION(
00343     n > m && var_dep.lbound() < var_indep.lbound() && (var_dep.lbound() != 1 || var_dep.ubound()+1 != var_indep.lbound())
00344     , std::invalid_argument
00345     ,"BasisSystemComposite::initialize(...): Error!" );
00346   TEST_FOR_EXCEPTION(
00347     n > m && var_dep.lbound() >= var_indep.lbound() && (var_indep.lbound() != 1 || var_indep.ubound()+1 != var_dep.lbound())
00348     , std::invalid_argument
00349     ,"BasisSystemComposite::initialize(...): Error!" );
00350   TEST_FOR_EXCEPTION(
00351     factory_C.get() == NULL, std::invalid_argument
00352     ,"BasisSystemComposite::initialize(...): Error!" );
00353 #endif
00354 
00355   space_x_         = space_x;
00356   var_dep_         = var_dep;
00357   var_indep_       = var_indep;
00358   space_c_         = space_c;
00359   factory_C_       = factory_C;
00360   factory_D_       = factory_D;
00361   if( factory_D_.get() == NULL ) {
00362     factory_D_ = Teuchos::abstractFactoryStd<MatrixOp,MultiVectorMutable>(
00363       AllocatorMultiVectorMutable(space_x_->sub_space(var_dep),var_indep.size() ) );
00364   }
00365   BasisSystem::initialize(factory_transDtD,factory_S);
00366 }
00367 
00368 void BasisSystemComposite::set_uninitialized()
00369 {
00370   namespace mmp = MemMngPack;
00371 
00372   space_x_         = Teuchos::null;
00373   var_dep_         = Range1D::Invalid;
00374   var_indep_       = Range1D::Invalid;
00375   factory_C_       = Teuchos::null;
00376   factory_D_       = Teuchos::null;
00377 }
00378 
00379 const VectorSpace::space_ptr_t&
00380 BasisSystemComposite::space_x() const
00381 {
00382   return space_x_;
00383 }
00384 
00385 const VectorSpace::space_ptr_t&
00386 BasisSystemComposite::space_c() const
00387 {
00388   return space_c_;
00389 }
00390 
00391 // To be overridden by subclasses
00392 
00393 void BasisSystemComposite::update_D(
00394   const MatrixOpNonsing       &C
00395   ,const MatrixOp             &N
00396   ,MatrixOp                   *D
00397   ,EMatRelations              mat_rel
00398   ) const
00399 {
00400   using LinAlgOpPack::M_StInvMtM;
00401   M_StInvMtM( D, -1.0, C, BLAS_Cpp::no_trans, N, BLAS_Cpp::no_trans ); // D = -inv(C)*N
00402 }
00403 
00404 // Overridden from BasisSytem
00405 
00406 const BasisSystem::mat_nonsing_fcty_ptr_t
00407 BasisSystemComposite::factory_C() const
00408 {
00409   return factory_C_;
00410 }
00411 
00412 const BasisSystem::mat_fcty_ptr_t
00413 BasisSystemComposite::factory_D() const
00414 {
00415   return factory_D_;
00416 }
00417 
00418 Range1D BasisSystemComposite::var_dep() const
00419 {
00420   return var_dep_;
00421 }
00422 
00423 Range1D BasisSystemComposite::var_indep() const
00424 {
00425   return var_indep_;
00426 }
00427 
00428 void BasisSystemComposite::update_basis(
00429   const MatrixOp          &Gc
00430   ,MatrixOpNonsing        *C
00431   ,MatrixOp               *D
00432   ,MatrixOp               *GcUP
00433   ,EMatRelations          mat_rel
00434   ,std::ostream           *out
00435   ) const
00436 {
00437   using Teuchos::dyn_cast;
00438   const index_type
00439     n  = var_dep_.size() + var_indep_.size(),
00440     m  = var_dep_.size();
00441   TEST_FOR_EXCEPTION(
00442     n == 0, std::logic_error
00443     ,"BasisSystemComposite::update_basis(...): Error, this must be initialized first!" );
00444   TEST_FOR_EXCEPTION(
00445     C == NULL && ( n > m ? D == NULL : false ), std::logic_error
00446     ,"BasisSystemComposite::update_basis(...): Error, C or D must be non-NULL!" );
00447   // Get references to the aggregate C and N matrices
00448   const MatrixOpNonsing
00449     *C_aggr = NULL;
00450   const MatrixOp
00451     *N_aggr = NULL;
00452   get_C_N( Gc, &C_aggr, &N_aggr );
00453   // Setup C
00454   if( C ) {
00455     *C = *C_aggr;  // Assignment had better work!
00456   }
00457   // Compute D
00458   if( D ) {
00459     update_D(*C_aggr,*N_aggr,D,mat_rel);
00460   }
00461 }
00462 
00463 } // end namespace AbstractLinAlgPack

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