AbstractLinAlgPack_MatrixNonsing.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_MatrixNonsing.hpp"
00035 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00036 #include "AbstractLinAlgPack_VectorSpace.hpp"
00037 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00038 #include "AbstractLinAlgPack_SpVectorView.hpp"
00039 #include "AbstractLinAlgPack_EtaVector.hpp"
00040 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00041 #include "Teuchos_TestForException.hpp"
00042 #include "Teuchos_dyn_cast.hpp"
00043 
00044 namespace AbstractLinAlgPack {
00045 
00046 // Clone
00047 
00048 MatrixNonsing::mat_mns_mut_ptr_t
00049 MatrixNonsing::clone_mns()
00050 {
00051   return Teuchos::null;
00052 }
00053 
00054 MatrixNonsing::mat_mns_ptr_t
00055 MatrixNonsing::clone_mns() const
00056 {
00057   return const_cast<MatrixNonsing*>(this)->clone_mns(); // Implicit conversion to const
00058 }
00059 
00060 // Level-2 BLAS
00061 
00062 void MatrixNonsing::V_InvMtV(
00063   VectorMutable* y, BLAS_Cpp::Transp M_trans, const SpVectorSlice& sx
00064   ) const
00065 {
00066   if( sx.nz() ) {
00067     VectorSpace::vec_mut_ptr_t
00068       x = (M_trans == BLAS_Cpp::no_trans
00069             ? this->space_cols()
00070             : this->space_rows()
00071         ).create_member();
00072     x->set_sub_vector(sub_vec_view(sx));
00073     this->V_InvMtV(y,M_trans,*x);
00074   }
00075   else {
00076     *y = 0.0;
00077   }
00078 }
00079 
00080 value_type MatrixNonsing::transVtInvMtV(
00081   const Vector& v_rhs1, BLAS_Cpp::Transp trans_rhs2, const Vector& v_rhs3
00082   ) const
00083 {
00084   VectorSpace::vec_mut_ptr_t
00085     v = (trans_rhs2 == BLAS_Cpp::no_trans
00086           ? this->space_rows()
00087           : this->space_cols()
00088       ).create_member();
00089   this->V_InvMtV( v.get(), trans_rhs2, v_rhs3 );
00090   return dot(v_rhs1,*v);
00091 }
00092 
00093 value_type MatrixNonsing::transVtInvMtV(
00094   const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2, const SpVectorSlice& sv_rhs3
00095   ) const
00096 {
00097   VectorSpace::vec_mut_ptr_t
00098     v = (trans_rhs2 == BLAS_Cpp::no_trans
00099           ? this->space_rows()
00100           : this->space_cols()
00101       ).create_member();
00102   this->V_InvMtV( v.get(), trans_rhs2, sv_rhs3 );
00103   return dot(sv_rhs1,*v);
00104 }
00105 
00106 // Level-3 BLAS
00107 
00108 void MatrixNonsing::M_StInvMtM(
00109   MatrixOp* C_lhs, value_type alpha
00110   ,BLAS_Cpp::Transp M_trans
00111   ,const MatrixOp& B, BLAS_Cpp::Transp B_trans
00112   ) const
00113 {
00114   //
00115   // C = a * inv(op(M)) * op(B)
00116   //
00117   using Teuchos::dyn_cast;
00118   using BLAS_Cpp::no_trans;
00119   using BLAS_Cpp::trans;
00120 #ifdef TEUCHOS_DEBUG
00121   TEST_FOR_EXCEPTION(
00122     C_lhs == NULL, std::invalid_argument
00123     ,"MatrixNonsing::M_StInvMtM(...) : Error!" );
00124   
00125 #endif
00126   const size_type
00127     C_rows = C_lhs->rows(),
00128     C_cols = C_lhs->cols();
00129   const size_type
00130     op_B_cols = BLAS_Cpp::cols( B.rows(), B.cols(), B_trans );
00131 #ifdef TEUCHOS_DEBUG
00132   // We can't check vector spaces since *this may not support MatrixOp
00133   // However, we could dynamic cast to see if MatrixOp is supported and then
00134   // be able to use Mp_MtM_assert_compatibility() but this is okay for now.
00135   const size_type
00136     M_rows    = this->rows(),
00137     M_cols    = this->cols(),
00138     op_B_rows = BLAS_Cpp::rows( B.rows(), B.cols(), B_trans );
00139   TEST_FOR_EXCEPTION(
00140     C_rows != M_rows || M_rows != M_cols || M_cols != op_B_rows || C_cols != op_B_cols
00141     , std::invalid_argument
00142     ,"MatrixNonsing::M_StInvMtM(...) : Error!" );
00143 #endif
00144   //
00145   // Compute C = a * inv(op(M)) * op(B) one column at a time:
00146   //
00147   // C(:,j) = inv(op(M)) * a * op(B) * e(j)    , for j = 1...C.cols()
00148   //                       \______________/    
00149   //                              t_j
00150   //
00151   MultiVectorMutable  &C = dyn_cast<MultiVectorMutable>(*C_lhs);
00152   VectorSpace::vec_mut_ptr_t
00153     t_j = ( B_trans == no_trans ? B.space_cols() : B.space_rows() ).create_member();
00154   for( size_type j = 1; j <= C_cols; ++j ) {
00155     // t_j = alpha * op(B) * e_j
00156     EtaVector e_j( j, op_B_cols );
00157     LinAlgOpPack::V_StMtV( t_j.get(), alpha, B, B_trans, e_j() );
00158     // C(:,j) = inv(op(M)) * t_j
00159     AbstractLinAlgPack::V_InvMtV( C.col(j).get(), *this, M_trans, *t_j );
00160   }
00161 }
00162 
00163 void MatrixNonsing::M_StMtInvM(
00164   MatrixOp* g_lhs, value_type alpha
00165   ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00166   ,BLAS_Cpp::Transp trans_rhs2
00167   ) const
00168 {
00169   TEST_FOR_EXCEPT(true); // ToDo: Implement!
00170 }
00171 
00172 } // end namespace AbstractLinAlgPack

Generated on Thu Sep 18 12:35:11 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1