MOOCHO (Single Doxygen Collection) Version of the Day
DenseLinAlgPack_DMatrixOpBLAS.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 "DenseLinAlgPack_DMatrixOp.hpp"
00043 #include "DenseLinAlgPack_BLAS_Cpp.hpp"
00044 
00045 namespace {
00046 
00047 using DenseLinAlgPack::DVector;
00048 using DenseLinAlgPack::DVectorSlice;
00049 using DenseLinAlgPack::DMatrix;
00050 using DenseLinAlgPack::DMatrixSlice;
00051 using DenseLinAlgPack::col;
00052 using DenseLinAlgPack::value_type;
00053 using DenseLinAlgPack::assert_gms_sizes;
00054 using DenseLinAlgPack::DMatrixSliceTriEle;
00055 using DenseLinAlgPack::DMatrixSliceTri;
00056 using DenseLinAlgPack::DMatrixSliceSym;
00057 using DenseLinAlgPack::assign;
00058 using BLAS_Cpp::Transp;
00059 using BLAS_Cpp::no_trans;
00060 using BLAS_Cpp::trans;
00061 using BLAS_Cpp::Uplo;
00062 using BLAS_Cpp::upper;
00063 using BLAS_Cpp::lower;
00064 
00065 } // end namespace
00066 
00067 // ///////////////////////////////////////////////////////////////////////////////////
00068 // Assignment Fucntions
00069 
00070 namespace {
00071 
00072 // implementation: gms_lhs = alpha (elementwise)
00073 inline void i_assign(DMatrixSlice* gms_lhs, value_type alpha) {
00074   for(DMatrixSlice::size_type i = 1; i <= gms_lhs->cols(); ++i)
00075     gms_lhs->col(i) = alpha;  
00076 }
00077 
00078 // implementaion: gms_lhs = gms_rhs. Most basic copy function for rectangular matrices.
00079 inline void i_assign_basic(DMatrixSlice* gms_lhs, const DMatrixSlice& gms_rhs
00080   ,  BLAS_Cpp::Transp trans_rhs)
00081 {
00082   for(DMatrixSlice::size_type i = 1; i <= gms_lhs->cols(); ++i)
00083       gms_lhs->col(i) = col(gms_rhs,trans_rhs,i);
00084 }
00085 
00086 // implementaion: gms_lhs = op(gms_rhs). Checks for overlap and creates temporaries accordingly.
00087 inline void i_assign(DMatrixSlice* gms_lhs, const DMatrixSlice& gms_rhs, BLAS_Cpp::Transp trans_rhs)
00088 {
00089   switch(gms_lhs->overlap(gms_rhs)) {
00090     case DenseLinAlgPack::NO_OVERLAP: // no overlap so just perform the copy
00091       i_assign_basic(gms_lhs,gms_rhs,trans_rhs);
00092       return;
00093     case DenseLinAlgPack::SAME_MEM:
00094       if(trans_rhs == BLAS_Cpp::no_trans) return; // assignment to self, nothing to do.
00095     default: // either same memory that needs to be transposed or some overlap so just generate temp.
00096       DMatrix temp = gms_rhs;
00097       i_assign_basic(gms_lhs,temp(),trans_rhs);
00098       return;
00099   }
00100 }
00101 
00102 // Copy one triangular region into another. Does not check sizes or aliasing of argument matrices.
00103 // A row of a upper triangular region corresponds to a col of a BLAS_Cpp::lower triangular region.
00104 inline void i_assign_basic(DMatrixSliceTriEle* tri_lhs, const DMatrixSliceTriEle& tri_rhs)
00105 {
00106   DMatrixSlice::size_type n = tri_lhs->gms().cols();
00107 
00108   // Access BLAS_Cpp::lower tri by col and upper tri by row
00109   BLAS_Cpp::Transp
00110     trans_lhs = (tri_lhs->uplo() == BLAS_Cpp::lower) ? BLAS_Cpp::no_trans : BLAS_Cpp::trans,
00111     trans_rhs = (tri_rhs.uplo() == BLAS_Cpp::lower) ? BLAS_Cpp::no_trans : BLAS_Cpp::trans;
00112   
00113   for(int i = 1; i <= n; ++i) { // Only copy the part of the vec in tri region.
00114     col(tri_lhs->gms(),trans_lhs,i)(i,n) = col(tri_rhs.gms(),trans_rhs,i)(i,n);
00115   }
00116 }
00117 
00118 } // end namespace
00119 
00120 // gm_lhs = alpha (elementwise)
00121 void DenseLinAlgPack::assign(DMatrix* gm_lhs, value_type alpha)
00122 {
00123   if(!gm_lhs->rows()) gm_lhs->resize(1,1,alpha);
00124   i_assign(&(*gm_lhs)(),alpha);
00125 }
00126 
00127 // gm_lhs = op(gms_rhs)
00128 void DenseLinAlgPack::assign(DMatrix* gm_lhs, const DMatrixSlice& gms_rhs, BLAS_Cpp::Transp trans_rhs)
00129 {
00130   if(gm_lhs->overlap(gms_rhs) == SAME_MEM && trans_rhs == BLAS_Cpp::no_trans) return; // assignment to self
00131   if(gm_lhs->overlap(gms_rhs) != NO_OVERLAP) {
00132     // some overlap so we must create a copy
00133     DMatrix tmp(gms_rhs);
00134     resize_gm_lhs(gm_lhs,gms_rhs.rows(),gms_rhs.cols(),trans_rhs);
00135     i_assign(&(*gm_lhs)(), tmp(), trans_rhs);
00136   }
00137   else {
00138     // no overlap so just assign
00139     resize_gm_lhs(gm_lhs,gms_rhs.rows(),gms_rhs.cols(),trans_rhs);
00140     i_assign(&(*gm_lhs)(), gms_rhs, trans_rhs);
00141   }
00142 }
00143 
00144 // gms_lhs = alpha (elementwise)
00145 void DenseLinAlgPack::assign(DMatrixSlice* gms_lhs, value_type alpha)
00146 {
00147   TEUCHOS_TEST_FOR_EXCEPTION(
00148     !gms_lhs->rows(), std::length_error
00149     ,"Error, assign(gms_lhs,alpha):  You can not assign Scalar to an unsized DMatrixSlice" );
00150   i_assign(gms_lhs, alpha);
00151 }
00152 
00153 // gms_lhs = op(gms_rhs)
00154 void DenseLinAlgPack::assign(DMatrixSlice* gms_lhs, const DMatrixSlice& gms_rhs, BLAS_Cpp::Transp trans_rhs)
00155 {
00156   assert_gms_lhs(*gms_lhs,gms_rhs.rows(),gms_rhs.cols(),trans_rhs);
00157   i_assign(gms_lhs, gms_rhs, trans_rhs);
00158 }
00159 
00160 // tri_lhs = alpha (elementwise)
00161 void DenseLinAlgPack::assign(DMatrixSliceTriEle* tri_lhs, value_type alpha)
00162 {
00163   TEUCHOS_TEST_FOR_EXCEPTION(
00164     !tri_lhs->gms().rows(), std::length_error
00165     ,"Error, assign(tri_lhs,alpha):  You can not assign a scalar to an unsized triangular DMatrixSlice" );
00166   assert_gms_square(tri_lhs->gms()); 
00167   DMatrixSlice::size_type n = tri_lhs->gms().cols();
00168   // access BLAS_Cpp::lower tri by col and upper tri by row
00169   BLAS_Cpp::Transp
00170     trans_lhs = (tri_lhs->uplo() == BLAS_Cpp::lower) ? BLAS_Cpp::no_trans : BLAS_Cpp::trans;
00171   for(int i = 1; i <= n; ++i)
00172     col( tri_lhs->gms(), trans_lhs , i )(i,n) = alpha;
00173 }
00174 
00175 // tri_lhs = tri_rhs
00176 void DenseLinAlgPack::assign(DMatrixSliceTriEle* tri_lhs, const DMatrixSliceTriEle& tri_rhs)
00177 {
00178   assert_gms_lhs(tri_lhs->gms(),tri_rhs.gms().rows(),tri_rhs.gms().cols(),BLAS_Cpp::no_trans);
00179     
00180   switch(tri_lhs->gms().overlap(tri_rhs.gms())) {
00181     case SAME_MEM:
00182       if(tri_lhs->uplo() == tri_rhs.uplo()) return; // Assignment to self to exit
00183 
00184     case SOME_OVERLAP:
00185       // Test for the special case where the upper tri region (above diagonal) of a
00186       // matrix is being copied into the BLAS_Cpp::lower tri region (below the diagonal) of
00187       // the same matrix or visa-versa.  No temp is need in this case
00188       if(tri_lhs->uplo() != tri_rhs.uplo()) {
00189         const DMatrixSlice  *up = (tri_lhs->uplo() == upper)
00190                         ? &tri_lhs->gms() : &tri_rhs.gms();
00191         const DMatrixSlice  *lo = (tri_rhs.uplo() == BLAS_Cpp::lower)
00192                         ? &tri_lhs->gms() : &tri_rhs.gms();
00193         if(lo->col_ptr(1) + lo->max_rows() - 1 == up->col_ptr(1)) { // false if SAME_MEM
00194           // One triangular region is being copied into another so no temp needed.
00195           i_assign_basic(tri_lhs, tri_rhs);
00196           return; 
00197         }
00198       }
00199       // Give up and copy the vs_rhs as a temp.
00200       {
00201         DMatrix temp(tri_rhs.gms());
00202         i_assign_basic(tri_lhs, tri_ele(temp(),tri_rhs.uplo()));
00203         return;
00204       }
00205 
00206     case NO_OVERLAP:  // no overlap so just perform the copy
00207       i_assign_basic(tri_lhs, tri_rhs);
00208       return;
00209   }
00210 }
00211 
00212 // /////////////////////////////////////////////////////////////////////////////////////////
00213 // Element-wise Algebraic Operations
00214 
00215 void DenseLinAlgPack::Mt_S(DMatrixSlice* gms_lhs, value_type alpha)
00216 {
00217   TEUCHOS_TEST_FOR_EXCEPTION(
00218     !gms_lhs->rows(), std::length_error
00219     ,"Error, Mt_S(gms_lhs,alpha):  You can not scale an unsized DMatrixSlice" );
00220   for(int j = 1; j <= gms_lhs->cols(); ++j)
00221     Vt_S(&gms_lhs->col(j), alpha);
00222 }
00223 
00224 void DenseLinAlgPack::M_diagVtM( DMatrixSlice* gms_lhs, const DVectorSlice& vs_rhs
00225                 , const DMatrixSlice& gms_rhs, BLAS_Cpp::Transp trans_rhs )
00226 {
00227   Mp_M_assert_sizes(gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00228     , gms_rhs.rows(), gms_rhs.cols(), trans_rhs);
00229   for(DMatrixSlice::size_type j = 1; j <= gms_lhs->cols(); ++j)
00230     prod( &gms_lhs->col(j), vs_rhs, col(gms_rhs,trans_rhs,j) );
00231 }
00232 
00233 void DenseLinAlgPack::Mt_S(DMatrixSliceTriEle* tri_lhs, value_type alpha) {
00234   TEUCHOS_TEST_FOR_EXCEPTION(
00235     !tri_lhs->gms().rows(), std::length_error
00236     ,"Error, Mt_S(tri_lhs,alpha):  You can not scale an unsized triangular DMatrixSlice" );
00237   BLAS_Cpp::Transp
00238     trans_lhs = (tri_lhs->uplo() == BLAS_Cpp::lower) ? BLAS_Cpp::no_trans : BLAS_Cpp::trans;
00239 
00240   DMatrixSlice::size_type n = tri_lhs->gms().cols();
00241   for(DMatrixSlice::size_type j = 1; j <= n; ++j)
00242     Vt_S( &col( tri_lhs->gms(), trans_lhs , j )(j,n), alpha );  
00243 }
00244 
00245 void DenseLinAlgPack::Mp_StM(DMatrixSliceTriEle* tri_lhs, value_type alpha, const DMatrixSliceTriEle& tri_ele_rhs)
00246 {
00247   Mp_M_assert_sizes(tri_lhs->gms().rows(), tri_lhs->gms().cols(), BLAS_Cpp::no_trans
00248     , tri_ele_rhs.gms().rows(), tri_ele_rhs.gms().cols(), BLAS_Cpp::no_trans);
00249 
00250   BLAS_Cpp::Transp
00251     trans_lhs = (tri_lhs->uplo() == BLAS_Cpp::lower) ? BLAS_Cpp::no_trans : BLAS_Cpp::trans,
00252     trans_rhs = (tri_ele_rhs.uplo() == BLAS_Cpp::lower) ? BLAS_Cpp::no_trans : BLAS_Cpp::trans;
00253 
00254   DMatrixSlice::size_type n = tri_lhs->gms().cols();
00255   for(DMatrixSlice::size_type j = 1; j <= n; ++j)
00256     Vp_StV( &col(tri_lhs->gms(),trans_lhs,j)(j,n), alpha, col(tri_ele_rhs.gms(),trans_rhs,j)(j,n) );  
00257 }
00258 
00259 // LinAlgOpPack Compatable (compile time polymorphism)
00260 
00261 void DenseLinAlgPack::Mp_StM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSlice& gms_rhs
00262   , BLAS_Cpp::Transp trans_rhs)
00263 {
00264   Mp_M_assert_sizes(gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00265     , gms_rhs.rows(), gms_rhs.cols(), trans_rhs);
00266   for(DMatrixSlice::size_type j = 1; j <= gms_lhs->cols(); ++j)
00267     Vp_StV( &gms_lhs->col(j), alpha, col(gms_rhs,trans_rhs,j) );  
00268 }
00269 
00270 namespace {
00271 // Implement upper and lower copies for a symmetric matrix
00272 // Does not check sizes.
00273 // inline
00274 void i_Mp_StSM( DenseLinAlgPack::DMatrixSlice* gms_lhs, DenseLinAlgPack::value_type alpha
00275   , const DenseLinAlgPack::DMatrixSliceTriEle& tri_ele_rhs)
00276 {
00277   using DenseLinAlgPack::Mp_StM;
00278   using DenseLinAlgPack::tri_ele;
00279   const DenseLinAlgPack::size_type n = gms_lhs->rows(); // same as cols
00280   Mp_StM( const_cast<DenseLinAlgPack::DMatrixSliceTriEle*>( 
00281         &tri_ele((*gms_lhs)(2,n,1,n-1), BLAS_Cpp::lower ) )
00282       , alpha, tri_ele_rhs );
00283   Mp_StM( const_cast<DenseLinAlgPack::DMatrixSliceTriEle*>(
00284         &tri_ele((*gms_lhs)(1,n-1,2,n), BLAS_Cpp::upper ) )
00285       , alpha, tri_ele_rhs );
00286 }
00287 } // end namespace
00288 
00289 void DenseLinAlgPack::Mp_StM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceSym& sym_rhs
00290   , BLAS_Cpp::Transp trans_rhs )
00291 {
00292   Mp_M_assert_sizes(gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00293     , sym_rhs.rows(), sym_rhs.cols(), trans_rhs);
00294   Vp_StV( &gms_lhs->diag(), alpha, sym_rhs.gms().diag() );
00295   const size_type n = gms_lhs->rows(); // same as cols
00296   switch(sym_rhs.uplo()) {
00297     case BLAS_Cpp::lower:
00298       i_Mp_StSM( gms_lhs, alpha, tri_ele(sym_rhs.gms()(2,n,1,n-1),BLAS_Cpp::lower) );
00299       break;
00300     case BLAS_Cpp::upper:
00301       i_Mp_StSM( gms_lhs, alpha, tri_ele(sym_rhs.gms()(1,n-1,2,n),BLAS_Cpp::upper) );
00302       break;
00303   }
00304 }
00305 
00306 void DenseLinAlgPack::Mp_StM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs
00307   , BLAS_Cpp::Transp trans_rhs)
00308 {
00309   using BLAS_Cpp::trans;
00310   using BLAS_Cpp::no_trans;
00311   using BLAS_Cpp::lower;
00312   using BLAS_Cpp::upper;
00313   Mp_M_assert_sizes(gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00314     , tri_rhs.rows(), tri_rhs.cols(), trans_rhs );
00315   // diagonal
00316   if( tri_rhs.diag() == BLAS_Cpp::nonunit )
00317     Vp_StV( &gms_lhs->diag(), alpha, tri_rhs.gms().diag() );
00318   else
00319     Vp_S( &gms_lhs->diag(), alpha );
00320   // upper or lower triangle
00321   const size_type n = gms_lhs->rows(); // same as cols
00322   if( n == 1 )
00323     return; // There is not lower or upper triangular region
00324   const DMatrixSliceTriEle
00325     tri = ( tri_rhs.uplo() == lower
00326         ?   tri_ele(tri_rhs.gms()(2,n,1,n-1),lower)
00327            :  tri_ele(tri_rhs.gms()(1,n-1,2,n),upper) );
00328   const BLAS_Cpp::Uplo
00329     as_uplo = (   ( tri_rhs.uplo() == lower && trans_rhs == no_trans )
00330           ||  ( tri_rhs.uplo() == upper && trans_rhs == trans )
00331           ? lower : upper                     );
00332   switch(as_uplo) {
00333     case lower:
00334       Mp_StM( const_cast<DenseLinAlgPack::DMatrixSliceTriEle*>(
00335             &tri_ele((*gms_lhs)(2,n,1,n-1),lower) )
00336           , alpha, tri );
00337       break;
00338     case upper:
00339       Mp_StM( const_cast<DenseLinAlgPack::DMatrixSliceTriEle*>(
00340             &tri_ele((*gms_lhs)(1,n-1,2,n),upper) )
00341           , alpha, tri );
00342       break;
00343   }
00344 }
00345 
00346 // /////////////////////////////////////////////////////////////////////////////////////
00347 // Level-2 BLAS (vector-matrtix) Liner Algebra Operations
00348 
00349 void DenseLinAlgPack::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00350   , BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2, value_type beta)
00351 {
00352   Vp_MtV_assert_sizes(vs_lhs->dim(), gms_rhs1.rows()  , gms_rhs1.cols(), trans_rhs1
00353     , vs_rhs2.dim());
00354   BLAS_Cpp::gemv(trans_rhs1,gms_rhs1.rows(),gms_rhs1.cols(),alpha,gms_rhs1.col_ptr(1)
00355     ,gms_rhs1.max_rows(), vs_rhs2.raw_ptr(),vs_rhs2.stride(),beta,vs_lhs->raw_ptr()
00356     ,vs_lhs->stride());
00357 }
00358 
00359 void DenseLinAlgPack::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, const DMatrixSliceSym& sym_rhs1
00360   , BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2, value_type beta)
00361 {
00362   assert_gms_square(sym_rhs1.gms());
00363   Vp_MtV_assert_sizes(vs_lhs->dim(), sym_rhs1.gms().rows(), sym_rhs1.gms().cols(), trans_rhs1
00364     , vs_rhs2.dim());
00365   BLAS_Cpp::symv(sym_rhs1.uplo(),sym_rhs1.gms().rows(),alpha,sym_rhs1.gms().col_ptr(1)
00366     ,sym_rhs1.gms().max_rows(),vs_rhs2.raw_ptr(),vs_rhs2.stride(),beta
00367     ,vs_lhs->raw_ptr(),vs_lhs->stride());
00368 }
00369 
00370 void DenseLinAlgPack::V_MtV(DVector* v_lhs, const DMatrixSliceTri& tri_rhs1, BLAS_Cpp::Transp trans_rhs1
00371   , const DVectorSlice& vs_rhs2)
00372 {
00373   assert_gms_square(tri_rhs1.gms());
00374   MtV_assert_sizes(tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1, vs_rhs2.dim());
00375   v_lhs->resize(tri_rhs1.gms().rows());
00376   (*v_lhs) = vs_rhs2;
00377   BLAS_Cpp::trmv(tri_rhs1.uplo(),trans_rhs1,tri_rhs1.diag(),tri_rhs1.gms().rows()
00378     ,tri_rhs1.gms().col_ptr(1),tri_rhs1.gms().max_rows(), v_lhs->raw_ptr(),1);
00379 }
00380 
00381 void DenseLinAlgPack::V_MtV(DVectorSlice* vs_lhs, const DMatrixSliceTri& tri_rhs1, BLAS_Cpp::Transp trans_rhs1
00382   , const DVectorSlice& vs_rhs2)
00383 {
00384   assert_gms_square(tri_rhs1.gms());
00385   MtV_assert_sizes(tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1, vs_rhs2.dim());
00386   Vp_V_assert_sizes( vs_lhs->dim(), tri_rhs1.gms().rows() );
00387   (*vs_lhs) = vs_rhs2;
00388   BLAS_Cpp::trmv(tri_rhs1.uplo(),trans_rhs1,tri_rhs1.diag(),tri_rhs1.gms().rows()
00389     ,tri_rhs1.gms().col_ptr(1),tri_rhs1.gms().max_rows(), vs_lhs->raw_ptr(),vs_lhs->stride());
00390 }
00391 
00392 void DenseLinAlgPack::Vp_StMtV(DVectorSlice* vs_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
00393   , BLAS_Cpp::Transp trans_rhs1, const DVectorSlice& vs_rhs2, value_type beta)
00394 {
00395   Vp_MtV_assert_sizes(vs_lhs->dim(),tri_rhs1.gms().rows(),tri_rhs1.gms().cols()
00396             ,trans_rhs1,vs_rhs2.dim() );
00397 
00398   // If op(gms_rhs2) == gms_lhs and beta = 0.0 then this is a direct call to the BLAS.
00399   if( vs_lhs->overlap(vs_rhs2) == SAME_MEM && beta == 0.0 )
00400   {
00401     V_MtV(vs_lhs, tri_rhs1, trans_rhs1, vs_rhs2);
00402     if(alpha != 1.0) Vt_S(vs_lhs,alpha);
00403   }
00404   else {
00405     // This is something else so the alias problem is not handled.
00406     if(beta != 1.0) Vt_S(vs_lhs,beta);
00407     DVector tmp;
00408     V_MtV(&tmp, tri_rhs1, trans_rhs1, vs_rhs2);
00409     Vp_StV(vs_lhs,alpha,tmp());
00410   }
00411 }
00412 
00413 void DenseLinAlgPack::V_InvMtV(DVector* v_lhs, const DMatrixSliceTri& tri_rhs1, BLAS_Cpp::Transp trans_rhs1
00414   , const DVectorSlice& vs_rhs2)
00415 {
00416   assert_gms_square(tri_rhs1.gms());
00417   MtV_assert_sizes(tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1, vs_rhs2.dim());
00418   v_lhs->resize(tri_rhs1.gms().rows());
00419   (*v_lhs) = vs_rhs2;
00420   BLAS_Cpp::trsv(tri_rhs1.uplo(),trans_rhs1,tri_rhs1.diag(),tri_rhs1.gms().rows()
00421     ,tri_rhs1.gms().col_ptr(1),tri_rhs1.gms().max_rows(), v_lhs->raw_ptr(),1);
00422 }
00423 
00424 void DenseLinAlgPack::V_InvMtV(DVectorSlice* vs_lhs, const DMatrixSliceTri& tri_rhs1, BLAS_Cpp::Transp trans_rhs1
00425   , const DVectorSlice& vs_rhs2)
00426 {
00427   assert_gms_square(tri_rhs1.gms());
00428   MtV_assert_sizes(tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1, vs_rhs2.dim());
00429   Vp_V_assert_sizes( vs_lhs->dim(), tri_rhs1.gms().rows() );
00430   (*vs_lhs) = vs_rhs2;
00431   BLAS_Cpp::trsv(tri_rhs1.uplo(),trans_rhs1,tri_rhs1.diag(),tri_rhs1.gms().rows()
00432     ,tri_rhs1.gms().col_ptr(1),tri_rhs1.gms().max_rows(), vs_lhs->raw_ptr(),vs_lhs->stride());
00433 }
00434 
00435 
00436 void DenseLinAlgPack::ger(
00437   value_type alpha, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2
00438   , DMatrixSlice* gms_lhs )
00439 {
00440   Vp_MtV_assert_sizes( vs_rhs2.dim(),  gms_lhs->rows(), gms_lhs->cols()
00441     , BLAS_Cpp::no_trans, vs_rhs1.dim() );
00442   BLAS_Cpp::ger(
00443     gms_lhs->rows(), gms_lhs->cols(), alpha
00444     ,vs_rhs1.raw_ptr(), vs_rhs1.stride()
00445     ,vs_rhs2.raw_ptr(), vs_rhs2.stride()
00446     ,gms_lhs->col_ptr(1), gms_lhs->max_rows() );
00447 }
00448 
00449 void DenseLinAlgPack::syr(value_type alpha, const DVectorSlice& vs_rhs, DMatrixSliceSym* sym_lhs)
00450 {
00451   assert_gms_square(sym_lhs->gms());
00452   MtV_assert_sizes( sym_lhs->gms().rows(), sym_lhs->gms().cols()
00453     , BLAS_Cpp::no_trans, vs_rhs.dim() );
00454   BLAS_Cpp::syr( sym_lhs->uplo(), vs_rhs.dim(), alpha, vs_rhs.raw_ptr()
00455     , vs_rhs.stride(), sym_lhs->gms().col_ptr(1), sym_lhs->gms().max_rows() );
00456 }
00457 
00458 void DenseLinAlgPack::syr2(value_type alpha, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2
00459   , DMatrixSliceSym* sym_lhs)
00460 {
00461   assert_gms_square(sym_lhs->gms());
00462   VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
00463   MtV_assert_sizes( sym_lhs->gms().rows(), sym_lhs->gms().cols()
00464     , BLAS_Cpp::no_trans, vs_rhs1.dim() );
00465   BLAS_Cpp::syr2( sym_lhs->uplo(), vs_rhs1.dim(), alpha, vs_rhs1.raw_ptr()
00466     , vs_rhs1.stride(), vs_rhs2.raw_ptr(), vs_rhs2.stride()
00467     , sym_lhs->gms().col_ptr(1), sym_lhs->gms().max_rows() );
00468 }
00469 
00470 // //////////////////////////////////////////////////////////////////////////////////////////
00471 // Level-3 BLAS (matrix-matrix) Linear Algebra Operations
00472 
00473 // ////////////////////////////
00474 // Rectangular Matrices
00475 
00476 void DenseLinAlgPack::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00477   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
00478   , BLAS_Cpp::Transp trans_rhs2, value_type beta)
00479 {
00480   Mp_MtM_assert_sizes(    gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00481               , gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1
00482               , gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2);
00483   BLAS_Cpp::gemm(trans_rhs1,trans_rhs2,gms_lhs->rows(),gms_lhs->cols()
00484     ,cols(gms_rhs1.rows(),gms_rhs1.cols(),trans_rhs1)
00485     ,alpha,gms_rhs1.col_ptr(1),gms_rhs1.max_rows()
00486     ,gms_rhs2.col_ptr(1),gms_rhs2.max_rows()
00487     ,beta,gms_lhs->col_ptr(1),gms_lhs->max_rows() );
00488 }
00489 
00490 // ////////////////////////////
00491 // Symmetric Matrices
00492 
00493 namespace {
00494 
00495 // implementation:  gms_lhs = alpha * sym_rhs * gms_rhs + beta * gms_lhs (left) (BLAS xSYMM).
00496 // or       gms_lhs = alpha * gms_rhs * sym_rhs + beta * gms_lhs (right).
00497 // does not check sizes.
00498 void i_symm(BLAS_Cpp::Side side, value_type alpha, const DMatrixSliceSym& sym_rhs
00499   , const DMatrixSlice& gms_rhs, value_type beta, DMatrixSlice* gms_lhs)
00500 {
00501   BLAS_Cpp::symm(side,sym_rhs.uplo(),gms_lhs->rows(),gms_lhs->cols()
00502     ,alpha,sym_rhs.gms().col_ptr(1),sym_rhs.gms().max_rows()
00503     ,gms_rhs.col_ptr(1),gms_rhs.max_rows()
00504     ,beta,gms_lhs->col_ptr(1),gms_lhs->max_rows() );
00505 }
00506 
00507 } // end namespace
00508 
00509 void DenseLinAlgPack::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceSym& sym_rhs1
00510   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
00511   , BLAS_Cpp::Transp trans_rhs2, value_type beta)
00512 {
00513   Mp_MtM_assert_sizes(    gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00514               , sym_rhs1.gms().rows(), sym_rhs1.gms().cols(), trans_rhs1
00515               , gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2);
00516   if(trans_rhs2 == BLAS_Cpp::no_trans) {
00517     i_symm(BLAS_Cpp::left,alpha,sym_rhs1,gms_rhs2,beta,gms_lhs);
00518   }
00519   else {
00520     // must make a temporary copy to call the BLAS
00521     DMatrix tmp;
00522     assign(&tmp,gms_rhs2,trans_rhs2);
00523     i_symm(BLAS_Cpp::left,alpha,sym_rhs1,tmp(),beta,gms_lhs);
00524   }
00525 }
00526 
00527 void DenseLinAlgPack::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00528   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceSym& sym_rhs2
00529   , BLAS_Cpp::Transp trans_rhs2, value_type beta)
00530 {
00531   Mp_MtM_assert_sizes(    gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00532               , gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1
00533               , sym_rhs2.gms().rows(), sym_rhs2.gms().cols(), trans_rhs2 );
00534   if(trans_rhs1 == BLAS_Cpp::no_trans) {
00535     i_symm(BLAS_Cpp::right,alpha,sym_rhs2,gms_rhs1,beta,gms_lhs);
00536   }
00537   else {
00538     // must make a temporary copy to call the BLAS
00539     DMatrix tmp;
00540     assign(&tmp,gms_rhs1,trans_rhs1);
00541     i_symm(BLAS_Cpp::right,alpha,sym_rhs2,tmp(),beta,gms_lhs);
00542   }
00543 }
00544 
00545 void DenseLinAlgPack::syrk(BLAS_Cpp::Transp trans, value_type alpha, const DMatrixSlice& gms_rhs
00546   , value_type beta, DMatrixSliceSym* sym_lhs)
00547 {
00548   Mp_MtM_assert_sizes(    sym_lhs->gms().rows(), sym_lhs->gms().cols(), BLAS_Cpp::no_trans
00549               , gms_rhs.rows(), gms_rhs.cols(), trans
00550               , gms_rhs.rows(), gms_rhs.cols(), trans_not(trans) );
00551   BLAS_Cpp::syrk(sym_lhs->uplo(),trans,sym_lhs->gms().rows()
00552     ,cols(gms_rhs.rows(), gms_rhs.cols(), trans),alpha
00553     ,gms_rhs.col_ptr(1),gms_rhs.max_rows(),beta
00554     ,sym_lhs->gms().col_ptr(1),sym_lhs->gms().max_rows() );
00555 }
00556 
00557 void DenseLinAlgPack::syr2k(BLAS_Cpp::Transp trans,value_type alpha, const DMatrixSlice& gms_rhs1
00558   , const DMatrixSlice& gms_rhs2, value_type beta, DMatrixSliceSym* sym_lhs)
00559 {
00560   Mp_MtM_assert_sizes(    sym_lhs->gms().rows(), sym_lhs->gms().cols(), BLAS_Cpp::no_trans
00561               , gms_rhs1.rows(), gms_rhs1.cols(), trans
00562               , gms_rhs2.rows(), gms_rhs2.cols(), trans_not(trans) );
00563   BLAS_Cpp::syr2k(sym_lhs->uplo(),trans,sym_lhs->gms().rows()
00564     ,cols(gms_rhs1.rows(), gms_rhs1.cols(), trans),alpha
00565     ,gms_rhs1.col_ptr(1),gms_rhs1.max_rows()
00566     ,gms_rhs2.col_ptr(1),gms_rhs2.max_rows(),beta
00567     ,sym_lhs->gms().col_ptr(1),sym_lhs->gms().max_rows() );
00568 }
00569 
00570 // ////////////////////////////
00571 // Triangular Matrices
00572 
00573 // ToDo: Finish the definitions below.
00574 
00575 namespace {
00576 
00577 // implementation:  gms_lhs = alpha * op(tri_rhs) * gms_lhs (left) (BLAS xTRMM).
00578 // or       gms_lhs = alpha * gms_rhs * op(tri_rhs) (right).
00579 // does not check sizes.
00580 inline
00581 void i_trmm(BLAS_Cpp::Side side, BLAS_Cpp::Transp trans, value_type alpha, const DMatrixSliceTri& tri_rhs
00582   , DMatrixSlice* gms_lhs)
00583 {
00584   BLAS_Cpp::trmm(side,tri_rhs.uplo(),trans,tri_rhs.diag(),gms_lhs->rows(),gms_lhs->cols()
00585     ,alpha,tri_rhs.gms().col_ptr(1),tri_rhs.gms().max_rows()
00586     ,gms_lhs->col_ptr(1),gms_lhs->max_rows() );
00587 }
00588 
00589 // implementation:  gms_lhs = alpha * op(tri_rhs) * op(gms_rhs) (left)
00590 // or       gms_lhs = alpha * op(gms_rhs) * op(tri_rhs) (right)
00591 // Takes care of temporaries but does not check sizes.
00592 inline
00593 void i_trmm_alt( BLAS_Cpp::Side side, value_type alpha, const DMatrixSliceTri& tri_rhs
00594   , BLAS_Cpp::Transp trans_tri_rhs, const DMatrixSlice& gms_rhs
00595   , BLAS_Cpp::Transp trans_gms_rhs, DMatrixSlice* gms_lhs )
00596 {
00597   assign(gms_lhs,gms_rhs,trans_gms_rhs);
00598   i_trmm( side, trans_tri_rhs, alpha, tri_rhs, gms_lhs );
00599 }
00600 
00601 // implementation:  gms_lhs = alpha * inv(op(tri_rhs)) * gms_lhs (left) (BLAS xTRSM).
00602 // or       gms_lhs = alpha * gms_rhs * inv(op(tri_rhs)) (right).
00603 // does not check sizes.
00604 inline
00605 void i_trsm(BLAS_Cpp::Side side, BLAS_Cpp::Transp trans, value_type alpha, const DMatrixSliceTri& tri_rhs
00606   , DMatrixSlice* gms_lhs)
00607 {
00608   BLAS_Cpp::trsm(side,tri_rhs.uplo(),trans,tri_rhs.diag(),gms_lhs->rows(),gms_lhs->cols()
00609     ,alpha,tri_rhs.gms().col_ptr(1),tri_rhs.gms().max_rows()
00610     ,gms_lhs->col_ptr(1),gms_lhs->max_rows() );
00611 }
00612 
00613 // implementation:  gms_lhs = alpha * inv(op(tri_rhs)) * op(gms_rhs) (left)
00614 // or       gms_lhs = alpha * op(gms_rhs) * inv(op(tri_rhs)) (right)
00615 // Takes care of temporaries but does not check sizes.
00616 inline
00617 void i_trsm_alt(BLAS_Cpp::Side side, value_type alpha, const DMatrixSliceTri& tri_rhs
00618   , BLAS_Cpp::Transp trans_tri_rhs, const DMatrixSlice& gms_rhs, BLAS_Cpp::Transp trans_gms_rhs
00619   , DMatrixSlice* gms_lhs)
00620 {
00621   assign(gms_lhs,gms_rhs,trans_gms_rhs);
00622   i_trsm( side, trans_tri_rhs, alpha, tri_rhs, gms_lhs );
00623 }
00624 
00625 } // end namespace
00626 
00627 void DenseLinAlgPack::M_StMtM(DMatrix* gm_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
00628   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
00629   , BLAS_Cpp::Transp trans_rhs2)
00630 {
00631   MtM_assert_sizes(  tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1
00632            , gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2       );
00633   gm_lhs->resize(   rows(tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1)
00634           , cols(gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2)        );
00635   i_trmm_alt(BLAS_Cpp::left,alpha,tri_rhs1,trans_rhs1,gms_rhs2,trans_rhs2,&(*gm_lhs)());
00636 }
00637 
00638 void DenseLinAlgPack::M_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
00639   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
00640   , BLAS_Cpp::Transp trans_rhs2)
00641 {
00642   Mp_MtM_assert_sizes(    gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00643               , tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1
00644               , gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2        );
00645   i_trmm_alt(BLAS_Cpp::left,alpha,tri_rhs1,trans_rhs1,gms_rhs2,trans_rhs2,gms_lhs);
00646 }
00647 
00648 void DenseLinAlgPack::M_StMtM(DMatrix* gm_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00649   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceTri& tri_rhs2
00650   , BLAS_Cpp::Transp trans_rhs2)
00651 {
00652   MtM_assert_sizes(  gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1
00653            , tri_rhs2.gms().rows(), tri_rhs2.gms().cols(), trans_rhs2   );
00654   gm_lhs->resize(   rows(gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1)
00655           , cols(tri_rhs2.gms().rows(), tri_rhs2.gms().cols(), trans_rhs2)    );
00656   i_trmm_alt(BLAS_Cpp::right,alpha,tri_rhs2,trans_rhs2,gms_rhs1,trans_rhs1,&(*gm_lhs)());
00657 }
00658 
00659 void DenseLinAlgPack::M_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00660   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceTri& tri_rhs2
00661   , BLAS_Cpp::Transp trans_rhs2)
00662 {
00663   Mp_MtM_assert_sizes(    gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00664               , gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1
00665               , tri_rhs2.gms().rows(), tri_rhs2.gms().cols(), trans_rhs2 );
00666   i_trmm_alt(BLAS_Cpp::right,alpha,tri_rhs2,trans_rhs2,gms_rhs1,trans_rhs1,gms_lhs);
00667 }
00668 
00669 void DenseLinAlgPack::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
00670   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
00671   , BLAS_Cpp::Transp trans_rhs2, value_type beta)
00672 {
00673   // If op(gms_rhs2) == gms_lhs and beta = 0.0 then this is a direct call to the BLAS.
00674   if(    gms_lhs->overlap(gms_rhs2) == SAME_MEM && trans_rhs2 == BLAS_Cpp::no_trans
00675     && beta == 0.0 )
00676   {
00677     i_trmm(BLAS_Cpp::left,trans_rhs1,alpha,tri_rhs1,gms_lhs);
00678   }
00679   else {
00680     // This is something else so the alias problem is not handled.
00681     if(beta != 1.0) Mt_S(gms_lhs,beta);
00682     DMatrix tmp;
00683     M_StMtM(&tmp,alpha,tri_rhs1,trans_rhs1,gms_rhs2,trans_rhs2);
00684     Mp_StM(gms_lhs,1.0,tmp(),BLAS_Cpp::no_trans);
00685   }
00686 }
00687 
00688 void DenseLinAlgPack::Mp_StMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00689   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceTri& tri_rhs2
00690   , BLAS_Cpp::Transp trans_rhs2, value_type beta)
00691 {
00692   // If op(gms_rhs1) == gms_lhs and beta = 0.0 then this is a direct call to the BLAS.
00693   if(    gms_lhs->overlap(gms_rhs1) == SAME_MEM && trans_rhs1 == BLAS_Cpp::no_trans
00694     && beta == 0.0 )
00695   {
00696     i_trmm(BLAS_Cpp::right,trans_rhs2,alpha,tri_rhs2,gms_lhs);
00697   }
00698   else {
00699     // This is something else so the alias problem is not handled.
00700     if(beta != 1.0) Mt_S(gms_lhs,beta);
00701     DMatrix tmp;
00702     M_StMtM(&tmp,alpha,gms_rhs1,trans_rhs1,tri_rhs2,trans_rhs2);
00703     Mp_StM(gms_lhs,1.0,tmp(),BLAS_Cpp::no_trans);
00704   }
00705 }
00706 
00707 void DenseLinAlgPack::M_StInvMtM(DMatrix* gm_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
00708   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
00709   , BLAS_Cpp::Transp trans_rhs2)
00710 {
00711   MtM_assert_sizes(  tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1
00712            , gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2         );
00713   gm_lhs->resize(   rows(tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1)
00714           , cols(gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2)        );
00715   i_trsm_alt(BLAS_Cpp::left,alpha,tri_rhs1,trans_rhs1,gms_rhs2,trans_rhs2,&(*gm_lhs)());
00716 }
00717 
00718 void DenseLinAlgPack::M_StInvMtM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
00719   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSlice& gms_rhs2
00720   , BLAS_Cpp::Transp trans_rhs2)
00721 {
00722   Mp_MtM_assert_sizes(    gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00723               , tri_rhs1.gms().rows(), tri_rhs1.gms().cols(), trans_rhs1
00724               , gms_rhs2.rows(), gms_rhs2.cols(), trans_rhs2 );
00725   i_trsm_alt(BLAS_Cpp::left,alpha,tri_rhs1,trans_rhs1,gms_rhs2,trans_rhs2,gms_lhs);
00726 }
00727 
00728 void DenseLinAlgPack::M_StMtInvM(DMatrix* gm_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00729   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceTri& tri_rhs2
00730   , BLAS_Cpp::Transp trans_rhs2)
00731 {
00732   MtM_assert_sizes(   gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1
00733             , tri_rhs2.gms().rows(), tri_rhs2.gms().cols(), trans_rhs2    );
00734   gm_lhs->resize(   rows(gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1)
00735           , cols(tri_rhs2.gms().rows(), tri_rhs2.gms().cols(), trans_rhs2)    );
00736   i_trsm_alt(BLAS_Cpp::right,alpha,tri_rhs2,trans_rhs2,gms_rhs1,trans_rhs1,&(*gm_lhs)());
00737 }
00738 
00739 void DenseLinAlgPack::M_StMtInvM(DMatrixSlice* gms_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00740   , BLAS_Cpp::Transp trans_rhs1, const DMatrixSliceTri& tri_rhs2
00741   , BLAS_Cpp::Transp trans_rhs2)
00742 {
00743   Mp_MtM_assert_sizes(    gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans
00744               , gms_rhs1.rows(), gms_rhs1.cols(), trans_rhs1
00745               , tri_rhs2.gms().rows(), tri_rhs2.gms().cols(), trans_rhs2 );
00746   i_trsm_alt(BLAS_Cpp::right,alpha,tri_rhs2,trans_rhs2,gms_rhs1,trans_rhs1,gms_lhs);
00747 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines