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