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