AbstractLinAlgPack_MatrixOpSubView.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 <typeinfo>
00032 #include <stdexcept>
00033 
00034 #include "AbstractLinAlgPack_MatrixOpSubView.hpp"
00035 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00036 #include "AbstractLinAlgPack_VectorSpace.hpp"
00037 #include "AbstractLinAlgPack_VectorMutable.hpp"
00038 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00039 #include "AbstractLinAlgPack_SpVectorView.hpp"
00040 #include "AbstractLinAlgPack_EtaVector.hpp"
00041 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00042 #include "Teuchos_RCP.hpp"
00043 #include "Teuchos_TestForException.hpp"
00044 
00045 namespace AbstractLinAlgPack {
00046 
00047 MatrixOpSubView::MatrixOpSubView(
00048   const mat_ptr_t& M_full
00049   ,const Range1D& rng_rows
00050   ,const Range1D& rng_cols
00051   ,BLAS_Cpp::Transp M_trans
00052   )
00053 {
00054   this->initialize(M_full,rng_rows,rng_cols,M_trans);
00055 }
00056     
00057 void MatrixOpSubView::initialize(
00058   const mat_ptr_t& M_full
00059   ,const Range1D& rng_rows_in
00060   ,const Range1D& rng_cols_in
00061   ,BLAS_Cpp::Transp M_trans
00062   )
00063 {
00064   namespace rcp = MemMngPack;
00065 
00066   if( M_full.get() ) {
00067     const index_type
00068       M_rows = M_full->rows(),
00069       M_cols = M_full->cols();
00070     const Range1D
00071       rng_rows = RangePack::full_range(rng_rows_in,1,M_rows),
00072       rng_cols = RangePack::full_range(rng_cols_in,1,M_cols);
00073     TEST_FOR_EXCEPTION(
00074       rng_rows.ubound() > M_rows, std::invalid_argument
00075       ,"MatrixOpSubView::initialize(...): Error, "
00076       "rng_rows = ["<<rng_rows.lbound()<<","<<rng_rows.ubound()<<"] is of range of "
00077       "[1,M_full->rows()] = [1,"<<M_rows<<"]" );
00078     TEST_FOR_EXCEPTION(
00079       rng_cols.ubound() > M_cols, std::invalid_argument
00080       ,"MatrixOpSubView::initialize(...): Error, "
00081       "rng_cols = ["<<rng_cols.lbound()<<","<<rng_cols.ubound()<<"] is of range of "
00082       "[1,M_full->cols()] = [1,"<<M_cols<<"]" );
00083     M_full_     = M_full;
00084     rng_rows_   = rng_rows;
00085     rng_cols_   = rng_cols;
00086     M_trans_    = M_trans;
00087     space_cols_ = ( M_trans == BLAS_Cpp::no_trans
00088             ? M_full->space_cols().sub_space(rng_rows)->clone()
00089             : M_full->space_rows().sub_space(rng_cols)->clone() );
00090     space_rows_ = ( M_trans == BLAS_Cpp::no_trans
00091             ? M_full->space_rows().sub_space(rng_cols)->clone()
00092             : M_full->space_cols().sub_space(rng_rows)->clone() );
00093   }
00094   else {
00095     M_full_     = Teuchos::null;
00096     rng_rows_   = Range1D::Invalid;
00097     rng_cols_   = Range1D::Invalid;
00098     M_trans_    = BLAS_Cpp::no_trans;
00099     space_cols_ = Teuchos::null;
00100     space_rows_ = Teuchos::null;
00101   }
00102 }
00103 
00104 // overridden from MatrixBase
00105 
00106 size_type MatrixOpSubView::rows() const
00107 {
00108   return ( M_full_.get() 
00109        ? BLAS_Cpp::rows( rng_rows_.size(), rng_cols_.size(), M_trans_ )
00110        : 0 );
00111 }
00112 
00113 size_type MatrixOpSubView::cols() const
00114 {
00115   return ( M_full_.get() 
00116        ? BLAS_Cpp::cols( rng_rows_.size(), rng_cols_.size(), M_trans_ )
00117        : 0 );
00118 }
00119 
00120 size_type MatrixOpSubView::nz() const
00121 {
00122   return ( M_full_.get()
00123        ? ( rng_rows_.full_range() && rng_cols_.full_range()
00124          ? M_full_->nz()
00125          : MatrixBase::nz() )
00126        : 0 );
00127 }
00128 
00129 // Overridden form MatrixOp
00130 
00131 const VectorSpace& MatrixOpSubView::space_cols() const
00132 {
00133   assert_initialized();
00134   return *space_cols_;
00135 }
00136 
00137 const VectorSpace& MatrixOpSubView::space_rows() const
00138 {
00139   assert_initialized();
00140   return *space_rows_;
00141 }
00142 
00143 MatrixOp::mat_ptr_t
00144 MatrixOpSubView::sub_view(const Range1D& row_rng, const Range1D& col_rng) const
00145 {
00146   assert_initialized();
00147   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00148   return Teuchos::null;
00149 }
00150 
00151 void MatrixOpSubView::zero_out()
00152 {
00153   assert_initialized();
00154   if( rng_rows_.full_range() && rng_cols_.full_range() ) {
00155     M_full_->zero_out();
00156     return;
00157   }
00158   TEST_FOR_EXCEPTION(
00159     true, std::logic_error, "MatrixOpSubView::zero_out(): "
00160     "Error, this method can not be implemented with a sub-view" );
00161 }
00162 
00163 void MatrixOpSubView::Mt_S( value_type alpha )
00164 {
00165   assert_initialized();
00166   if( rng_rows_.full_range() && rng_cols_.full_range() ) {
00167     M_full_->Mt_S(alpha);
00168     return;
00169   }
00170   TEST_FOR_EXCEPTION(
00171     true, std::logic_error, "MatrixOpSubView::Mt_S(alpha): "
00172     "Error, this method can not be implemented with a sub-view" );
00173 }
00174 
00175 MatrixOp& MatrixOpSubView::operator=(const MatrixOp& M)
00176 {
00177   assert_initialized();
00178   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00179   return *this;
00180 }
00181 
00182 std::ostream& MatrixOpSubView::output(std::ostream& out) const
00183 {
00184   assert_initialized();
00185   return MatrixOp::output(out); // ToDo: Specialize if needed?
00186 }
00187 
00188 // Level-1 BLAS
00189 
00190 // rhs matrix argument
00191 
00192 bool MatrixOpSubView::Mp_StM(
00193   MatrixOp* m_lhs, value_type alpha
00194   , BLAS_Cpp::Transp trans_rhs) const
00195 {
00196   assert_initialized();
00197   return MatrixOp::Mp_StM(m_lhs,alpha,trans_rhs); // ToDo: Specialize?
00198 }
00199 
00200 bool MatrixOpSubView::Mp_StMtP(
00201   MatrixOp* m_lhs, value_type alpha
00202   , BLAS_Cpp::Transp M_trans
00203   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00204   ) const
00205 {
00206   assert_initialized();
00207   return MatrixOp::Mp_StMtP(m_lhs,alpha,M_trans,P_rhs,P_rhs_trans); // ToDo: Specialize?
00208 }
00209 
00210 bool MatrixOpSubView::Mp_StPtM(
00211   MatrixOp* m_lhs, value_type alpha
00212   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00213   , BLAS_Cpp::Transp M_trans
00214   ) const
00215 {
00216   assert_initialized();
00217   return MatrixOp::Mp_StPtM(m_lhs,alpha,P_rhs,P_rhs_trans,M_trans); // ToDo: Specialize?
00218 }
00219 
00220 bool MatrixOpSubView::Mp_StPtMtP(
00221   MatrixOp* m_lhs, value_type alpha
00222   , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00223   , BLAS_Cpp::Transp M_trans
00224   , const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
00225   ) const
00226 {
00227   assert_initialized();
00228   return MatrixOp::Mp_StPtMtP(
00229     m_lhs,alpha,P_rhs1,P_rhs1_trans,M_trans,P_rhs2,P_rhs2_trans); // ToDo: Specialize?
00230 }
00231 
00232 // lhs matrix argument
00233 
00234 bool MatrixOpSubView::Mp_StM(
00235   value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs)
00236 {
00237   assert_initialized();
00238   return MatrixOp::Mp_StM(alpha,M_rhs,trans_rhs); // ToDo: Specialize?
00239 }
00240 
00241 bool MatrixOpSubView::Mp_StMtP(
00242   value_type alpha
00243   ,const MatrixOp& M_rhs, BLAS_Cpp::Transp M_trans
00244   ,const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00245   )
00246 {
00247   assert_initialized();
00248   return MatrixOp::Mp_StMtP(alpha,M_rhs,M_trans,P_rhs,P_rhs_trans); // ToDo: Specialize?
00249 }
00250 
00251 bool MatrixOpSubView::Mp_StPtM(
00252   value_type alpha
00253   ,const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00254   ,const MatrixOp& M_rhs, BLAS_Cpp::Transp M_trans
00255   )
00256 {
00257   assert_initialized();
00258   return MatrixOp::Mp_StPtM(
00259     alpha,P_rhs,P_rhs_trans,M_rhs,M_trans); // ToDo: Specialize?
00260 }
00261 
00262 bool MatrixOpSubView::Mp_StPtMtP(
00263   value_type alpha
00264   ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00265   ,const MatrixOp& M_rhs, BLAS_Cpp::Transp M_trans
00266   ,const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
00267   )
00268 {
00269   assert_initialized();
00270   return MatrixOp::Mp_StPtMtP(
00271     alpha,P_rhs1,P_rhs1_trans,M_rhs,M_trans,P_rhs2,P_rhs2_trans); // ToDo: Specialize?
00272 }
00273 
00274 // Level-2 BLAS
00275 
00276 void MatrixOpSubView::Vp_StMtV(
00277   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans_in
00278   , const Vector& x, value_type b
00279   ) const
00280 {
00281   using BLAS_Cpp::no_trans;
00282   using BLAS_Cpp::trans;
00283 
00284   assert_initialized();
00285 
00286   BLAS_Cpp::Transp
00287     M_trans_trans = ( M_trans_==no_trans ? M_trans_in : BLAS_Cpp::trans_not(M_trans_in) );
00288 
00289   // ToDo: Assert compatibility!
00290 
00291   if( rng_rows_.full_range() && rng_cols_.full_range() ) {
00292     AbstractLinAlgPack::Vp_StMtV(  // The matrix is just transposed
00293       y, a
00294       ,*M_full_, M_trans_trans
00295       ,x, b );
00296     return;
00297   }
00298   // y = b*y
00299   Vt_S( y, b );
00300   //
00301   // xt1                      = 0.0
00302   // xt3 = xt(op_op_rng_cols) = x
00303   // xt3                      = 0.0
00304   //
00305   // [ yt1 ]                        [ xt1 ]
00306   // [ yt2 ] = a * op(op(M_full)) * [ xt3 ]
00307   // [ yt3 ]                        [ xt3 ]
00308   //
00309   // =>
00310   //
00311   // y += yt2 = yt(op_op_rng_rows)
00312   //
00313   const Range1D
00314     op_op_rng_rows = ( M_trans_trans == no_trans ? rng_rows_ : rng_cols_ ),
00315     op_op_rng_cols = ( M_trans_trans == no_trans ? rng_cols_ : rng_rows_ );
00316   VectorSpace::vec_mut_ptr_t
00317     xt = ( M_trans_trans == no_trans ? M_full_->space_rows() : M_full_->space_cols() ).create_member(),
00318     yt = ( M_trans_trans == no_trans ? M_full_->space_cols() : M_full_->space_rows() ).create_member();
00319   *xt = 0.0;
00320   *xt->sub_view(op_op_rng_cols) = x;
00321     LinAlgOpPack::V_StMtV( yt.get(), a, *M_full_, M_trans_trans, *xt );
00322   LinAlgOpPack::Vp_V( y, *yt->sub_view(op_op_rng_rows) );
00323 }
00324 
00325 void MatrixOpSubView::Vp_StMtV(
00326   VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00327   , const SpVectorSlice& sv_rhs2, value_type beta) const
00328 {
00329   assert_initialized();
00330   MatrixOp::Vp_StMtV(v_lhs,alpha,trans_rhs1,sv_rhs2,beta); // ToDo: Specialize?
00331 }
00332 
00333 void MatrixOpSubView::Vp_StPtMtV(
00334   VectorMutable* v_lhs, value_type alpha
00335   , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00336   , BLAS_Cpp::Transp M_rhs2_trans
00337   , const Vector& v_rhs3, value_type beta) const
00338 {
00339   assert_initialized();
00340   MatrixOp::Vp_StPtMtV(
00341     v_lhs,alpha,P_rhs1,P_rhs1_trans,M_rhs2_trans,v_rhs3,beta); // ToDo: Specialize?
00342 }
00343 
00344 void MatrixOpSubView::Vp_StPtMtV(
00345   VectorMutable* v_lhs, value_type alpha
00346   , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00347   , BLAS_Cpp::Transp M_rhs2_trans
00348   , const SpVectorSlice& sv_rhs3, value_type beta) const
00349 {
00350   assert_initialized();
00351   MatrixOp::Vp_StPtMtV(
00352     v_lhs,alpha,P_rhs1,P_rhs1_trans,M_rhs2_trans,sv_rhs3,beta); // ToDo: Specialize?
00353 }
00354 
00355 value_type MatrixOpSubView::transVtMtV(
00356   const Vector& v_rhs1, BLAS_Cpp::Transp trans_rhs2
00357   , const Vector& v_rhs3) const
00358 {
00359   assert_initialized();
00360   return MatrixOp::transVtMtV(v_rhs1,trans_rhs2,v_rhs3); // ToDo: Specialize?
00361 }
00362 
00363 value_type MatrixOpSubView::transVtMtV(
00364   const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2
00365   , const SpVectorSlice& sv_rhs3) const
00366 {
00367   assert_initialized();
00368   return MatrixOp::transVtMtV(sv_rhs1,trans_rhs2,sv_rhs3); // ToDo: Specialize?
00369 }
00370 
00371 void MatrixOpSubView::syr2k(
00372   BLAS_Cpp::Transp M_trans, value_type alpha
00373   , const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
00374   , const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
00375   , value_type beta, MatrixSymOp* sym_lhs ) const
00376 {
00377   assert_initialized();
00378   MatrixOp::syr2k(
00379     M_trans,alpha,P1,P1_trans,P2,P2_trans,beta,sym_lhs); // ToDo: Specialize?
00380 }
00381 
00382 // Level-3 BLAS
00383 
00384 bool MatrixOpSubView::Mp_StMtM(
00385   MatrixOp* m_lhs, value_type alpha
00386   , BLAS_Cpp::Transp trans_rhs1, const MatrixOp& mwo_rhs2
00387   , BLAS_Cpp::Transp trans_rhs2, value_type beta) const
00388 {
00389   assert_initialized();
00390   return MatrixOp::Mp_StMtM(
00391     m_lhs,alpha,trans_rhs1,mwo_rhs2,trans_rhs2,beta); // ToDo: Specialize?
00392 }
00393 
00394 bool MatrixOpSubView::Mp_StMtM(
00395   MatrixOp* m_lhs, value_type alpha
00396   , const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00397   , BLAS_Cpp::Transp trans_rhs2, value_type beta ) const
00398 {
00399   return MatrixOp::Mp_StMtM(
00400     m_lhs,alpha,mwo_rhs1,trans_rhs1,trans_rhs2,beta); // ToDo: Specialize?
00401 }
00402 
00403 bool MatrixOpSubView::Mp_StMtM(
00404   value_type alpha
00405   ,const MatrixOp& mvw_rhs1, BLAS_Cpp::Transp trans_rhs1
00406   ,const MatrixOp& mwo_rhs2,BLAS_Cpp::Transp trans_rhs2
00407   ,value_type beta )
00408 {
00409   assert_initialized();
00410   return MatrixOp::Mp_StMtM(
00411     alpha,mvw_rhs1,trans_rhs1,mwo_rhs2,trans_rhs2,beta); // ToDo: Specialize?
00412 }
00413 
00414 bool MatrixOpSubView::syrk(
00415   BLAS_Cpp::Transp M_trans, value_type alpha
00416   , value_type beta, MatrixSymOp* sym_lhs ) const
00417 {
00418   assert_initialized();
00419   return MatrixOp::syrk(M_trans,alpha,beta,sym_lhs); // ToDo: Specialize?
00420 }
00421 
00422 // private
00423 
00424 void MatrixOpSubView::assert_initialized() const {
00425   TEST_FOR_EXCEPTION(
00426     M_full_.get() == NULL, std::logic_error
00427     ,"Error, the MatrixOpSubView object has not been initialize!" );
00428 }
00429 
00430 } // end namespace AbstractLinAlgPack

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