AbstractLinAlgPack_MatrixComposite.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 "AbstractLinAlgPack_MatrixComposite.hpp"
00030 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00031 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00032 #include "AbstractLinAlgPack_VectorMutableBlocked.hpp"
00033 #include "AbstractLinAlgPack_AssertOp.hpp"
00034 //#include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00035 #include "Teuchos_Workspace.hpp"
00036 #include "Teuchos_TestForException.hpp"
00037 #include "ProfileHackPack_profile_hack.hpp"
00038 
00039 namespace {
00040 
00041 // Get an element from a vector (two versions)
00042 
00043 inline
00044 AbstractLinAlgPack::value_type
00045 get_element( const AbstractLinAlgPack::Vector& v, AbstractLinAlgPack::index_type i )
00046 {
00047    return v.get_ele(i);
00048 }
00049 
00050 inline
00051 AbstractLinAlgPack::value_type
00052 get_element( const AbstractLinAlgPack::SpVectorSlice& v, AbstractLinAlgPack::index_type i )
00053 {
00054   const AbstractLinAlgPack::SpVectorSlice::element_type
00055     *ele = v.lookup_element(i);
00056   return ele != NULL ? ele->value() : 0.0;
00057 }
00058 
00059 // Get a view of a vector (two versions)
00060 
00061 inline
00062 Teuchos::RCP<const AbstractLinAlgPack::Vector>
00063 get_view(
00064   const AbstractLinAlgPack::Vector& v
00065   ,AbstractLinAlgPack::index_type l
00066   ,AbstractLinAlgPack::index_type u
00067   )
00068 {
00069    return v.sub_view(l,u);
00070 }
00071 
00072 inline
00073 Teuchos::RCP<const AbstractLinAlgPack::SpVectorSlice>
00074 get_view(
00075   const AbstractLinAlgPack::SpVectorSlice& v
00076   ,AbstractLinAlgPack::index_type l
00077   ,AbstractLinAlgPack::index_type u
00078   )
00079 {
00080   return Teuchos::rcp( new AbstractLinAlgPack::SpVectorSlice( v(l,u) ) );
00081 }
00082 
00083 // Perform a matrix vector multiplication
00084 
00085 template<class V>
00086 void Vp_StMtV_imp(
00087   AbstractLinAlgPack::VectorMutable* y, AbstractLinAlgPack::value_type a
00088   ,AbstractLinAlgPack::size_type M_rows, AbstractLinAlgPack::size_type M_cols
00089   ,const AbstractLinAlgPack::MatrixComposite::matrix_list_t& mat_list
00090   ,const AbstractLinAlgPack::MatrixComposite::vector_list_t& vec_list
00091   ,BLAS_Cpp::Transp M_trans
00092   ,const V& x, AbstractLinAlgPack::value_type b
00093   )
00094 {
00095   using BLAS_Cpp::no_trans;
00096   using BLAS_Cpp::trans;
00097   using BLAS_Cpp::rows;
00098   using BLAS_Cpp::cols;
00099   using AbstractLinAlgPack::dot;  // We should not have to do this but some compiles &%!#$
00100   typedef AbstractLinAlgPack::MatrixComposite::matrix_list_t  mat_list_t;
00101   typedef AbstractLinAlgPack::MatrixComposite::vector_list_t  vec_list_t;
00102 
00103     AbstractLinAlgPack::Vt_S( y, b );  // Will take care of b == 0.0
00104 
00105   if( vec_list.size() ) {
00106     for( vec_list_t::const_iterator itr = vec_list.begin(); itr != vec_list.end(); ++itr ) {
00107       const BLAS_Cpp::Transp
00108         op_v_trans = ( M_trans == itr->v_trans_ ? no_trans : trans );
00109       const AbstractLinAlgPack::index_type
00110         r = ( M_trans == no_trans ? itr->r_l_ : itr->c_l_ ),
00111         c = ( M_trans == no_trans ? itr->c_l_ : itr->r_l_ );
00112       if( itr->rng_G_.full_range() ) { // op(v)
00113         if( op_v_trans == no_trans ) {
00114           //
00115           //         [ y1 ]         [              ]   [ x1 ]
00116           // r:r+n-1 [ y2 ] +=  a * [   beta * v   ] * [ x2 ] c:c 
00117           //         [ y3 ]         [              ]   [ x3 ]
00118           //                            \______/        
00119           //                              c:c
00120           // =>
00121           //
00122           // y(r:r+n-1) += (a * beta * x(c)) * v
00123           // 
00124           AbstractLinAlgPack::Vp_StV( y->sub_view(r,r+itr->v_->dim()-1).get(), a * itr->beta_ * get_element(x,c), *itr->v_ );
00125         }
00126         else {
00127           //
00128           //     [ y1 ]        [                ]   [ x1 ]
00129           // r:r [ y2 ] += a * [     beta*v'    ] * [ x2 ] c:c+n-1
00130           //     [ y3 ]        [                ]   [ x3 ]
00131           //                         \_____/  
00132           //                         c:c+n-1
00133           // =>
00134           //
00135           // y(r) += a * beta * v'*x(c,c+n-1)
00136           //
00137 //          y->set_ele( r, y->get_ele(r) + a * itr->beta_ * dot( *itr->v_, *get_view(x,c,c+itr->v_->dim()-1) ) );
00138           TEST_FOR_EXCEPT(true); // ToDo: Implement the above method in VectorStdOps for Vector,SpVectorSlice!
00139         }
00140       }
00141       else { // op(op(G)*v) or op(v(rng_G))
00142         TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
00143       }
00144     }
00145   }
00146   if( mat_list.size() ) {
00147     for( mat_list_t::const_iterator itr = mat_list.begin(); itr != mat_list.end(); ++itr ) {
00148       const AbstractLinAlgPack::index_type
00149         rl = rows(itr->r_l_,itr->c_l_,M_trans),
00150         ru = rows(itr->r_u_,itr->c_u_,M_trans),
00151         cl = cols(itr->r_l_,itr->c_l_,M_trans),
00152         cu = cols(itr->r_u_,itr->c_u_,M_trans);
00153       const BLAS_Cpp::Transp
00154         op_P_trans = ( M_trans == itr->P_trans_ ? no_trans : trans ),
00155         op_A_trans = ( M_trans == itr->A_trans_ ? no_trans : trans ),
00156         op_Q_trans = ( M_trans == itr->Q_trans_ ? no_trans : trans );
00157       if( itr->rng_P_.full_range() && itr->rng_Q_.full_range() ) { // op(A)
00158         //
00159         //       [ y1 ]        [                        ]   [ x1 ]
00160         // rl:ru [ y2 ] += a * [    alpha * op(op(A))   ] * [ x2 ] cl:cu
00161         //       [ y3 ]        [                        ]   [ x3 ]
00162         //                          \_______________/
00163         //                               cl:cu
00164         // =>
00165         //
00166         // y(rl:ru) += a * alpha * op(op(A)) * x(cl:cu)
00167         //
00168         AbstractLinAlgPack::Vp_StMtV( y->sub_view(rl,ru).get(), a * itr->alpha_, *itr->A_, op_A_trans, *get_view(x,cl,cu) );
00169       }
00170       else {
00171         if( itr->A_ == NULL ) { // op(P)
00172           TEST_FOR_EXCEPT( !(  itr->P_.get() && !itr->P_->is_identity()  ) );
00173           //
00174           //       [ y1 ]        [                        ]   [ x1 ]
00175           // rl:ru [ y2 ] += a * [    alpha * op(op(P))   ] * [ x2 ] cl:cu
00176           //       [ y3 ]        [                        ]   [ x3 ]
00177           //                          \_______________/
00178           //                               cl:cu
00179           // =>
00180           //
00181           // y(rl:ru) += a * alpha * op(op(P)) * x(cl:cu)
00182           //
00183 //          AbstractLinAlgPack::Vp_StMtV( y->sub_view(rl,ru).get(), a * itr->alpha_, itr->P_, op_P_trans, *get_view(x,cl,cu) );
00184           TEST_FOR_EXCEPT(true); // ToDo: Implement the above method properly!
00185         }
00186         else { // op(P)*op(A)*op(Q)  [or some simplification]
00187           //
00188           //       [ y1 ]        [                                ]   [ x1 ]
00189           // rl:ru [ y2 ] += a * [    alpha * op(P)*op(A)*op(Q)   ] * [ x2 ] cl:cu
00190           //       [ y3 ]        [                                ]   [ x3 ]
00191           //                          \______________________/
00192           //                                    cl:cu
00193           // =>
00194           //
00195           // y(rl:ru) += a * alpha * op(P)*op(A)*op(Q) * x(cl:cu)
00196           //
00197           if( !itr->rng_P_.full_range() && !itr->rng_Q_.full_range() ) { // op(A)(rng_Q,rng_Q)
00198             AbstractLinAlgPack::Vp_StMtV(
00199               y->sub_view(rl,ru).get(), a * itr->alpha_
00200               ,*itr->A_->sub_view(
00201                 itr->A_trans_==no_trans  ? itr->rng_P_ : itr->rng_Q_
00202                 ,itr->A_trans_==no_trans ? itr->rng_Q_ : itr->rng_P_ )
00203               ,op_A_trans
00204               ,*get_view(x,cl,cu)
00205               );
00206           }
00207           else {
00208             TEST_FOR_EXCEPT(true); // ToDo: Implement when needed!
00209           }
00210         }
00211       }
00212     }
00213   }
00214 }
00215 
00216 // Define a comparison object for comparing SubVectorEntry or SubMatrixEntry
00217 // objects for storting by row or by column
00218 template<class T>
00219 class CompSubEntry {
00220 public:
00221   enum EByRowCol { BY_ROW, BY_COL };
00222   CompSubEntry(enum EByRowCol by_row_col)
00223     : by_row_col_(by_row_col)
00224     {}
00225   bool operator()( const T& e1, const T& e2 )
00226     {
00227       return
00228         ( by_row_col_ == BY_ROW
00229           ? e1.r_l_ < e2.r_l_
00230           : e1.c_l_ < e2.c_l_
00231           );
00232     }
00233 private:
00234   EByRowCol  by_row_col_;
00235   CompSubEntry(); // not defined and not to be called
00236 }; // end class CompSubVectorEntry
00237 
00238 
00239 } // end namespace
00240 
00241 namespace AbstractLinAlgPack {
00242 
00243 MatrixComposite::MatrixComposite( size_type rows, size_type cols )
00244 {
00245   reinitialize(rows,cols);
00246 }
00247 
00248 void MatrixComposite::reinitialize( size_type rows, size_type cols )
00249 {
00250   namespace rcp = MemMngPack;
00251   
00252   fully_constructed_ = true;
00253   rows_ = rows;
00254   cols_ = cols;
00255   if(matrix_list_.size()) {
00256     matrix_list_.erase(matrix_list_.begin(),matrix_list_.end());
00257   }
00258   if(vector_list_.size()) {
00259     vector_list_.erase(vector_list_.begin(),vector_list_.end());
00260   }
00261   space_rows_  = Teuchos::null;
00262   space_cols_  = Teuchos::null;
00263 }
00264 
00265 void MatrixComposite::add_vector(
00266   size_type                      row_offset
00267   ,size_type                     col_offset
00268   ,value_type                    beta
00269   ,const GenPermMatrixSlice      *G
00270   ,const release_resource_ptr_t  &G_release
00271   ,BLAS_Cpp::Transp              G_trans
00272   ,const Vector            *v
00273   ,const release_resource_ptr_t  &v_release
00274   ,BLAS_Cpp::Transp              v_trans
00275   )
00276 {
00277   fully_constructed_ = false;
00278   TEST_FOR_EXCEPT(true); // ToDo: Finish!
00279 }
00280 
00281 void MatrixComposite::add_vector(
00282   size_type                      row_offset
00283   ,size_type                     col_offset
00284   ,value_type                    beta
00285   ,const Range1D                 &rng_G
00286   ,const Vector            *v
00287   ,const release_resource_ptr_t  &v_release
00288   ,BLAS_Cpp::Transp              v_trans
00289   )
00290 {
00291   fully_constructed_ = false;
00292   TEST_FOR_EXCEPT(true); // ToDo: Finish!
00293 }
00294 
00295 void MatrixComposite::add_vector(
00296   size_type                      row_offset
00297   ,size_type                     col_offset
00298   ,value_type                    beta
00299   ,const Vector            *v
00300   ,const release_resource_ptr_t  &v_release
00301   ,BLAS_Cpp::Transp              v_trans
00302   )
00303 {
00304   namespace rcp = MemMngPack;
00305 
00306   TEST_FOR_EXCEPT( !(  beta != 0.0  ) );
00307   TEST_FOR_EXCEPT( !(  v != NULL  ) );
00308   fully_constructed_ = false;
00309   if( v_trans == BLAS_Cpp::no_trans ) {
00310     TEST_FOR_EXCEPT( !(  row_offset + v->dim() <= rows_  ) );
00311     TEST_FOR_EXCEPT( !(  col_offset + 1 <= cols_  ) );
00312   }
00313   else {
00314     TEST_FOR_EXCEPT( !(  row_offset + 1 <= rows_  ) );
00315     TEST_FOR_EXCEPT( !(  col_offset + v->dim() <= cols_  ) );
00316   }
00317 
00318   vector_list_.push_back(
00319     SubVectorEntry(
00320       row_offset+1,col_offset+1,beta,Range1D()
00321       ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
00322       ,v,v_release,v_trans ) );
00323 }
00324 
00325 void MatrixComposite::remove_vector( vector_list_t::iterator itr )
00326 {
00327   fully_constructed_ = false;
00328   vector_list_.erase(itr);
00329 }
00330 
00331 void MatrixComposite::add_matrix(
00332   size_type                      row_offset
00333   ,size_type                     col_offset
00334   ,value_type                    alpha
00335   ,const GenPermMatrixSlice      *P
00336   ,const release_resource_ptr_t  &P_release
00337   ,BLAS_Cpp::Transp              P_trans
00338   ,const MatrixOp            *A
00339   ,const release_resource_ptr_t  &A_release
00340   ,BLAS_Cpp::Transp              A_trans
00341   ,const GenPermMatrixSlice      *Q
00342   ,const release_resource_ptr_t  &Q_release
00343   ,BLAS_Cpp::Transp              Q_trans
00344   )
00345 {
00346   fully_constructed_ = false;
00347   TEST_FOR_EXCEPT(true); // ToDo: Finish!
00348 }
00349 
00350 void MatrixComposite::add_matrix(
00351   size_type                      row_offset
00352   ,size_type                     col_offset
00353   ,value_type                    alpha
00354   ,const Range1D                 &rng_P
00355   ,const MatrixOp            *A
00356   ,const release_resource_ptr_t  &A_release
00357   ,BLAS_Cpp::Transp              A_trans
00358   ,const GenPermMatrixSlice      *Q
00359   ,const release_resource_ptr_t  &Q_release
00360   ,BLAS_Cpp::Transp              Q_trans
00361   )
00362 {
00363   fully_constructed_ = false;
00364   TEST_FOR_EXCEPT(true); // ToDo: Finish!
00365 }
00366 
00367 void MatrixComposite::add_matrix(
00368   size_type                      row_offset
00369   ,size_type                     col_offset
00370   ,value_type                    alpha
00371   ,const GenPermMatrixSlice      *P
00372   ,const release_resource_ptr_t  &P_release
00373   ,BLAS_Cpp::Transp              P_trans
00374   ,const MatrixOp            *A
00375   ,const release_resource_ptr_t  &A_release
00376   ,BLAS_Cpp::Transp              A_trans
00377   ,const Range1D                 &rng_Q
00378   )
00379 {
00380   fully_constructed_ = false;
00381   TEST_FOR_EXCEPT(true); // ToDo: Finish!
00382 }
00383 
00384 void MatrixComposite::add_matrix(
00385   size_type                      row_offset
00386   ,size_type                     col_offset
00387   ,value_type                    alpha
00388   ,const Range1D                 &rng_P_in
00389   ,const MatrixOp            *A
00390   ,const release_resource_ptr_t  &A_release
00391   ,BLAS_Cpp::Transp              A_trans
00392   ,const Range1D                 &rng_Q_in
00393   )
00394 {
00395   namespace rcp = MemMngPack;
00396   using BLAS_Cpp::rows;
00397   using BLAS_Cpp::cols;
00398   using RangePack::full_range;
00399 
00400   TEST_FOR_EXCEPTION(
00401     alpha == 0.0, std::invalid_argument
00402     ,"MatrixComposite::add_matrix(...) : Error!" );
00403   TEST_FOR_EXCEPTION(
00404     A == NULL, std::invalid_argument
00405     ,"MatrixComposite::add_matrix(...) : Error!" );
00406 
00407   const size_type
00408     A_rows   = A->rows(),
00409     A_cols   = A->cols(),
00410     opA_rows = rows(A_rows,A_cols,A_trans),
00411     opA_cols = cols(A_rows,A_cols,A_trans);
00412   const Range1D
00413     rng_P = full_range(rng_P_in,1,opA_rows),
00414     rng_Q = full_range(rng_Q_in,1,opA_cols);
00415   const size_type
00416     opPopAopQ_rows = rng_P.size(),
00417     opPopAopQ_cols = rng_Q.size();
00418 
00419   TEST_FOR_EXCEPTION(
00420     row_offset + opPopAopQ_rows > rows_, std::invalid_argument
00421     ,"MatrixComposite::add_matrix(...) : Error!" );
00422   TEST_FOR_EXCEPTION(
00423     col_offset + opPopAopQ_cols > cols_, std::invalid_argument
00424     ,"MatrixComposite::add_matrix(...) : Error!" );
00425 
00426   fully_constructed_ = false;
00427   
00428   matrix_list_.push_back(
00429     SubMatrixEntry(
00430       row_offset+1,row_offset+opPopAopQ_rows
00431       ,col_offset+1,col_offset+opPopAopQ_cols
00432       ,alpha
00433       ,rng_P
00434       ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
00435       ,A,A_release,A_trans
00436       ,rng_Q
00437       ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
00438       )
00439     );
00440 
00441   // ToDo: We could create a sub-view of the matrix using rng_P and rng_Q
00442   // and then store the sub-view in the data structure?  However, this
00443   // would confuse an external client who wanted to iterate through
00444   // the matrix list using the public iterators.
00445 }
00446 
00447 void MatrixComposite::add_matrix(
00448   size_type                      row_offset
00449   ,size_type                     col_offset
00450   ,value_type                    alpha
00451   ,const MatrixOp            *A
00452   ,const release_resource_ptr_t  &A_release
00453   ,BLAS_Cpp::Transp              A_trans
00454   )
00455 {
00456   namespace rcp = MemMngPack;
00457   using BLAS_Cpp::rows;
00458   using BLAS_Cpp::cols;
00459 
00460   TEST_FOR_EXCEPTION(
00461     alpha == 0.0, std::invalid_argument
00462     ,"MatrixComposite::add_matrix(...) : Error!" );
00463   TEST_FOR_EXCEPTION(
00464     A == NULL, std::invalid_argument
00465     ,"MatrixComposite::add_matrix(...) : Error!" );
00466 
00467   const size_type
00468     A_rows   = A->rows(),
00469     A_cols   = A->cols(),
00470     opA_rows = rows(A_rows,A_cols,A_trans),
00471     opA_cols = cols(A_rows,A_cols,A_trans);
00472 
00473   TEST_FOR_EXCEPTION(
00474     row_offset + opA_rows > rows_, std::invalid_argument
00475     ,"MatrixComposite::add_matrix(...) : Error!" );
00476   TEST_FOR_EXCEPTION(
00477     col_offset + opA_cols > cols_, std::invalid_argument
00478     ,"MatrixComposite::add_matrix(...) : Error!" );
00479 
00480   fully_constructed_ = false;
00481 
00482   matrix_list_.push_back(
00483     SubMatrixEntry(
00484       row_offset+1,row_offset+opA_rows
00485       ,col_offset+1,col_offset+opA_cols
00486       ,alpha
00487       ,Range1D()
00488       ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
00489       ,A,A_release,A_trans
00490       ,Range1D()
00491       ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
00492       )
00493     );
00494 }
00495 
00496 void MatrixComposite::add_matrix(
00497   size_type                      row_offset
00498   ,size_type                     col_offset
00499   ,value_type                    alpha
00500   ,const GenPermMatrixSlice      *P
00501   ,const release_resource_ptr_t  &P_release
00502   ,BLAS_Cpp::Transp              P_trans
00503   )
00504 {
00505   namespace rcp = MemMngPack;
00506   using BLAS_Cpp::rows;
00507   using BLAS_Cpp::cols;
00508 
00509   TEST_FOR_EXCEPT( !(  alpha != 0.0  ) );
00510   TEST_FOR_EXCEPT( !(  P != NULL  ) );
00511 
00512   fully_constructed_ = false;
00513 
00514   const size_type
00515     P_rows   = P->rows(),
00516     P_cols   = P->cols(),
00517     opP_rows = rows(P_rows,P_cols,P_trans),
00518     opP_cols = cols(P_rows,P_cols,P_trans);
00519 
00520   TEST_FOR_EXCEPT( !(  row_offset + opP_rows <= rows_  ) );
00521   TEST_FOR_EXCEPT( !(  col_offset + opP_cols <= cols_  ) );
00522 
00523   matrix_list_.push_back(
00524     SubMatrixEntry(
00525       row_offset+1,row_offset+opP_rows,col_offset+1,col_offset+opP_cols,alpha
00526       ,Range1D::Invalid
00527       ,Teuchos::rcp(new GenPermMatrixSlice(*P)),P_release,P_trans
00528       ,NULL,Teuchos::null,BLAS_Cpp::no_trans
00529       ,Range1D()
00530       ,Teuchos::null,Teuchos::null,BLAS_Cpp::no_trans
00531       )
00532     );
00533 }
00534 
00535 void MatrixComposite::remove_matrix( matrix_list_t::iterator itr )
00536 {
00537   fully_constructed_ = false;
00538   matrix_list_.erase(itr);
00539 }
00540 
00541 void MatrixComposite::finish_construction(
00542   const VectorSpace::space_ptr_t&  space_cols
00543   ,const VectorSpace::space_ptr_t& space_rows
00544   )
00545 {
00546   TEST_FOR_EXCEPTION(
00547     !space_cols.get(), std::invalid_argument
00548     ,"MatrixComposite::finish_construction(...): Error, space_cols.get() can not be NULL" );
00549   TEST_FOR_EXCEPTION(
00550     !space_rows.get(), std::invalid_argument
00551     ,"MatrixComposite::finish_construction(...): Error, space_rows.get() can not be NULL" );
00552   TEST_FOR_EXCEPTION(
00553     space_cols->dim() != rows_, std::invalid_argument
00554     ,"MatrixComposite::finish_construction(...): Error, space_colss->dim() = " << space_cols->dim()
00555     << " != rows = " << rows_ << " where cols was passed to this->reinitialize(...)" );
00556   TEST_FOR_EXCEPTION(
00557     space_rows->dim() != cols_, std::invalid_argument
00558     ,"MatrixComposite::finish_construction(...): Error, space_rows->dim() = " << space_rows->dim()
00559     << " != cols = " << cols_ << " where cols was passed to this->reinitialize(...)" );
00560 
00561   space_rows_ = space_rows;
00562   space_cols_ = space_cols;
00563   fully_constructed_ = true;
00564 }
00565 
00566 // Member access
00567 
00568 int MatrixComposite::num_vectors() const
00569 {
00570   return vector_list_.size();
00571 }
00572 
00573 MatrixComposite::vector_list_t::iterator
00574 MatrixComposite::vectors_begin()
00575 {
00576   return vector_list_.begin();
00577 }
00578 
00579 MatrixComposite::vector_list_t::iterator
00580 MatrixComposite::vectors_end()
00581 {
00582   return vector_list_.end();
00583 }
00584 
00585 MatrixComposite::vector_list_t::const_iterator
00586 MatrixComposite::vectors_begin() const
00587 {
00588   return vector_list_.begin();
00589 }
00590 
00591 MatrixComposite::vector_list_t::const_iterator
00592 MatrixComposite::vectors_end() const
00593 {
00594   return vector_list_.end();
00595 }
00596 
00597 int MatrixComposite::num_matrices() const
00598 {
00599   return matrix_list_.size();
00600 }
00601 
00602 MatrixComposite::matrix_list_t::iterator
00603 MatrixComposite::matrices_begin()
00604 {
00605   return matrix_list_.begin();
00606 }
00607 
00608 MatrixComposite::matrix_list_t::iterator
00609 MatrixComposite::matrices_end()
00610 {
00611   return matrix_list_.end();
00612 }
00613 
00614 MatrixComposite::matrix_list_t::const_iterator
00615 MatrixComposite::matrices_begin() const
00616 {
00617   return matrix_list_.begin();
00618 }
00619 
00620 MatrixComposite::matrix_list_t::const_iterator
00621 MatrixComposite::matrices_end() const
00622 {
00623   return matrix_list_.end();
00624 }
00625 
00626 // Overridden from Matrix
00627 
00628 size_type MatrixComposite::rows() const
00629 {
00630   return fully_constructed_ ? rows_ : 0;
00631 }
00632 
00633 size_type MatrixComposite::cols() const
00634 {
00635   return fully_constructed_ ? cols_ : 0;
00636 }
00637 
00638 size_type MatrixComposite::nz() const
00639 {
00640   if(fully_constructed_) {
00641     size_type nz = 0;
00642     {for(matrix_list_t::const_iterator itr = matrix_list_.begin(); itr != matrix_list_.end(); ++itr) {
00643       if( itr->A_ != NULL )
00644         nz += itr->A_->nz();
00645       else
00646         nz += (*itr->P_).nz();
00647     }}
00648     {for(vector_list_t::const_iterator itr = vector_list_.begin(); itr != vector_list_.end(); ++itr) {
00649       nz += itr->v_->nz();
00650     }}
00651     // ToDo: Update the above code to consider DMatrixSlice and Range1D views
00652     return nz;
00653   }
00654   return 0;
00655 }
00656 
00657 // Overridden from MatrixOp
00658 
00659 const VectorSpace& MatrixComposite::space_rows() const
00660 {
00661   assert_fully_constructed();
00662   return *space_rows_;
00663 }
00664 
00665 const VectorSpace& MatrixComposite::space_cols() const
00666 {
00667   assert_fully_constructed();
00668   return *space_cols_;
00669 }
00670 
00671 MatrixOp::mat_ptr_t
00672 MatrixComposite::sub_view(const Range1D& row_rng, const Range1D& col_rng) const
00673 {
00674   assert_fully_constructed();
00675   // For this initial implementation we will just look through the list of matrices
00676   // and find an exact match.  If we don't find it we will just return the result
00677   // from the default implementation of this method.
00678 
00679   // ToDo: Implement!
00680 
00681   // Just return the default sub-view
00682   return MatrixOp::sub_view(row_rng,col_rng);
00683 }
00684 
00685 void MatrixComposite::Vp_StMtV(
00686   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00687   , const Vector& x, value_type b
00688   ) const
00689 {
00690 #ifdef PROFILE_HACK_ENABLED
00691   ProfileHackPack::ProfileTiming profile_timing( "MatrixComposite::Vp_StMtV(...DVectorSlice...)" );
00692 #endif
00693   assert_fully_constructed();
00694   AbstractLinAlgPack::Vp_MtV_assert_compatibility( y, *this, M_trans, x );
00695   Vp_StMtV_imp(y,a,rows_,cols_,matrix_list_,vector_list_,M_trans,x,b);
00696 }
00697 
00698 void MatrixComposite::Vp_StMtV(
00699   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00700   , const SpVectorSlice& x, value_type b
00701   ) const
00702 {
00703 #ifdef PROFILE_HACK_ENABLED
00704   ProfileHackPack::ProfileTiming profile_timing( "MatrixComposite::Vp_StMtV(...SpVectorSlice...)" );
00705 #endif
00706   assert_fully_constructed();
00707   AbstractLinAlgPack::Vp_MtV_assert_compatibility( y, *this, M_trans, x );
00708   Vp_StMtV_imp(y,a,rows_,cols_,matrix_list_,vector_list_,M_trans,x,b);
00709 }
00710 
00711 void MatrixComposite::Vp_StPtMtV(
00712   VectorMutable* y, value_type a
00713   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00714   , BLAS_Cpp::Transp M_trans
00715   , const Vector& x, value_type b
00716   ) const
00717 {
00718 #ifdef PROFILE_HACK_ENABLED
00719   ProfileHackPack::ProfileTiming profile_timing( "MatrixComposite::Vp_StPtMtV(...DVectorSlice...)" );
00720 #endif
00721   assert_fully_constructed();
00722   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement when needed!
00723 }
00724 
00725 void MatrixComposite::Vp_StPtMtV(
00726   VectorMutable* y, value_type a
00727   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00728   , BLAS_Cpp::Transp M_trans
00729   , const SpVectorSlice& x, value_type b
00730   ) const
00731 { 
00732 #ifdef PROFILE_HACK_ENABLED
00733   ProfileHackPack::ProfileTiming profile_timing( "MatrixComposite::Vp_StPtMtV(...SpVectorSlice...)" );
00734 #endif
00735   assert_fully_constructed();
00736   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement when needed!
00737 }
00738 
00739 // private
00740 
00741 void MatrixComposite::assert_fully_constructed() const
00742 {
00743   const bool fully_constructed = fully_constructed_;
00744   TEST_FOR_EXCEPTION(
00745     !fully_constructed, std::logic_error
00746     ,"MatrixComposite::assert_fully_constructed() : Error, not fully constructed!");
00747 }
00748 
00749 } // end namespace AbstractLinAlgPack

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