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