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