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