AbstractLinAlgPack_MatrixVectorTemplateOpDef.hpp

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

Generated on Wed May 12 21:52:27 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7