AbstractLinAlgPack_MatrixOpSerial.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 // ToDo: 7/27/99: Give these default implementations and test them.
00030 
00031 #include <typeinfo>
00032 
00033 #include "AbstractLinAlgPack_MatrixOpSerial.hpp"
00034 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00035 #include "AbstractLinAlgPack_MatrixOpGetGMSMutable.hpp"
00036 #include "AbstractLinAlgPack_MatrixOpGetGMSTri.hpp"
00037 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
00038 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00039 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00040 #include "AbstractLinAlgPack_EtaVector.hpp"
00041 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00042 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00043 #include "DenseLinAlgPack_DMatrixClass.hpp"
00044 #include "DenseLinAlgPack_DMatrixOut.hpp"
00045 #include "DenseLinAlgPack_AssertOp.hpp"
00046 #include "Teuchos_Workspace.hpp"
00047 #include "Teuchos_dyn_cast.hpp"
00048 
00049 namespace LinAlgOpPack {
00050   using AbstractLinAlgPack::Vp_StV;
00051   using AbstractLinAlgPack::Vp_StMtV;
00052   using AbstractLinAlgPack::Mp_StM;
00053 }
00054 
00055 namespace AbstractLinAlgPack {
00056 
00057 // Level-1 BLAS
00058 
00059 void MatrixOpSerial::Mp_StM(DMatrixSlice* gms_lhs, value_type alpha
00060   , BLAS_Cpp::Transp trans_rhs) const
00061 {
00062   DenseLinAlgPack::Mp_M_assert_sizes( gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00063     , rows(), cols(), trans_rhs );
00064   const size_type
00065     m = gms_lhs->rows(),
00066     n = gms_lhs->cols();
00067   //
00068   // Use sparse matrix-vector multiplication to perform this operation.
00069   // C += a * B = a * B * I = [ a*B*e(1), a*B*e(2), ..., a*B*e(m) ]
00070   //
00071   SpVector rhs;
00072   rhs.uninitialized_resize( n, 1, 1 );
00073   for( size_type j = 1; j <=n; ++j ) {
00074     rhs.begin()->initialize( j, 1.0 );  // e(j)
00075     this->Vp_StMtV( &gms_lhs->col(j), alpha, trans_rhs, rhs(), 1.0 );
00076   }
00077 }
00078 
00079 void MatrixOpSerial::Mp_StMtP(DMatrixSlice* C, value_type a
00080   , BLAS_Cpp::Transp M_trans
00081   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00082   ) const 
00083 {
00084   // C += a * op(M) * op(P)
00085   TEST_FOR_EXCEPT(true);  // Implement this!
00086 }
00087 
00088 void MatrixOpSerial::Mp_StPtM(DMatrixSlice* C, value_type a
00089   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00090   , BLAS_Cpp::Transp M_trans
00091   ) const 
00092 {
00093   // C += a * op(P) * op(M)
00094   TEST_FOR_EXCEPT(true);  // Implement this!
00095 }
00096 
00097 void MatrixOpSerial::Mp_StPtMtP( DMatrixSlice* C, value_type a
00098   , const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
00099   , BLAS_Cpp::Transp M_trans
00100   , const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
00101   ) const
00102 {
00103   // C += a * op(P1) * op(M) * op(P2)
00104   TEST_FOR_EXCEPT(true);  // Implement this!
00105 }
00106 
00107 // Level-2 BLAS
00108 
00109 void MatrixOpSerial::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha
00110   , BLAS_Cpp::Transp trans_rhs1, const SpVectorSlice& sv_rhs2, value_type beta) const
00111 {
00112   Vp_MtV_assert_sizes( vs_lhs->dim(), rows(), cols(), trans_rhs1, sv_rhs2.dim() );
00113   if( !sv_rhs2.nz() ) {
00114     // vs_lhs = beta * vs_lhs
00115     if(beta==0.0)      *vs_lhs = 0.0;
00116     else if(beta!=1.0) DenseLinAlgPack::Vt_S(vs_lhs,beta);
00117   }
00118   else {
00119     // Convert to dense by default.
00120     if( sv_rhs2.dim() == sv_rhs2.nz() && sv_rhs2.is_sorted() ) {
00121       const DVectorSlice vs_rhs2 = AbstractLinAlgPack::dense_view(sv_rhs2);
00122       this->Vp_StMtV( vs_lhs, alpha, trans_rhs1, vs_rhs2, beta );
00123     }
00124     else {
00125       DVector v_rhs2;
00126       LinAlgOpPack::assign( &v_rhs2, sv_rhs2 );
00127       this->Vp_StMtV( vs_lhs, alpha, trans_rhs1, v_rhs2(), beta );
00128     }
00129   }
00130 }
00131 
00132 void MatrixOpSerial::Vp_StPtMtV(DVectorSlice* y, value_type a
00133   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00134   , BLAS_Cpp::Transp M_trans
00135   , const DVectorSlice& x, value_type b) const
00136 {
00137   using Teuchos::Workspace;
00138   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00139 
00140   Workspace<value_type> t_ws(wss,BLAS_Cpp::cols(P.rows(),P.cols(),P_trans));
00141   DVectorSlice                t(&t_ws[0],t_ws.size());
00142     LinAlgOpPack::V_StMtV(&t,a,*this,M_trans,x);
00143   LinAlgOpPack::Vp_MtV( y, P, P_trans, t, b ); 
00144 }
00145 
00146 void MatrixOpSerial::Vp_StPtMtV(DVectorSlice* y, value_type a
00147   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00148   , BLAS_Cpp::Transp M_trans
00149   , const SpVectorSlice& x, value_type b) const
00150 {
00151   using Teuchos::Workspace;
00152   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00153 
00154   Workspace<value_type> t_ws(wss,BLAS_Cpp::cols(P.rows(),P.cols(),P_trans));
00155   DVectorSlice                t(&t_ws[0],t_ws.size());
00156     LinAlgOpPack::V_StMtV(&t,a,*this,M_trans,x);
00157   LinAlgOpPack::Vp_MtV( y, P, P_trans, t, b ); 
00158 }
00159 
00160 value_type MatrixOpSerial::transVtMtV(const DVectorSlice& x1
00161   , BLAS_Cpp::Transp M_trans, const DVectorSlice& x2) const
00162 {
00163   DenseLinAlgPack::Vp_MtV_assert_sizes( x1.dim(), rows(), cols(), M_trans, x2.dim() );
00164   DVector tmp(x1.dim());
00165   this->Vp_StMtV( &tmp(), 1.0, M_trans, x2, 0.0 );
00166   return DenseLinAlgPack::dot( x1, tmp() );
00167 }
00168 
00169 value_type MatrixOpSerial::transVtMtV(const SpVectorSlice& x1
00170   , BLAS_Cpp::Transp M_trans, const SpVectorSlice& x2) const
00171 {
00172   DenseLinAlgPack::Vp_MtV_assert_sizes( x1.dim(), rows(), cols(), M_trans, x2.dim() );
00173   if( !x1.nz() || !x2.nz() ) {
00174     return 0.0;
00175   }
00176   else {
00177     if( x1.overlap(x2) == DenseLinAlgPack::SAME_MEM && x1.dim() == x1.nz() && x1.is_sorted()  ) {
00178       const DVectorSlice x1_d = AbstractLinAlgPack::dense_view(x1);
00179       return this->transVtMtV( x1_d, M_trans, x1_d );
00180     }
00181     DVector tmp(x1.dim());
00182     this->Vp_StMtV( &tmp(), 1.0, M_trans, x2, 0.0 );
00183     return AbstractLinAlgPack::dot( x1, tmp() );
00184   }
00185 }
00186 
00187 void MatrixOpSerial::syr2k(
00188    BLAS_Cpp::Transp M_trans_in, value_type a
00189   , const GenPermMatrixSlice& P1_in, BLAS_Cpp::Transp P1_trans
00190   , const GenPermMatrixSlice& P2_in, BLAS_Cpp::Transp P2_trans
00191   , value_type b, DMatrixSliceSym* S ) const
00192 {
00193   using BLAS_Cpp::no_trans;
00194   using BLAS_Cpp::trans;
00195   using BLAS_Cpp::trans_not;
00196   using BLAS_Cpp::rows;
00197   using BLAS_Cpp::cols;
00198   //
00199   // S = b * S
00200   //
00201   // S += a*op(P1')*op(M)*op(P2) + a*op(P2')*op(M')*op(P1)
00202   //
00203   // We will start by renaming P1 and P2 such that op(P1).rows() >= op(P2).rows().
00204   // This is because we are going to store some temparary vectors and we don't
00205   // want them to be too big.
00206   //
00207   // We will perform the above operation by working with columns of:
00208   //
00209   //    op(P1)(:,j(k)) = e(i(k)) <: R^n
00210   //
00211   // Then for each column in op(P1) we will perform:
00212   //
00213   //
00214   // for k = 1...P1.nz()
00215   //
00216   //              [    .     ]
00217   //     S += a * [ e(i(k))' ] * op(M)*op(P2) + a * op(P2') * op(M') * [  ...  e(i(k))  ...  ] 
00218   //              [    .     ]
00219   //                row j(k)                                                    col j(k)
00220   //     =>
00221   //              [  .   ]
00222   //     S += a * [ y_k' ] + a * [  ...  y_k  ...  ] 
00223   //              [  .   ]
00224   //               row j(k)            col j(k)
00225   //
00226   //     where: y_k = a * op(P2') * op(M') * e(i(k)) <: R^m
00227   // 
00228   // Of course above we only need to set the row and column elements for S(j(k),:) and S(:,j(k))
00229   // for the portion of the symmetric S that is being stored.
00230   //
00231   const size_type
00232     M_rows  = this->rows(),
00233     M_cols  = this->cols(),
00234     P1_cols = cols( P1_in.rows(), P1_in.cols(), P1_trans );
00235   DenseLinAlgPack::MtM_assert_sizes(
00236     M_rows, M_cols, trans_not(M_trans_in)
00237     , P1_in.rows(), P1_in.cols(), P1_trans );
00238   DenseLinAlgPack::MtM_assert_sizes(
00239     M_rows, M_cols, M_trans_in
00240     , P2_in.rows(), P2_in.cols(), P2_trans );
00241   DenseLinAlgPack::Mp_M_assert_sizes(
00242     S->rows(), S->cols(), no_trans
00243     , P1_cols, P1_cols, no_trans );
00244   // Rename P1 and P2 so that op(P1).rows() >= op(P2).rows()
00245   const bool
00246     reorder_P1_P2 = ( rows( P1_in.rows(), P1_in.cols(), P1_trans ) 
00247               < rows( P2_in.rows(), P2_in.cols(), P2_trans ) );
00248   const GenPermMatrixSlice
00249     &P1 = reorder_P1_P2 ? P2_in : P1_in,
00250     &P2 = reorder_P1_P2 ? P1_in : P2_in;
00251   const BLAS_Cpp::Transp
00252     M_trans  = reorder_P1_P2 ? trans_not(M_trans_in) : M_trans_in;
00253   // Set rows and columns of S
00254   const size_type
00255     m  = S->rows(),
00256     n = rows( P1.rows(), P1.cols(), P1_trans );
00257   DVector y_k_store(m); // ToDo: use workspace allocator!
00258   DVectorSlice y_k = y_k_store();
00259   for( GenPermMatrixSlice::const_iterator P1_itr = P1.begin(); P1_itr != P1.end(); ++P1_itr )
00260   {
00261     const size_type
00262       i_k = P1_trans == no_trans ? P1_itr->row_i() : P1_itr->col_j(),
00263       j_k = P1_trans == no_trans ? P1_itr->col_j() : P1_itr->row_i();
00264     // e(i(k))
00265     EtaVector
00266       e_i_k(i_k,n);
00267     // y_k = op(P2')*op(M')*e(i(k))
00268     AbstractLinAlgPack::Vp_StPtMtV( &y_k, 1.0, P2, trans_not(P2_trans), *this, trans_not(M_trans), e_i_k(), 0.0 );
00269     // S(j(k),:) += a*y_k'
00270     if( S->uplo() == BLAS_Cpp::upper )
00271       DenseLinAlgPack::Vp_StV( &S->gms().row(j_k)(1,j_k), a, y_k(1,j_k) );
00272     else
00273       DenseLinAlgPack::Vp_StV( &S->gms().row(j_k)(j_k,m), a, y_k(j_k,m) );
00274     // S(:,j(k)) += a*y_k
00275     if( S->uplo() == BLAS_Cpp::upper )
00276       DenseLinAlgPack::Vp_StV( &S->gms().col(j_k)(1,j_k), a, y_k(1,j_k) );
00277     else
00278       DenseLinAlgPack::Vp_StV( &S->gms().col(j_k)(j_k,m), a, y_k(j_k,m) );
00279   }
00280 }
00281 
00282 // Level-3 BLAS
00283 
00284 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* C, value_type a
00285   , BLAS_Cpp::Transp A_trans, const DMatrixSlice& B
00286   , BLAS_Cpp::Transp B_trans, value_type b) const
00287 {
00288   DenseLinAlgPack::Mp_MtM_assert_sizes( C->rows(), C->cols(), BLAS_Cpp::no_trans
00289     , rows(), cols(), A_trans, B.rows(), B.cols(), B_trans );
00290   //
00291   // C = b*C + a*op(A)*op(B)
00292   // 
00293   // C.col(j) = b*col(j) + a*op(A)*op(B).col(j)
00294   //
00295 
00296   // ToDo: Add the ability to also perform by row if faster
00297 
00298   for( size_type j = 1; j <= C->cols(); ++j )
00299     AbstractLinAlgPack::Vp_StMtV( &C->col(j), a, *this, A_trans, DenseLinAlgPack::col( B, B_trans, j ), b );
00300 }
00301 
00302 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* C, value_type a, const DMatrixSlice& A
00303   , BLAS_Cpp::Transp A_trans, BLAS_Cpp::Transp B_trans, value_type b) const
00304 {
00305   DenseLinAlgPack::Mp_MtM_assert_sizes( C->rows(), C->cols(), BLAS_Cpp::no_trans
00306     , A.rows(), A.cols(), A_trans, rows(), cols(), B_trans );
00307   //
00308   // C = b*C + a*op(A)*op(B)
00309   // 
00310   // C.row(i) = b*row(i) + a*op(B)'*op(A).row(i)
00311   //
00312 
00313   // ToDo: Add the ability to also perform by column if faster
00314 
00315   for( size_type i = 1; i <= C->rows(); ++i )
00316     AbstractLinAlgPack::Vp_StMtV( &C->row(i), a, *this, BLAS_Cpp::trans_not(A_trans)
00317       , DenseLinAlgPack::row(A,A_trans,i) , b );
00318 }
00319 
00320 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* C, value_type a
00321   , BLAS_Cpp::Transp A_trans, const MatrixOpSerial& B
00322   , BLAS_Cpp::Transp B_trans, value_type b) const
00323 {
00324   using LinAlgOpPack::assign;
00325   // C = b*C + a*op(A)*op(B)
00326   DenseLinAlgPack::Mp_MtM_assert_sizes( C->rows(), C->cols(), BLAS_Cpp::no_trans
00327     , rows(), cols(), A_trans, B.rows(), B.cols(), B_trans );
00328   // Convert one of the matrices to dense, which ever one is the smallest!
00329   const size_type
00330     size_A = rows() * cols(),
00331     size_B = B.rows() * B.cols();
00332   if( size_A < size_B ) {
00333     DMatrix A_dense;
00334     assign( &A_dense, *this, BLAS_Cpp::no_trans );
00335     AbstractLinAlgPack::Mp_StMtM( C, a, A_dense(), A_trans, B, B_trans, b );
00336   }
00337   else {
00338     DMatrix B_dense;
00339     assign( &B_dense, B, BLAS_Cpp::no_trans );
00340     AbstractLinAlgPack::Mp_StMtM( C, a, *this, A_trans, B_dense(), B_trans, b );
00341   }
00342 }
00343 
00344 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha
00345   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceSym& sym_rhs2
00346   , BLAS_Cpp::Transp trans_rhs2, value_type beta) const
00347 {
00348   TEST_FOR_EXCEPT(true); // Todo: Implement!
00349 }
00350 
00351 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceSym& sym_rhs1
00352   , BLAS_Cpp::Transp trans_rhs1, BLAS_Cpp::Transp trans_rhs2, value_type beta) const
00353 {
00354   TEST_FOR_EXCEPT(true); // Todo: Implement!
00355 }
00356 
00357 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha
00358   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceTri& tri_rhs2
00359   , BLAS_Cpp::Transp trans_rhs2, value_type beta) const
00360 {
00361   TEST_FOR_EXCEPT(true); // Todo: Implement!
00362 }
00363 
00364 void MatrixOpSerial::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
00365   , BLAS_Cpp::Transp trans_rhs1, BLAS_Cpp::Transp trans_rhs2, value_type beta) const
00366 {
00367   TEST_FOR_EXCEPT(true); // Todo: Implement!
00368 }
00369 
00370 void MatrixOpSerial:: syrk(
00371   BLAS_Cpp::Transp M_trans, value_type a
00372   , value_type b, DMatrixSliceSym* S ) const
00373 {
00374   using BLAS_Cpp::no_trans;
00375   using BLAS_Cpp::trans;
00376   using BLAS_Cpp::trans_not;
00377   using BLAS_Cpp::rows;
00378   using BLAS_Cpp::cols;
00379   using Teuchos::Workspace;
00380   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00381   //
00382   // S = b*S + a*op(M)*op(M')
00383   //
00384   const size_type
00385     M_rows = this->rows(),
00386     M_cols = this->cols(),
00387     opM_rows = rows( M_rows, M_cols, M_trans ),
00388     opM_cols = cols( M_rows, M_cols, M_trans ),
00389     m = opM_rows;
00390   DenseLinAlgPack::Mp_MtM_assert_sizes(
00391     S->rows(), S->cols(), no_trans
00392     ,M_rows, M_cols, M_trans
00393     ,M_rows, M_cols, trans_not(M_trans) );
00394   // S = b*S
00395   DenseLinAlgPack::Mt_S( &DenseLinAlgPack::nonconst_tri_ele(S->gms(),S->uplo()), b );
00396   //
00397   // Form this matrix one column at a time by multiplying by e(j):
00398   //
00399   // S(:,j) += a*op(M)*(op(M')*e(j))
00400   //
00401   //    j = 1 ... opM_rows
00402   //
00403   Workspace<value_type> t1_ws(wss,opM_cols),
00404                            t2_ws(wss,opM_rows);
00405   DVectorSlice                t1(&t1_ws[0],t1_ws.size()),
00406                            t2(&t2_ws[0],t2_ws.size());
00407   for( size_type j = 1; j <= opM_rows; ++j ) {
00408     EtaVector e_j(j,opM_rows);
00409     LinAlgOpPack::V_MtV(&t1,*this,trans_not(M_trans),e_j()); // t1 = op(M')*e(j)
00410     LinAlgOpPack::V_MtV(&t2,*this,M_trans,t1);               // t2 = op(M)*t1
00411     // S(j,:) += a*t2' 
00412     if( S->uplo() == BLAS_Cpp::upper )
00413       DenseLinAlgPack::Vp_StV( &S->gms().row(j)(1,j), a, t2(1,j) );
00414     else
00415       DenseLinAlgPack::Vp_StV( &S->gms().row(j)(j,m), a, t2(j,m) );
00416     // S(:,j) += a*t2
00417     if( S->uplo() == BLAS_Cpp::upper )
00418       DenseLinAlgPack::Vp_StV( &S->gms().col(j)(1,j), a, t2(1,j) );
00419     else
00420       DenseLinAlgPack::Vp_StV( &S->gms().col(j)(j,m), a, t2(j,m) );
00421   }
00422 }
00423 
00424 // Overridden from MatrixOp
00425   
00426 const VectorSpace& MatrixOpSerial::space_cols() const
00427 {
00428   const size_type rows = this->rows();
00429   if(space_cols_.dim() != rows)
00430     space_cols_.initialize(rows);
00431   return space_cols_;
00432 }
00433   
00434 const VectorSpace& MatrixOpSerial::space_rows() const
00435 {
00436   const size_type cols = this->cols();
00437   if(space_rows_.dim() != cols)
00438     space_rows_.initialize(cols);
00439   return space_rows_;
00440 }
00441   
00442 std::ostream& MatrixOpSerial::output(std::ostream& out) const {
00443   DMatrix tmp( 0.0, rows(), cols() );
00444   this->Mp_StM( &tmp(), 1.0 , BLAS_Cpp::no_trans );
00445   return out << tmp();
00446 }
00447 
00448 bool MatrixOpSerial::Mp_StM(
00449   MatrixOp* mwo_lhs, value_type alpha
00450   ,BLAS_Cpp::Transp trans_rhs
00451   ) const
00452 {
00453   MatrixOpGetGMSMutable
00454     *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
00455   if(!mwo_gms_lhs)
00456     return MatrixOp::Mp_StM(mwo_lhs,alpha,trans_rhs); // boot it!
00457   this->Mp_StM( &MatrixDenseMutableEncap(mwo_gms_lhs)(), alpha, trans_rhs );
00458   return true;
00459 }
00460   
00461 bool MatrixOpSerial::Mp_StMtP(
00462   MatrixOp* mwo_lhs, value_type alpha
00463   , BLAS_Cpp::Transp M_trans
00464   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00465   ) const
00466 {
00467   MatrixOpGetGMSMutable
00468     *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
00469   if(!mwo_gms_lhs)
00470     return MatrixOp::Mp_StMtP(mwo_lhs,alpha,M_trans,P_rhs,P_rhs_trans); // boot it!
00471   this->Mp_StMtP(&MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,M_trans,P_rhs,P_rhs_trans);
00472   return true;
00473 }
00474   
00475 bool MatrixOpSerial::Mp_StPtM(
00476   MatrixOp* mwo_lhs, value_type alpha
00477   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00478   , BLAS_Cpp::Transp M_trans
00479   ) const
00480 {
00481   MatrixOpGetGMSMutable
00482     *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
00483   if(!mwo_gms_lhs)
00484     return MatrixOp::Mp_StPtM(mwo_lhs,alpha,P_rhs,P_rhs_trans,M_trans); // boot it!
00485   this->Mp_StPtM(&MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,P_rhs,P_rhs_trans,M_trans);
00486   return true;
00487 }
00488   
00489 bool MatrixOpSerial::Mp_StPtMtP(
00490   MatrixOp* mwo_lhs, value_type alpha
00491   ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00492   ,BLAS_Cpp::Transp M_trans
00493   ,const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
00494   ) const
00495 {
00496   MatrixOpGetGMSMutable
00497     *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
00498   if(!mwo_gms_lhs)
00499     return MatrixOp::Mp_StPtMtP(mwo_lhs,alpha,P_rhs1,P_rhs1_trans,M_trans,P_rhs2,P_rhs2_trans); // boot it!
00500   this->Mp_StPtMtP(&MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,P_rhs1,P_rhs1_trans,M_trans,P_rhs2,P_rhs2_trans);
00501   return true;
00502 }
00503   
00504 void MatrixOpSerial::Vp_StMtV(
00505   VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00506   , const Vector& v_rhs2, value_type beta) const
00507 {
00508   VectorDenseMutableEncap       vs_lhs(*v_lhs);
00509   VectorDenseEncap              vs_rhs2(v_rhs2);
00510   this->Vp_StMtV( &vs_lhs(), alpha, trans_rhs1, vs_rhs2(), beta );  
00511 }
00512 
00513 void MatrixOpSerial::Vp_StMtV(
00514   VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00515   , const SpVectorSlice& sv_rhs2, value_type beta) const
00516 {
00517   VectorDenseMutableEncap       vs_lhs(*v_lhs);
00518   this->Vp_StMtV( &vs_lhs(), alpha, trans_rhs1, sv_rhs2, beta );  
00519 }
00520   
00521 void MatrixOpSerial::Vp_StPtMtV(
00522   VectorMutable* v_lhs, value_type alpha
00523   , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00524   , BLAS_Cpp::Transp M_rhs2_trans
00525   , const Vector& v_rhs3, value_type beta) const
00526 {
00527   VectorDenseMutableEncap       vs_lhs(*v_lhs);
00528   VectorDenseEncap              vs_rhs3(v_rhs3);
00529   this->Vp_StPtMtV( &vs_lhs(), alpha, P_rhs1, P_rhs1_trans, M_rhs2_trans, vs_rhs3(), beta );  
00530 }
00531   
00532 void MatrixOpSerial::Vp_StPtMtV(
00533   VectorMutable* v_lhs, value_type alpha
00534   , const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00535   , BLAS_Cpp::Transp M_rhs2_trans
00536   , const SpVectorSlice& sv_rhs3, value_type beta) const
00537 {
00538   VectorDenseMutableEncap       vs_lhs(*v_lhs);
00539   this->Vp_StPtMtV( &vs_lhs(), alpha, P_rhs1, P_rhs1_trans, M_rhs2_trans, sv_rhs3, beta );  
00540 }
00541   
00542 value_type MatrixOpSerial::transVtMtV(
00543   const Vector& v_rhs1, BLAS_Cpp::Transp trans_rhs2
00544   , const Vector& v_rhs3) const
00545 {
00546   VectorDenseEncap              vs_rhs1(v_rhs1);
00547   VectorDenseEncap              vs_rhs3(v_rhs3);
00548   return this->transVtMtV(vs_rhs1(),trans_rhs2,vs_rhs3());
00549 }
00550   
00551 void MatrixOpSerial::syr2k(
00552   BLAS_Cpp::Transp M_trans, value_type alpha
00553   , const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
00554   , const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
00555   , value_type beta, MatrixSymOp* symwo_lhs ) const
00556 {
00557   MatrixSymOpGetGMSSymMutable
00558     *symwo_gms_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(symwo_lhs);
00559   if(!symwo_gms_lhs) {
00560     MatrixOp::syr2k(M_trans,alpha,P1,P1_trans,P2,P2_trans,beta,symwo_lhs); // Boot it
00561     return;
00562   }
00563   this->syr2k(
00564     M_trans,alpha,P1,P1_trans,P2,P2_trans,beta
00565     ,&MatrixDenseSymMutableEncap(symwo_gms_lhs)()
00566      );
00567 }
00568 
00569 bool MatrixOpSerial::Mp_StMtM(
00570   MatrixOp* mwo_lhs, value_type alpha
00571   , BLAS_Cpp::Transp trans_rhs1, const MatrixOp& mwo_rhs2
00572   , BLAS_Cpp::Transp trans_rhs2, value_type beta ) const
00573 {
00574   MatrixOpGetGMSMutable
00575     *mwo_gms_lhs = dynamic_cast<MatrixOpGetGMSMutable*>(mwo_lhs);
00576   if(mwo_gms_lhs) {
00577     // Try to match the rhs arguments to known serial interfaces
00578     if(const MatrixOpGetGMS* mwo_gms_rhs2 = dynamic_cast<const MatrixOpGetGMS*>(&mwo_rhs2)) {
00579       this->Mp_StMtM(
00580         &MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,trans_rhs1
00581         ,MatrixDenseEncap(*mwo_gms_rhs2)(),trans_rhs2,beta );
00582       return true;
00583     }
00584     if(const MatrixSymOpGetGMSSym* mwo_sym_gms_rhs2 = dynamic_cast<const MatrixSymOpGetGMSSym*>(&mwo_rhs2)) {
00585       this->Mp_StMtM(
00586         &MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,trans_rhs1
00587         ,MatrixDenseEncap(*mwo_sym_gms_rhs2)(),trans_rhs2,beta );
00588       return true;
00589     }
00590     if(const MatrixOpGetGMSTri* mwo_tri_gms_rhs2 = dynamic_cast<const MatrixOpGetGMSTri*>(&mwo_rhs2)) {
00591       this->Mp_StMtM(
00592         &MatrixDenseMutableEncap(mwo_gms_lhs)(),alpha,trans_rhs1
00593         ,MatrixDenseEncap(*mwo_tri_gms_rhs2)(),trans_rhs2,beta );
00594       return true;
00595     }
00596     // If we get here, the matrix arguments did not match up so we have to give up (I think?)
00597   }
00598   // Let the default implementation try to find matrix arguments that can handle this!
00599   return MatrixOp::Mp_StMtM(mwo_lhs,alpha,trans_rhs1,mwo_rhs2,trans_rhs2,beta); // Boot it!
00600 }
00601 
00602 bool MatrixOpSerial::syrk(
00603   BLAS_Cpp::Transp M_trans, value_type alpha
00604   , value_type beta, MatrixSymOp* symwo_lhs ) const
00605 {
00606   MatrixSymOpGetGMSSymMutable
00607     *symwo_gms_lhs = dynamic_cast<MatrixSymOpGetGMSSymMutable*>(symwo_lhs);
00608   if(!symwo_gms_lhs) {
00609     return MatrixOp::syrk(M_trans,alpha,beta,symwo_lhs); // Boot it
00610   }
00611   this->syrk(M_trans,alpha,beta,&MatrixDenseSymMutableEncap(symwo_gms_lhs)());
00612   return true;
00613 }
00614 
00615 } // end namespace AbstractLinAlgPack

Generated on Tue Jul 13 09:30:50 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7