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