AbstractLinAlgPack_SparseVectorOpDef.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 // The declarations for these functions is in the file AbstractLinAlgPack_SparseVectorOpDecl.hpp
00030 // but because of a bug with the MS VC++ 5.0 compiler you can not use
00031 // namespace qualification with definitions of previously declared
00032 // nonmember template funcitons.  By not including the declarations
00033 // and by including this file for automatic instantiation, then
00034 // if the function prototypes are not the same then a compile
00035 // time error is more likely to occur.  Otherwise you could have
00036 // to settle for a compile-time warning that the funciton has
00037 // not been defined or a link-time error that the definition
00038 // could not be found which will be the case when explicit
00039 // instantiation is used.
00040 
00041 // ToDo: 6/9/98 Finish upgrade
00042 
00043 #ifndef SPARSE_VECTOR_OP_DEF_H
00044 #define SPARSE_VECTOR_OP_DEF_H
00045 
00046 #include "AbstractLinAlgPack_Types.hpp"
00047 #include "AbstractLinAlgPack_SparseVectorClass.hpp"
00048 #include "DenseLinAlgPack_DVectorOp.hpp"
00049 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"  // also included in AbstractLinAlgPack_SparseVectorOpDef.hpp
00050 #include "DenseLinAlgPack_DMatrixClass.hpp"
00051 #include "DenseLinAlgPack_AssertOp.hpp"
00052 
00053 namespace {
00054 template< class T >
00055 inline
00056 T my_my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00057 template< class T >
00058 inline
00059 T my_my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00060 } // end namespace
00061 
00062 namespace AbstractLinAlgPack {
00063 
00064 using DenseLinAlgPack::VopV_assert_sizes;
00065 using DenseLinAlgPack::Vp_V_assert_sizes;
00066 using DenseLinAlgPack::Vp_MtV_assert_sizes;
00067 
00068 using DenseLinAlgPack::row;
00069 using DenseLinAlgPack::col;
00070 
00071 namespace SparseVectorUtilityPack {
00072 template<class T_SpVec>
00073 value_type imp_dot2_V_V_SV(const DVectorSlice& vs1, const DVectorSlice& vs2, const T_SpVec& sv);
00074 }
00075 
00076 // result = dot(vs_rhs1,sv_rhs2)
00077 template<class T_SpVec>
00078 value_type dot_V_SV(const DVectorSlice& vs_rhs1, const T_SpVec& sv_rhs2) {
00079   VopV_assert_sizes(vs_rhs1.dim(),sv_rhs2.dim());
00080   value_type result = 0.0;
00081   typename T_SpVec::difference_type offset = sv_rhs2.offset();
00082   for(typename T_SpVec::const_iterator iter = sv_rhs2.begin(); iter != sv_rhs2.end(); ++iter)
00083     result += vs_rhs1(iter->index()+offset) * iter->value();
00084   return result;
00085 }
00086 
00087 // result = dot(sv_rhs1,vs_rhs2).  Just call the above in reverse order
00088 template<class T_SpVec>
00089 value_type dot_SV_V(const T_SpVec& sv_rhs1, const DVectorSlice& vs_rhs2) {
00090   return dot_V_SV(vs_rhs2,sv_rhs1);
00091 }
00092 
00093 // result = ||sv_rhs||1
00094 template<class T_SpVec>
00095 value_type norm_1_SV(const T_SpVec& sv_rhs) {
00096   typename T_SpVec::element_type::value_type result = 0.0;
00097   for(typename T_SpVec::const_iterator iter = sv_rhs.begin(); iter != sv_rhs.end(); ++iter)
00098     result += ::fabs(iter->value());
00099   return result;
00100 }
00101 
00102 // result = ||sv_rhs||2
00103 template<class T_SpVec>
00104 value_type norm_2_SV(const T_SpVec& sv_rhs) {
00105   typename T_SpVec::element_type::value_type result = 0.0;
00106   for(typename T_SpVec::const_iterator iter = sv_rhs.begin(); iter != sv_rhs.end(); ++iter)
00107     result += (iter->value()) * (iter->value());
00108   return result;
00109 }
00110 
00111 // result = ||sv_rhs||inf
00112 template<class T_SpVec>
00113 value_type norm_inf_SV(const T_SpVec& sv_rhs) {
00114   typename T_SpVec::element_type::value_type result = 0.0;
00115   for(typename T_SpVec::const_iterator iter = sv_rhs.begin(); iter != sv_rhs.end(); ++iter)
00116     result = my_my_max(result,::fabs(iter->value()));
00117   return result;
00118 }
00119 
00120 // result = max(sv_rhs)
00121 template<class T_SpVec>
00122 value_type max_SV(const T_SpVec& sv_rhs) {
00123   typename T_SpVec::element_type::value_type result = 0.0;
00124   for(typename T_SpVec::const_iterator iter = sv_rhs.begin(); iter != sv_rhs.end(); ++iter)
00125     result = my_my_max(iter->value(),result);
00126   return result;
00127 }
00128 
00129 // result = min(sv_rhs)
00130 template<class T_SpVec>
00131 value_type min_SV(const T_SpVec& sv_rhs) {
00132   typename T_SpVec::element_type::value_type result = 0.0;
00133   for(typename T_SpVec::const_iterator iter = sv_rhs.begin(); iter != sv_rhs.end(); ++iter)
00134     result = my_my_min(result,iter->value());
00135   return result;
00136 }
00137 
00138 // vs_lhs += alpha * sv_rhs (BLAS xAXPY)
00139 template<class T_SpVec>
00140 void Vt_S( T_SpVec* sv_lhs, value_type alpha )
00141 {
00142   if( alpha == 1.0 ) return;
00143   for(typename T_SpVec::iterator iter = sv_lhs->begin(); iter != sv_lhs->end(); ++iter)
00144     iter->value() *= alpha;
00145 }
00146 
00147 // vs_lhs += alpha * sv_rhs (BLAS xAXPY)
00148 template<class T_SpVec>
00149 void Vp_StSV(DVectorSlice* vs_lhs, value_type alpha, const T_SpVec& sv_rhs)
00150 {
00151   Vp_V_assert_sizes(vs_lhs->dim(),sv_rhs.dim());
00152   typename T_SpVec::difference_type offset = sv_rhs.offset();
00153   for(typename T_SpVec::const_iterator iter = sv_rhs.begin(); iter != sv_rhs.end(); ++iter)
00154     (*vs_lhs)(iter->index() + offset) += alpha * iter->value();
00155 }
00156 
00157 // vs_lhs += alpha * op(gms_rhs1) * sv_rhs2 (BLAS xGEMV) (time = O(sv_rhs2.nz() * vs_lhs.dim())
00158 template<class T_SpVec>
00159 void Vp_StMtSV(DVectorSlice* pvs_lhs, value_type alpha, const DMatrixSlice& gms_rhs1
00160   , BLAS_Cpp::Transp trans_rhs1, const T_SpVec& sv_rhs2)
00161 {
00162 #ifdef _WINDOWS
00163   using DenseLinAlgPack::Vp_StV;  // MS VC++ 6.0 needs help with the name lookups
00164 #endif
00165   DVectorSlice& vs_lhs = *pvs_lhs;
00166 
00167   Vp_MtV_assert_sizes(vs_lhs.dim(),gms_rhs1.rows(),gms_rhs1.cols(),trans_rhs1
00168             , sv_rhs2.dim());
00169   
00170   // Perform the operation by iterating through the sparse vector and performing
00171   // all of the operations on it.
00172   //
00173   // For sparse element e we do the following:
00174   //
00175   // vs_lhs += alpha * e.value() * gms_rhs1.col(e.index());
00176 
00177   typename T_SpVec::difference_type offset = sv_rhs2.offset();
00178 
00179   for(typename T_SpVec::const_iterator sv_rhs2_itr = sv_rhs2.begin(); sv_rhs2_itr != sv_rhs2.end(); ++sv_rhs2_itr)
00180     DenseLinAlgPack::Vp_StV( &vs_lhs, alpha * sv_rhs2_itr->value()
00181             , col( gms_rhs1, trans_rhs1, sv_rhs2_itr->index() + offset ) );
00182 }
00183 
00184 // vs_lhs += alpha * op(tri_rhs1) * sv_rhs2 (BLAS xTRMV)
00185 template<class T_SpVec>
00186 void Vp_StMtSV(DVectorSlice* pvs_lhs, value_type alpha, const DMatrixSliceTri& tri_rhs1
00187   , BLAS_Cpp::Transp trans_rhs1, const T_SpVec& sv_rhs2)
00188 {
00189   DVectorSlice &vs_lhs = *pvs_lhs;
00190 
00191   Vp_MtV_assert_sizes(vs_lhs.dim(),tri_rhs1.rows(),tri_rhs1.cols(),trans_rhs1
00192             , sv_rhs2.dim());
00193 
00194   // Get the effective matrix
00195   BLAS_Cpp::Uplo effective_uplo;
00196   if( (tri_rhs1.uplo() == BLAS_Cpp::lower && trans_rhs1 == BLAS_Cpp::no_trans) ||
00197     (tri_rhs1.uplo() == BLAS_Cpp::upper && trans_rhs1 == BLAS_Cpp::trans)      )
00198   {
00199     effective_uplo = BLAS_Cpp::lower;
00200   }
00201   else {  // must be effective upper
00202     effective_uplo = BLAS_Cpp::upper;
00203   }
00204 
00205   size_type n = tri_rhs1.gms().rows();  // should be same as cols()
00206 
00207   // Implement the operation by looping through the sparse vector only once
00208   // and performing the row operations.  This gives a time = O(n * sv_rhs2.nz())
00209   typename T_SpVec::difference_type offset = sv_rhs2.offset();
00210   for(typename T_SpVec::const_iterator sv_itr = sv_rhs2.begin(); sv_itr != sv_rhs2.end(); ++sv_itr)
00211   {
00212     size_type j = sv_itr->index() + offset;
00213 
00214     // For the nonzero element j = sv_itr->index() we perfom the following
00215     // operations.
00216     //
00217     // Lower:
00218     //    [\]   [\ 0 0 0] [\]
00219     //    [#]  += [\ # 0 0] * [#] jth element
00220     //    [#]   [\ # \ 0] [\]
00221     //    [#]   [\ # \ \] [\]
00222     //          jth
00223     //          col
00224     //
00225     // Upper:
00226     //    [#]   [\ # \ \] [\]
00227     //    [#]  += [0 # \ \] * [#] jth element
00228     //    [\]   [0 0 \ \] [\]
00229     //    [\]   [0 0 0 \] [\]
00230     //          jth
00231     //          col
00232     //
00233     // If we were told that is it is unit diagonal then we will adjust
00234     // accordingly.
00235 
00236     size_type j_adjusted = j; // will be adjusted for unit diagonal
00237     
00238     switch(effective_uplo) {
00239       case BLAS_Cpp::lower: {
00240         if(tri_rhs1.diag() == BLAS_Cpp::unit)
00241         {
00242           // Make the adjustment for unit diaganal
00243           ++j_adjusted;
00244           vs_lhs(j) += alpha * sv_itr->value(); // diagonal element is one
00245         }
00246         // vs_lhs(j,n) = vs_lhs(j,n) + alpha * sv_itr->value() * tri_rhs1.col(j)(j,n)
00247         if(j_adjusted <= n)
00248         {
00249           DenseLinAlgPack::Vp_StV( &vs_lhs(j_adjusted,n), alpha * sv_itr->value()
00250             ,col(tri_rhs1.gms(),trans_rhs1,j)(j_adjusted,n) );
00251         }
00252         break;
00253       }
00254       case BLAS_Cpp::upper: {
00255         if(tri_rhs1.diag() == BLAS_Cpp::unit)
00256         {
00257           // Make the adjustment for unit diaganal
00258           --j_adjusted;
00259           vs_lhs(j) += alpha * sv_itr->value(); // diagonal element is one
00260         }
00261         // vs_lhs(1,j) = vs_lhs(1,j) + alpha * sv_itr->value() * tri_rhs1.col(j)(1,j)
00262         if(j_adjusted > 0)
00263         {
00264           DenseLinAlgPack::Vp_StV( &vs_lhs(1,j_adjusted), alpha * sv_itr->value()
00265             ,col(tri_rhs1.gms(),trans_rhs1,j)(1,j_adjusted) );
00266         }
00267         break;
00268       }
00269     }
00270   }
00271 }
00272 
00273 // vs_lhs += alpha * op(sym_rhs1) * sv_rhs2 (BLAS xSYMV)
00274 template<class T_SpVec>
00275 void Vp_StMtSV(DVectorSlice* pvs_lhs, value_type alpha, const DMatrixSliceSym& sym_rhs1
00276   , BLAS_Cpp::Transp trans_rhs1, const T_SpVec& sv_rhs2)
00277 {
00278   DVectorSlice& vs_lhs = *pvs_lhs;
00279 
00280   Vp_MtV_assert_sizes(vs_lhs.dim(),sym_rhs1.rows(),sym_rhs1.cols(),trans_rhs1
00281             , sv_rhs2.dim());
00282 
00283   size_type size = sv_rhs2.dim();
00284   switch(sym_rhs1.uplo()) {
00285     case BLAS_Cpp::lower: {
00286       DVectorSlice::iterator vs_lhs_itr; size_type i;
00287       for(vs_lhs_itr = vs_lhs.begin(), i = 1; i <= size; ++i)
00288       {
00289         if(i < size) {
00290           *vs_lhs_itr++ +=
00291             alpha *
00292             SparseVectorUtilityPack::imp_dot2_V_V_SV(
00293                      sym_rhs1.gms().row(i)(1,i)
00294                     ,sym_rhs1.gms().col(i)(i+1,size)
00295                     ,sv_rhs2);
00296         }
00297         else
00298           *vs_lhs_itr++ += alpha *
00299             dot_V_SV(sym_rhs1.gms().row(i),sv_rhs2);
00300       }
00301       break;      
00302     }
00303     case BLAS_Cpp::upper: {
00304       DVectorSlice::iterator vs_lhs_itr; size_type i;
00305       for(vs_lhs_itr = vs_lhs.begin(), i = 1; i <= size; ++i)
00306       {
00307         if(i > 1) {
00308           *vs_lhs_itr++ +=
00309             alpha *
00310             SparseVectorUtilityPack::imp_dot2_V_V_SV(
00311                      sym_rhs1.gms().col(i)(1,i-1)
00312                     ,sym_rhs1.gms().row(i)(i,size)
00313                     ,sv_rhs2);
00314         }
00315         else
00316         *vs_lhs_itr++ += alpha * dot_V_SV(sym_rhs1.gms().row(i),sv_rhs2);
00317       }
00318       break;      
00319     }
00320   } 
00321 }
00322 
00323 namespace SparseVectorUtilityPack {
00324 
00325 // Implementation for the product of a concatonated dense vector with a
00326 // sparse vector.  Used for symetric matrix mulitplication.
00327 // In Matlab notation: result = [vs1' , vs2' ] * sv
00328 // where split = vs1.dim(), vs2.dim() == sv.dim() - split
00329 //
00330 // time = O(sv.nz()), space = O(1)
00331 //
00332 template<class T_SpVec>
00333 value_type imp_dot2_V_V_SV(const DVectorSlice& vs1, const DVectorSlice& vs2, const T_SpVec& sv)
00334 {
00335   size_type split = vs1.dim();
00336   value_type result = 0;
00337   typename T_SpVec::difference_type offset = sv.offset();
00338   for(typename T_SpVec::const_iterator sv_itr = sv.begin(); sv_itr != sv.end(); ++sv_itr) {
00339     typename T_SpVec::element_type::indice_type curr_indice = sv_itr->index()+offset;
00340     if(curr_indice <= split)
00341       result += vs1(curr_indice) * sv_itr->value();
00342     else
00343       result += vs2(curr_indice - split) * sv_itr->value();
00344   }
00345   return result;
00346 }
00347 
00348 } // end namespace SparseVectorUtilityPack
00349 
00350 } // end namespace AbstractLinAlgPack
00351 
00352 #endif // SPARSE_VECTOR_OP_DEF_H

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