AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_MatrixVectorTemplateOpDef.hpp
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 // Definitions of template functions declared in AbstractLinAlgPack_MatrixVectorTemplateOp.hpp.
00043 
00044 #ifndef MATRIX_VECTOR_TEMPLATE_OP_DEF_H
00045 #define MATRIX_VECTOR_TEMPLATE_OP_DEF_H
00046 
00047 #include "AbstractLinAlgPack_MatrixVectorTemplateOp.hpp"
00048 #include "DenseLinAlgPack_DMatrixClass.hpp"
00049 
00050 // ///////////////////////////////////
00051 // Matrix assignment
00052 
00053 namespace {
00054 // Typedef for vector returning functions (row or col but not diagonal)
00055 typedef AbstractLinAlgPack::DVectorSlice (AbstractLinAlgPack::DMatrixSlice::*Pvec_func)
00056   (AbstractLinAlgPack::DMatrixSlice::size_type);
00057 // Implement sparse matrix, dense matrix assignment.  Sizes are not checked.
00058 template<class T_Matrix>
00059 void imp_assign(AbstractLinAlgPack::DMatrixSlice& gms_lhs, const T_Matrix& gm_rhs
00060   , BLAS_Cpp::Transp trans_rhs)
00061 {
00062   // If trans, then copy col into row, otherwise copy col into col.
00063   Pvec_func vec_func;
00064   if(trans_rhs == BLAS_Cpp::no_trans)   vec_func = &AbstractLinAlgPack::DMatrixSlice::col;
00065   else                  vec_func = &AbstractLinAlgPack::DMatrixSlice::row;
00066   for(int k = 1; k <= gm_rhs.cols(); ++k)
00067     AbstractLinAlgPack::assign((gms_lhs.*vec_func)(k), gm_rhs.col(k));
00068 }
00069 } // end namespace
00070 
00071 // Definitions of template functions for matrix-matrix assignment
00072 
00074 template<class T_Matrix>
00075 void AbstractLinAlgPack::assign(DMatrix& gm_lhs, const T_Matrix& gm_rhs
00076   , BLAS_Cpp::Transp trans_rhs)
00077 {
00078   DenseLinAlgPack::resize_gm_lhs(gm_lhs,gm_rhs.rows(),gm_rhs.cols(),trans_rhs);
00079   DMatrixSlice gms_lhs = gm_lhs;
00080   imp_assign(gms_lhs,gm_rhs,trans_rhs);
00081 }
00082 
00084 template<class T_Matrix>
00085 void AbstractLinAlgPack::assign(DMatrixSlice& gms_lhs, const T_Matrix& gm_rhs
00086   , BLAS_Cpp::Transp trans_rhs)
00087 {
00088   DenseLinAlgPack::assert_gms_lhs(gms_lhs,gm_rhs.rows(),gm_rhs.cols(),trans_rhs);
00089   imp_assign(gms_lhs,gm_rhs,trans_rhs);
00090 }
00091 
00092 // //////////////////////////////////
00093 // Matrix-DVector multiplication
00094 
00095 namespace {
00096 // Throw an exeption of the rhs arguments do not match
00097 template<class T_Matrix>
00098 void imp_assert_V_MtV_rhs_sizes(const T_Matrix& gm_rhs1, BLAS_Cpp::Transp trans_rhs1
00099   , const AbstractLinAlgPack::DVectorSlice& vs_rhs2)
00100 {
00101   typename T_Matrix::size_type
00102     cols = (trans_rhs1 == BLAS_Cpp::no_trans) ? gm_rhs1.cols() : gm_rhs1.rows();
00103 
00104   if(cols != vs_rhs2.size())
00105     throw std::length_error("V_MtV: The sizes of the rhs expression do not match");
00106 }
00107 
00108 // Implementation of matrix-vector multiply (no transpose).  Sizes are not checked.
00109 template<class T_Matrix>
00110 void imp_V_MtV_no_trans(AbstractLinAlgPack::DVectorSlice& vs_lhs, const T_Matrix& gm_rhs1
00111   , const AbstractLinAlgPack::DVectorSlice& vs_rhs2)
00112 {
00113   typedef typename T_Matrix::size_type size_type;
00114   size_type rows = gm_rhs1.rows();
00115   AbstractLinAlgPack::DVectorSlice::iterator itr_v_lhs = vs_lhs.begin();
00116   for(size_type i = 1; i <= rows; ++i)
00117     *itr_v_lhs++ = AbstractLinAlgPack::dot(vs_rhs2,gm_rhs1.row(i));
00118 }
00119 // Implementation of matrix-vector multiply (transpose).  Sizes are not checked.
00120 template<class T_Matrix>
00121 void imp_V_MtV_trans(AbstractLinAlgPack::DVectorSlice& vs_lhs, const T_Matrix& gm_rhs1
00122   , const AbstractLinAlgPack::DVectorSlice& vs_rhs2)
00123 {
00124   typedef typename T_Matrix::size_type size_type;
00125   size_type cols = gm_rhs1.cols();
00126   AbstractLinAlgPack::DVectorSlice::iterator itr_v_lhs = vs_lhs.begin();
00127   for(size_type j = 1; j <= cols; ++j)
00128     *itr_v_lhs++ = AbstractLinAlgPack::dot(vs_rhs2,gm_rhs1.col(j));
00129 }
00130 
00131 } // end namespace
00132 
00133 // Definitions of template functions for matrix-vector multiplication
00134 
00135 template<class T_Matrix>
00136 void AbstractLinAlgPack::V_MtV(DVector& v_lhs, const T_Matrix& gm_rhs1, BLAS_Cpp::Transp trans_rhs1
00137   , const DVectorSlice& vs_rhs2)
00138 {
00139   imp_assert_V_MtV_rhs_sizes(gm_rhs1,trans_rhs1,vs_rhs2);
00140   v_lhs.resize( (trans_rhs1==BLAS_Cpp::no_trans) ? gm_rhs1.rows() : gm_rhs1.cols() );
00141   DVectorSlice vs_lhs = v_lhs;
00142   if(trans_rhs1 == BLAS_Cpp::no_trans)
00143     imp_V_MtV_no_trans(vs_lhs,gm_rhs1,vs_rhs2);
00144   else
00145     imp_V_MtV_trans(vs_lhs,gm_rhs1,vs_rhs2);
00146 }
00147 
00148 template<class T_Matrix>
00149 void AbstractLinAlgPack::V_MtV(DVectorSlice& v_lhs, const T_Matrix& gm_rhs1, BLAS_Cpp::Transp trans_rhs1
00150   , const DVectorSlice& vs_rhs2)
00151 {
00152   imp_assert_V_MtV_rhs_sizes(gm_rhs1,trans_rhs1,vs_rhs2);
00153   DenseLinAlgPack::assert_resize_vs_lhs(v_lhs, (trans_rhs1==BLAS_Cpp::no_trans) ? gm_rhs1.rows() : gm_rhs1.cols());
00154   if(trans_rhs1 == BLAS_Cpp::no_trans)
00155     imp_V_MtV_no_trans(v_lhs,gm_rhs1,vs_rhs2);
00156   else
00157     imp_V_MtV_trans(v_lhs,gm_rhs1,vs_rhs2);
00158 }
00159 
00160 #endif // MATRIX_VECTOR_TEMPLATE_OP_DEF_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends