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