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