AbstractLinAlgPack_MatrixNonsingSerial.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: 3/6/00: Provide default implementations for these
00030 // operations.
00031 
00032 #include <assert.h>
00033 
00034 #include "AbstractLinAlgPack_MatrixNonsingSerial.hpp"
00035 #include "AbstractLinAlgPack_MatrixOpSerial.hpp"
00036 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00037 #include "AbstractLinAlgPack_MatrixOpGetGMSMutable.hpp"
00038 #include "AbstractLinAlgPack_MatrixOpGetGMSTri.hpp"
00039 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
00040 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00041 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00042 #include "DenseLinAlgPack_DMatrixClass.hpp"
00043 #include "DenseLinAlgPack_DVectorClass.hpp"
00044 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00045 #include "DenseLinAlgPack_AssertOp.hpp"
00046 
00047 namespace LinAlgOpPack {
00048   using AbstractLinAlgPack::Vp_StV;
00049   using AbstractLinAlgPack::Mp_StM;
00050   using AbstractLinAlgPack::Vp_StMtV;
00051 }
00052 
00053 namespace AbstractLinAlgPack {
00054 
00055 //  Level-2 BLAS
00056 
00057 void MatrixNonsingSerial::V_InvMtV(
00058   DVector* v_lhs, BLAS_Cpp::Transp trans_rhs1,const DVectorSlice& vs_rhs2
00059   ) const
00060 {
00061   const size_type n = rows();
00062   DenseLinAlgPack::MtV_assert_sizes( n, n, trans_rhs1, vs_rhs2.dim() );
00063   v_lhs->resize(n);
00064   this->V_InvMtV( &(*v_lhs)(), trans_rhs1, vs_rhs2 );
00065 }
00066 
00067 void MatrixNonsingSerial::V_InvMtV(
00068   DVector* v_lhs, BLAS_Cpp::Transp trans_rhs1, const SpVectorSlice& sv_rhs2
00069   ) const
00070 {
00071   const size_type n = rows();
00072   DenseLinAlgPack::MtV_assert_sizes( n, n, trans_rhs1, sv_rhs2.dim() );
00073   v_lhs->resize(n);
00074   DVector v_rhs2;
00075   LinAlgOpPack::assign( &v_rhs2, sv_rhs2 );
00076   this->V_InvMtV( &(*v_lhs)(), trans_rhs1, v_rhs2() );
00077 }
00078 
00079 void MatrixNonsingSerial::V_InvMtV(
00080   DVectorSlice* vs_lhs, BLAS_Cpp::Transp trans_rhs1, const SpVectorSlice& sv_rhs2
00081   ) const
00082 {
00083   const size_type n = rows();
00084   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_lhs->dim(), n, n, trans_rhs1, sv_rhs2.dim() );
00085   DVector v_rhs2;
00086   LinAlgOpPack::assign( &v_rhs2, sv_rhs2 );
00087   this->V_InvMtV( vs_lhs, trans_rhs1, v_rhs2() );
00088 }
00089 
00090 value_type MatrixNonsingSerial::transVtInvMtV(
00091   const DVectorSlice& vs_rhs1, BLAS_Cpp::Transp trans_rhs2, const DVectorSlice& vs_rhs3
00092   ) const
00093 {
00094   const size_type n = rows();
00095   DenseLinAlgPack::Vp_MtV_assert_sizes( vs_rhs1.dim(), n, n, trans_rhs2, vs_rhs3.dim() );
00096   DVector tmp;
00097   this->V_InvMtV( &tmp, trans_rhs2, vs_rhs3 );
00098   return DenseLinAlgPack::dot( vs_rhs1, tmp() );
00099 }
00100 
00101 value_type MatrixNonsingSerial::transVtInvMtV(
00102   const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2, const SpVectorSlice& sv_rhs3
00103   ) const
00104 {
00105   const size_type n = rows();
00106   DenseLinAlgPack::Vp_MtV_assert_sizes( sv_rhs1.dim(), n, n, trans_rhs2, sv_rhs3.dim() );
00107   DVector tmp;
00108   this->V_InvMtV( &tmp, trans_rhs2, sv_rhs3 );
00109   return AbstractLinAlgPack::dot( sv_rhs1, tmp() );
00110 }
00111 
00112 // Level-3 BLAS
00113 
00114 void MatrixNonsingSerial::M_StInvMtM(
00115   DMatrix* C, value_type a
00116   ,BLAS_Cpp::Transp A_trans
00117   ,const DMatrixSlice& B, BLAS_Cpp::Transp B_trans
00118   ) const
00119 {
00120   DenseLinAlgPack::MtM_assert_sizes( rows(), cols(), A_trans, B.rows(), B.cols(), B_trans );
00121   C->resize(
00122       BLAS_Cpp::rows( rows(), cols(), A_trans )
00123     , BLAS_Cpp::cols( B.rows(), B.cols(), B_trans )
00124     );
00125   this->M_StInvMtM( &(*C)(), a, A_trans, B, B_trans );
00126 }
00127 
00128 void MatrixNonsingSerial::M_StInvMtM(
00129   DMatrixSlice* C, value_type a
00130   ,BLAS_Cpp::Transp A_trans
00131   ,const DMatrixSlice& B, BLAS_Cpp::Transp B_trans
00132   ) const
00133 {
00134   DenseLinAlgPack::Mp_MtM_assert_sizes( C->rows(), C->cols(), BLAS_Cpp::no_trans
00135     , rows(), cols(), A_trans, B.rows(), B.cols(), B_trans );
00136   //
00137   // C = a * inv(op(A)) * op(B)
00138   //
00139   // C.col(j) = a * inv(op(A)) * op(B).col(j)
00140   //
00141 
00142   for( size_type j = 1; j <= C->cols(); ++j )
00143     AbstractLinAlgPack::V_InvMtV( &C->col(j), *this, A_trans
00144       , DenseLinAlgPack::col( B, B_trans, j ) );
00145   if( a != 1.0 )
00146     LinAlgOpPack::Mt_S( C, a );
00147 }
00148 
00149 void MatrixNonsingSerial::M_StMtInvM(
00150   DMatrix* gm_lhs, value_type alpha
00151   ,const DMatrixSlice& gms_rhs1, BLAS_Cpp::Transp trans_rhs1
00152   ,BLAS_Cpp::Transp trans_rhs2
00153   ) const
00154 {
00155   TEST_FOR_EXCEPT(true);  // ToDo: Implement this!
00156 }
00157 
00158 void MatrixNonsingSerial::M_StMtInvM(
00159   DMatrixSlice* gms_lhs, value_type alpha
00160   ,const DMatrixSlice& gms_rhs1, BLAS_Cpp::Transp trans_rhs1
00161   ,BLAS_Cpp::Transp trans_rhs2
00162   ) const
00163 {
00164   TEST_FOR_EXCEPT(true);  // ToDo: Implement this!
00165 }
00166 
00167 void MatrixNonsingSerial::M_StInvMtM(
00168   DMatrix* C, value_type a
00169   ,BLAS_Cpp::Transp A_trans
00170   ,const MatrixOpSerial& B, BLAS_Cpp::Transp B_trans
00171   ) const
00172 {
00173   DenseLinAlgPack::MtM_assert_sizes( rows(), cols(), A_trans, B.rows(), B.cols(), B_trans );
00174   C->resize(
00175       BLAS_Cpp::rows( rows(), cols(), A_trans )
00176     , BLAS_Cpp::cols( B.rows(), B.cols(), B_trans )
00177     );
00178   AbstractLinAlgPack::M_StInvMtM( &(*C)(), a, *this, A_trans, B, B_trans );
00179 }
00180 
00181 void MatrixNonsingSerial::M_StInvMtM(
00182   DMatrixSlice* C, value_type a
00183   ,BLAS_Cpp::Transp A_trans
00184   ,const MatrixOpSerial& B, BLAS_Cpp::Transp B_trans
00185   ) const
00186 {
00187   using LinAlgOpPack::assign;
00188   DenseLinAlgPack::Mp_MtM_assert_sizes( C->rows(), C->cols(), BLAS_Cpp::no_trans
00189     , rows(), cols(), A_trans, B.rows(), B.cols(), B_trans );
00190   DMatrix B_dense;
00191   assign( &B_dense, B, BLAS_Cpp::no_trans );
00192   AbstractLinAlgPack::M_StInvMtM( C, a, *this, A_trans, B_dense(), B_trans );
00193 }
00194 
00195 void MatrixNonsingSerial::M_StMtInvM(
00196   DMatrix* gm_lhs, value_type alpha
00197   ,const MatrixOpSerial& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00198   ,BLAS_Cpp::Transp trans_rhs2
00199   ) const
00200 {
00201   TEST_FOR_EXCEPT(true);  // ToDo: Implement this!
00202 }
00203 
00204 void MatrixNonsingSerial::M_StMtInvM(
00205   DMatrixSlice* gms_lhs, value_type alpha
00206   ,const MatrixOpSerial& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00207   ,BLAS_Cpp::Transp trans_rhs2
00208   ) const
00209 {
00210   TEST_FOR_EXCEPT(true);  // ToDo: Implement this!
00211 }
00212 
00213 // Overridden from MatrixNonsing
00214 
00215 void MatrixNonsingSerial::V_InvMtV(
00216   VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
00217   ,const Vector& v_rhs2) const
00218 {
00219   VectorDenseMutableEncap       vs_lhs(*v_lhs);
00220   VectorDenseEncap              vs_rhs2(v_rhs2);
00221   this->V_InvMtV( &vs_lhs(), trans_rhs1, vs_rhs2() ); 
00222 }
00223 
00224 void MatrixNonsingSerial::V_InvMtV(
00225   VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
00226   ,const SpVectorSlice& sv_rhs2) const
00227 {
00228   this->V_InvMtV( &VectorDenseMutableEncap(*v_lhs)(), trans_rhs1, sv_rhs2 );
00229 }
00230 
00231 value_type MatrixNonsingSerial::transVtInvMtV(
00232   const Vector& v_rhs1
00233   ,BLAS_Cpp::Transp trans_rhs2, const Vector& v_rhs3) const
00234 {
00235   VectorDenseEncap              vs_rhs1(v_rhs1);
00236   VectorDenseEncap              vs_rhs3(v_rhs3);
00237   return this->transVtInvMtV(vs_rhs1(),trans_rhs2,vs_rhs3());
00238 }
00239 
00240 void MatrixNonsingSerial::M_StInvMtM(
00241   MatrixOp* m_lhs, value_type alpha
00242   ,BLAS_Cpp::Transp trans_rhs1
00243   ,const MatrixOp& mwo_rhs2,BLAS_Cpp::Transp trans_rhs2
00244   ) const
00245 {
00246   using Teuchos::dyn_cast;
00247   MatrixDenseMutableEncap
00248     gms_lhs(m_lhs);      // Warning!  This may throw an exception!
00249   if(const MatrixOpGetGMS* mwo_gms_rhs2 = dynamic_cast<const MatrixOpGetGMS*>(&mwo_rhs2)) {
00250     this->M_StInvMtM(&gms_lhs(),alpha,trans_rhs1,MatrixDenseEncap(*mwo_gms_rhs2)(),trans_rhs2);
00251     return;
00252   }
00253   this->M_StInvMtM(&gms_lhs(),alpha,trans_rhs1,dyn_cast<const MatrixOpSerial>(mwo_rhs2),trans_rhs2);
00254 }
00255 
00256 void MatrixNonsingSerial::M_StMtInvM(
00257   MatrixOp* m_lhs, value_type alpha
00258   ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00259   ,BLAS_Cpp::Transp trans_rhs2
00260   ) const
00261 {
00262   using Teuchos::dyn_cast;
00263   MatrixDenseMutableEncap
00264     gms_lhs(m_lhs);      // Warning!  This may throw an exception!
00265   if(const MatrixOpGetGMS* mwo_gms_rhs1 = dynamic_cast<const MatrixOpGetGMS*>(&mwo_rhs1)) {
00266     this->M_StMtInvM(&gms_lhs(),alpha,MatrixDenseEncap(*mwo_gms_rhs1)(),trans_rhs1,trans_rhs2);
00267     return;
00268   }
00269   this->M_StMtInvM(&gms_lhs(),alpha,dyn_cast<const MatrixOpSerial>(mwo_rhs1),trans_rhs1,trans_rhs2);
00270 }
00271 
00272 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:10:56 2011 for MOOCHO (Single Doxygen Collection) by  doxygen 1.6.3