MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_MatrixOp.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 #include <math.h>
00031 
00032 #include <typeinfo>
00033 #include <stdexcept>
00034 
00035 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
00036 #include "AbstractLinAlgPack_MatrixOpSubView.hpp"
00037 #include "AbstractLinAlgPack_MatrixPermAggr.hpp"
00038 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00039 #include "AbstractLinAlgPack_VectorSpace.hpp"
00040 #include "AbstractLinAlgPack_VectorMutable.hpp"
00041 #include "AbstractLinAlgPack_Permutation.hpp"
00042 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00043 #include "AbstractLinAlgPack_SpVectorView.hpp"
00044 #include "AbstractLinAlgPack_EtaVector.hpp"
00045 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00046 #include "Teuchos_TestForException.hpp"
00047 #include "Teuchos_FancyOStream.hpp"
00048 
00049 namespace AbstractLinAlgPack {
00050 
00051 void MatrixOp::zero_out()
00052 {
00053   TEST_FOR_EXCEPTION(
00054     true, std::logic_error, "MatrixOp::zero_out(): "
00055     "Error, this method as not been defined by the subclass \'"
00056     <<typeName(*this)<<"\'" );
00057 }
00058 
00059 void MatrixOp::Mt_S(value_type alpha)
00060 {
00061   TEST_FOR_EXCEPTION(
00062     true, std::logic_error, "MatrixOp::Mt_S(): "
00063     "Error, this method as not been defined by the subclass \'"
00064     <<typeName(*this)<<"\'" );
00065 }
00066 
00067 MatrixOp& MatrixOp::operator=(const MatrixOp& M)
00068 {
00069   const bool assign_to_self = dynamic_cast<const void*>(this) == dynamic_cast<const void*>(&M);
00070   TEST_FOR_EXCEPTION(
00071     !assign_to_self, std::logic_error
00072     ,"MatrixOp::operator=(M) : Error, this is not assignment "
00073     "to self and this method is not overridden for the subclass \'"
00074     << typeName(*this) << "\'" );
00075   return *this; // assignment to self
00076 }
00077 
00078 std::ostream& MatrixOp::output(std::ostream& out_arg) const
00079 {
00080   Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
00081   Teuchos::OSTab tab(out);
00082   const size_type m = this->rows(), n = this->cols();
00083   VectorSpace::vec_mut_ptr_t
00084     row_vec = space_rows().create_member(); // dim() == n
00085   *out << m << " " << n << std::endl;
00086   for( size_type i = 1; i <= m; ++i ) {
00087     LinAlgOpPack::V_StMtV( &*row_vec, 1.0, *this, BLAS_Cpp::trans, EtaVector(i,m)() );
00088     row_vec->output(*out,false,true);
00089   }
00090   return out_arg;
00091 }
00092 
00093 // Clone
00094 
00095 MatrixOp::mat_mut_ptr_t
00096 MatrixOp::clone()
00097 {
00098   return Teuchos::null;
00099 }
00100 
00101 MatrixOp::mat_ptr_t
00102 MatrixOp::clone() const
00103 {
00104   return Teuchos::null;
00105 }
00106 
00107 // Norms
00108 
00109 const MatrixOp::MatNorm
00110 MatrixOp::calc_norm(
00111   EMatNormType  requested_norm_type
00112   ,bool         allow_replacement
00113   ) const
00114 {
00115   using BLAS_Cpp::no_trans;
00116   using BLAS_Cpp::trans;
00117   using LinAlgOpPack::V_MtV;
00118   const VectorSpace
00119     &space_cols = this->space_cols(),
00120     &space_rows = this->space_rows();
00121   const index_type
00122     num_rows = space_cols.dim(),
00123     num_cols = space_rows.dim();
00124   TEST_FOR_EXCEPTION(
00125     !(requested_norm_type == MAT_NORM_1 || requested_norm_type == MAT_NORM_INF), MethodNotImplemented
00126     ,"MatrixOp::calc_norm(...): Error, This default implemenation can only "
00127     "compute the one norm or the infinity norm!"
00128     );
00129   //
00130   // Here we implement Algorithm 2.5 in "Applied Numerical Linear Algebra", Demmel (1997)
00131   // using the momenclature in the text.
00132   //
00133   const MatrixOp
00134     &B = *this;
00135   bool
00136     do_trans = requested_norm_type == MAT_NORM_INF;
00137   VectorSpace::vec_mut_ptr_t
00138     x    = (!do_trans ? space_rows : space_cols).create_member(1.0/(!do_trans ? num_cols : num_rows)),
00139     w    = (!do_trans ? space_cols : space_rows).create_member(),
00140     zeta = (!do_trans ? space_cols : space_rows).create_member(),
00141     z    = (!do_trans ? space_rows : space_cols).create_member();
00142   const index_type max_iter = 5; // Recommended by Highman 1988, (see Demmel's reference)
00143   value_type w_nrm = 0.0;
00144   for( index_type k = 0; k <= max_iter; ++k ) {
00145     V_MtV( w.get(), B, !do_trans ? no_trans : trans, *x );    // w = B*x
00146     sign( *w, zeta.get() );                                   // zeta = sign(w)
00147     V_MtV( z.get(), B, !do_trans ? trans : no_trans, *zeta ); // z = B'*zeta
00148     value_type  z_j = 0.0;                                    // max |z(j)| = ||z||inf
00149     index_type  j   = 0;
00150     max_abs_ele( *z, &z_j, &j );
00151     const value_type zTx = dot(*z,*x);                        // z'*x
00152     w_nrm = w->norm_1();                                      // ||w||1
00153     if( ::fabs(z_j) <= zTx ) {                                // Update
00154       break;
00155     }
00156     else {
00157       *x = 0.0;
00158       x->set_ele(j,1.0);
00159     }
00160   }
00161   return MatNorm(w_nrm,requested_norm_type);
00162 }
00163 
00164 // Subview
00165 
00166 MatrixOp::mat_ptr_t
00167 MatrixOp::sub_view(const Range1D& row_rng, const Range1D& col_rng) const
00168 {
00169   namespace rcp = MemMngPack;
00170 
00171   if( 
00172     ( ( row_rng.lbound() == 1 && row_rng.ubound() == this->rows() )
00173       || row_rng.full_range() )
00174     &&
00175     ( ( col_rng.lbound() == 1 && col_rng.ubound() == this->cols() )
00176       || row_rng.full_range() )
00177     ) 
00178   {
00179     return Teuchos::rcp(this,false); // don't clean up memory
00180   }
00181   return Teuchos::rcp(
00182     new MatrixOpSubView(
00183       Teuchos::rcp(const_cast<MatrixOp*>(this),false) // don't clean up memory
00184       ,row_rng,col_rng ) );
00185 }
00186 
00187 // Permuted views
00188 
00189 MatrixOp::mat_ptr_t
00190 MatrixOp::perm_view(
00191   const Permutation          *P_row
00192   ,const index_type          row_part[]
00193   ,int                       num_row_part
00194   ,const Permutation         *P_col
00195   ,const index_type          col_part[]
00196   ,int                       num_col_part
00197   ) const
00198 {
00199   namespace rcp = MemMngPack;
00200   return Teuchos::rcp(
00201     new MatrixPermAggr(
00202       Teuchos::rcp(this,false)
00203       ,Teuchos::rcp(P_row,false)
00204       ,Teuchos::rcp(P_col,false)
00205       ,Teuchos::null
00206       ) );
00207 }
00208 
00209 MatrixOp::mat_ptr_t
00210 MatrixOp::perm_view_update(
00211   const Permutation          *P_row
00212   ,const index_type          row_part[]
00213   ,int                       num_row_part
00214   ,const Permutation         *P_col
00215   ,const index_type          col_part[]
00216   ,int                       num_col_part
00217   ,const mat_ptr_t           &perm_view
00218   ) const
00219 {
00220   return this->perm_view(
00221     P_row,row_part,num_row_part
00222     ,P_col,col_part,num_col_part );
00223 }
00224 
00225 // Level-1 BLAS
00226 
00227 bool MatrixOp::Mp_StM(
00228   MatrixOp* m_lhs, value_type alpha
00229   , BLAS_Cpp::Transp trans_rhs) const
00230 {
00231   return false;
00232 }
00233 
00234 bool MatrixOp::Mp_StM(
00235   value_type alpha,const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs)
00236 {
00237   return false;
00238 }
00239 
00240 bool MatrixOp::Mp_StMtP(
00241   MatrixOp* m_lhs, value_type alpha
00242   , BLAS_Cpp::Transp M_trans
00243   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00244   ) const
00245 {
00246   return false;
00247 }
00248 
00249 bool MatrixOp::Mp_StMtP(
00250   value_type alpha
00251   ,const MatrixOp& M_rhs, BLAS_Cpp::Transp M_trans
00252   ,const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00253   )
00254 {
00255   return false;
00256 }
00257 
00258 bool MatrixOp::Mp_StPtM(
00259   MatrixOp* m_lhs, value_type alpha
00260   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00261   , BLAS_Cpp::Transp M_trans
00262   ) const
00263 {
00264   return false;
00265 }
00266 
00267 bool MatrixOp::Mp_StPtM(
00268   value_type alpha
00269   ,const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00270   ,const MatrixOp& M_rhs, BLAS_Cpp::Transp M_trans
00271   )
00272 {
00273   return false;
00274 }
00275 
00276 bool MatrixOp::Mp_StPtMtP(
00277   MatrixOp* m_lhs, value_type alpha
00278   , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00279   , BLAS_Cpp::Transp M_trans
00280   , const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
00281   ) const
00282 {
00283   return false;
00284 }
00285 
00286 bool MatrixOp::Mp_StPtMtP(
00287   value_type alpha
00288   ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00289   ,const MatrixOp& M_rhs, BLAS_Cpp::Transp M_trans
00290   ,const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
00291   )
00292 {
00293   return false;
00294 }
00295 
00296 // Level-2 BLAS
00297 
00298 void MatrixOp::Vp_StMtV(
00299   VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00300   , const SpVectorSlice& sv_rhs2, value_type beta) const
00301 {
00302   Vp_MtV_assert_compatibility(v_lhs,*this,trans_rhs1,sv_rhs2 );
00303   if( sv_rhs2.nz() ) {
00304     VectorSpace::vec_mut_ptr_t
00305       v_rhs2 = (trans_rhs1 == BLAS_Cpp::no_trans
00306             ? this->space_rows()
00307             : this->space_cols()
00308         ).create_member();
00309     v_rhs2->set_sub_vector(sub_vec_view(sv_rhs2));
00310     this->Vp_StMtV(v_lhs,alpha,trans_rhs1,*v_rhs2,beta);
00311   }
00312   else {
00313     Vt_S( v_lhs, beta );
00314   }
00315 }
00316 
00317 void MatrixOp::Vp_StPtMtV(
00318   VectorMutable* y, value_type a
00319   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00320   ,BLAS_Cpp::Transp M_trans
00321   ,const Vector& x, value_type b
00322   ) const
00323 {
00324   VectorSpace::vec_mut_ptr_t
00325     t = ( M_trans == BLAS_Cpp::no_trans
00326         ? this->space_cols()
00327         : this->space_rows() ).create_member();
00328   LinAlgOpPack::V_MtV( t.get(), *this, M_trans, x );
00329   AbstractLinAlgPack::Vp_StMtV( y, a, P, P_trans, *t, b );
00330 }
00331 
00332 void MatrixOp::Vp_StPtMtV(
00333   VectorMutable* y, value_type a
00334   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00335   ,BLAS_Cpp::Transp M_trans
00336   ,const SpVectorSlice& x, value_type b
00337   ) const
00338 {
00339   VectorSpace::vec_mut_ptr_t
00340     t = ( M_trans == BLAS_Cpp::no_trans
00341         ? this->space_cols()
00342         : this->space_rows() ).create_member();
00343   LinAlgOpPack::V_MtV( t.get(), *this, M_trans, x );
00344   AbstractLinAlgPack::Vp_StMtV( y, a, P, P_trans, *t, b );
00345 }
00346 
00347 value_type MatrixOp::transVtMtV(
00348   const Vector& vs_rhs1, BLAS_Cpp::Transp trans_rhs2
00349   , const Vector& vs_rhs3) const
00350 {
00351   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00352   return 0.0;
00353 }
00354 
00355 value_type MatrixOp::transVtMtV(
00356   const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2
00357   , const SpVectorSlice& sv_rhs3) const
00358 {
00359   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00360   return 0.0;
00361 }
00362 
00363 void MatrixOp::syr2k(
00364   BLAS_Cpp::Transp M_trans, value_type alpha
00365   , const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
00366   , const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
00367   , value_type beta, MatrixSymOp* sym_lhs ) const
00368 {
00369   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00370 }
00371 
00372 // Level-3 BLAS
00373 
00374 bool MatrixOp::Mp_StMtM(
00375   MatrixOp* C, value_type a
00376   ,BLAS_Cpp::Transp A_trans, const MatrixOp& B
00377   ,BLAS_Cpp::Transp B_trans, value_type b) const
00378 {
00379   return false;
00380 }
00381 
00382 bool MatrixOp::Mp_StMtM(
00383   MatrixOp* m_lhs, value_type alpha
00384   ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00385   ,BLAS_Cpp::Transp trans_rhs2, value_type beta ) const
00386 {
00387   return false;
00388 }
00389 
00390 bool MatrixOp::Mp_StMtM(
00391   value_type alpha
00392   ,const MatrixOp& mvw_rhs1, BLAS_Cpp::Transp trans_rhs1
00393   ,const MatrixOp& mwo_rhs2,BLAS_Cpp::Transp trans_rhs2
00394   ,value_type beta )
00395 {
00396   return false;
00397 }
00398 
00399 bool MatrixOp::syrk(
00400   BLAS_Cpp::Transp   M_trans
00401   ,value_type        alpha
00402   ,value_type        beta
00403   ,MatrixSymOp   *sym_lhs
00404   ) const
00405 {
00406   return false;
00407 }
00408 
00409 bool MatrixOp::syrk(
00410   const MatrixOp  &mwo_rhs
00411   ,BLAS_Cpp::Transp   M_trans
00412   ,value_type         alpha
00413   ,value_type         beta
00414   )
00415 {
00416   return false;
00417 }
00418 
00419 } // end namespace AbstractLinAlgPack
00420 
00421 // Non-member functions
00422 
00423 // level-1 BLAS
00424 
00425 void AbstractLinAlgPack::Mt_S( MatrixOp* m_lhs, value_type alpha )
00426 {
00427   if(alpha == 0.0)
00428     m_lhs->zero_out();
00429   else if( alpha != 1.0 )
00430     m_lhs->Mt_S(alpha);
00431 }
00432 
00433 void AbstractLinAlgPack::Mp_StM(
00434   MatrixOp* mwo_lhs, value_type alpha, const MatrixOp& M_rhs
00435   , BLAS_Cpp::Transp trans_rhs)
00436 {
00437   using BLAS_Cpp::no_trans;
00438   using BLAS_Cpp::trans;
00439 
00440   // Give the rhs argument a chance to implement the operation
00441   if(M_rhs.Mp_StM(mwo_lhs,alpha,trans_rhs))
00442     return;
00443 
00444   // Give the lhs argument a change to implement the operation
00445   if(mwo_lhs->Mp_StM(alpha,M_rhs,trans_rhs))
00446     return;
00447 
00448   // We must try to implement the method
00449   MultiVectorMutable
00450     *m_mut_lhs = dynamic_cast<MultiVectorMutable*>(mwo_lhs);
00451   TEST_FOR_EXCEPTION(
00452     !m_mut_lhs || !(m_mut_lhs->access_by() & MultiVector::COL_ACCESS)
00453     ,MatrixOp::MethodNotImplemented
00454     ,"MatrixOp::Mp_StM(...) : Error, mwo_lhs of type \'"
00455     << typeName(*mwo_lhs) << "\' could not implement the operation "
00456     "and does not support the "
00457     "\'MultiVectorMutable\' interface.  Furthermore, "
00458     "the rhs matix argument of type \'" << typeName(*mwo_lhs)
00459     << "\' could not implement the operation!" );
00460     
00461 #ifdef TEUCHOS_DEBUG
00462   TEST_FOR_EXCEPTION(
00463     !mwo_lhs->space_rows().is_compatible(
00464       trans_rhs == no_trans ? M_rhs.space_rows() : M_rhs.space_cols() )
00465     || !mwo_lhs->space_cols().is_compatible(
00466       trans_rhs == no_trans ? M_rhs.space_cols() : M_rhs.space_rows() )
00467     , MatrixOp::IncompatibleMatrices
00468     ,"MatrixOp::Mp_StM(mwo_lhs,...): Error, mwo_lhs of type \'"<<typeName(*mwo_lhs)<<"\' "
00469     <<"is not compatible with M_rhs of type \'"<<typeName(M_rhs)<<"\'" );
00470 #endif
00471 
00472   const size_type
00473     rows = BLAS_Cpp::rows( mwo_lhs->rows(), mwo_lhs->cols(), trans_rhs ),
00474     cols = BLAS_Cpp::cols( mwo_lhs->rows(), mwo_lhs->cols(), trans_rhs );
00475   for( size_type j = 1; j <= cols; ++j )
00476     AbstractLinAlgPack::Vp_StMtV( m_mut_lhs->col(j).get(), alpha, M_rhs, trans_rhs, EtaVector(j,cols)() );
00477   // ToDo: consider row access?
00478 }
00479 
00480 void AbstractLinAlgPack::Mp_StMtP(
00481   MatrixOp* mwo_lhs, value_type alpha
00482   , const MatrixOp& M_rhs, BLAS_Cpp::Transp M_trans
00483   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00484   )
00485 {
00486 
00487   // Give the rhs argument a chance to implement the operation
00488   if(M_rhs.Mp_StMtP(mwo_lhs,alpha,M_trans,P_rhs,P_rhs_trans))
00489     return;
00490 
00491   // Give the lhs argument a change to implement the operation
00492   if(mwo_lhs->Mp_StMtP(alpha,M_rhs,M_trans,P_rhs,P_rhs_trans))
00493     return;
00494 
00495   // We must try to implement the method
00496   MultiVectorMutable
00497     *m_mut_lhs = dynamic_cast<MultiVectorMutable*>(mwo_lhs);
00498   TEST_FOR_EXCEPTION(
00499     !m_mut_lhs, MatrixOp::MethodNotImplemented
00500     ,"MatrixOp::Mp_StMtP(...) : Error, mwo_lhs of type \'"
00501     << typeName(*mwo_lhs) << "\' does not support the "
00502     "\'MultiVectorMutable\' interface!" );
00503 
00504   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00505 }
00506 
00507 void AbstractLinAlgPack::Mp_StPtM(
00508   MatrixOp* mwo_lhs, value_type alpha
00509   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00510   , const MatrixOp& M_rhs, BLAS_Cpp::Transp M_trans
00511   )
00512 {
00513 
00514   // Give the rhs argument a chance to implement the operation
00515   if(M_rhs.Mp_StPtM(mwo_lhs,alpha,P_rhs,P_rhs_trans,M_trans))
00516     return;
00517 
00518   // Give the lhs argument a change to implement the operation
00519   if(mwo_lhs->Mp_StPtM(alpha,P_rhs,P_rhs_trans,M_rhs,M_trans))
00520     return;
00521 
00522   // We must try to implement the method
00523   MultiVectorMutable
00524     *m_mut_lhs = dynamic_cast<MultiVectorMutable*>(mwo_lhs);
00525   TEST_FOR_EXCEPTION(
00526     !m_mut_lhs, MatrixOp::MethodNotImplemented
00527     ,"MatrixOp::Mp_StPtM(...) : Error, mwo_lhs of type \'"
00528     << typeName(*mwo_lhs) << "\' does not support the "
00529     "\'MultiVectorMutable\' interface!" );
00530 
00531   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00532 
00533 }
00534 
00535 void AbstractLinAlgPack::Mp_StPtMtP(
00536   MatrixOp* mwo_lhs, value_type alpha
00537   , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00538   , const MatrixOp& M_rhs, BLAS_Cpp::Transp trans_rhs
00539   , const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
00540   )
00541 {
00542 
00543   // Give the rhs argument a chance to implement the operation
00544   if(M_rhs.Mp_StPtMtP(mwo_lhs,alpha,P_rhs1,P_rhs1_trans,trans_rhs,P_rhs2,P_rhs2_trans))
00545     return;
00546 
00547   // Give the lhs argument a change to implement the operation
00548   if(mwo_lhs->Mp_StPtMtP(alpha,P_rhs1,P_rhs1_trans,M_rhs,trans_rhs,P_rhs2,P_rhs2_trans))
00549     return;
00550 
00551   // We must try to implement the method
00552   MultiVectorMutable
00553     *m_mut_lhs = dynamic_cast<MultiVectorMutable*>(mwo_lhs);
00554   TEST_FOR_EXCEPTION(
00555     !m_mut_lhs, MatrixOp::MethodNotImplemented
00556     ,"MatrixOp::Mp_StPtMtP(...) : Error, mwo_lhs of type \'"
00557     << typeName(*mwo_lhs) << "\' does not support the "
00558     "\'MultiVectorMutable\' interface!" );
00559 
00560   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00561 
00562 }
00563 
00564 // level-3 blas
00565 
00566 void AbstractLinAlgPack::Mp_StMtM(
00567   MatrixOp* C, value_type a
00568   ,const MatrixOp& A, BLAS_Cpp::Transp A_trans
00569   ,const MatrixOp& B, BLAS_Cpp::Transp B_trans
00570   ,value_type b
00571   )
00572 {
00573   
00574   // Give A a chance
00575   if(A.Mp_StMtM(C,a,A_trans,B,B_trans,b))
00576     return;
00577   // Give B a chance
00578   if(B.Mp_StMtM(C,a,A,A_trans,B_trans,b))
00579     return;
00580   // Give C a chance
00581   if(C->Mp_StMtM(a,A,A_trans,B,B_trans,b))
00582     return;
00583 
00584   //
00585   // C = b*C + a*op(A)*op(B)
00586   //
00587   // We will perform this by column as:
00588   //
00589   //   C(:,j) = b*C(:,j) + a*op(A)*op(B)*e(j), for j = 1...C.cols()
00590   //
00591   //   by performing:
00592   //
00593   //        t = op(B)*e(j)
00594   //        C(:,j) = b*C(:,j) + a*op(A)*t
00595   //
00596   Mp_MtM_assert_compatibility(C,BLAS_Cpp::no_trans,A,A_trans,B,B_trans);
00597   MultiVectorMutable *Cmv = dynamic_cast<MultiVectorMutable*>(C);
00598   TEST_FOR_EXCEPTION(
00599     !Cmv || !(Cmv->access_by() & MultiVector::COL_ACCESS)
00600     ,MatrixOp::MethodNotImplemented
00601     ,"AbstractLinAlgPack::Mp_StMtM(...) : Error, mwo_lhs of type \'"
00602     << typeName(*C) << "\' does not support the "
00603     "\'MultiVectorMutable\' interface or does not support column access!" );
00604   // ToDo: We could do this by row also!
00605   VectorSpace::vec_mut_ptr_t
00606     t = ( B_trans == BLAS_Cpp::no_trans ? B.space_cols() : B.space_rows() ).create_member();
00607   const index_type
00608     C_rows = Cmv->rows(),
00609     C_cols = Cmv->cols();
00610   for( index_type j = 1; j <= C_cols; ++j ) {
00611     // t = op(B)*e(j)
00612     LinAlgOpPack::V_MtV( t.get(), B, B_trans, EtaVector(j,C_cols) );  
00613     // C(:,j) = a*op(A)*t + b*C(:,j)
00614     Vp_StMtV( Cmv->col(j).get(), a, A, A_trans, *t, b );
00615   }
00616 }
00617 
00618 void AbstractLinAlgPack::syrk(
00619   const MatrixOp  &A
00620   ,BLAS_Cpp::Transp   A_trans
00621   ,value_type         a
00622   ,value_type         b
00623   ,MatrixSymOp    *B
00624   )
00625 {
00626   // Give A a chance
00627   if(A.syrk(A_trans,a,b,B))
00628     return;
00629   // Give B a chance
00630   if(B->syrk(A,A_trans,a,b))
00631     return;
00632   
00633   TEST_FOR_EXCEPTION(
00634     true, MatrixOp::MethodNotImplemented
00635     ,"AbstractLinAlgPack::syrk(...) : Error, neither the right-hand-side matrix "
00636     "argument mwo_rhs of type \'" << typeName(A) << " nore the left-hand-side matrix "
00637     "argument sym_lhs of type \'" << typeName(*B) << "\' could implement this operation!"
00638     );
00639 
00640 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines