DenseLinAlgPack_DVectorOpBLAS.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 // 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 #include <algorithm>
00030 #include <functional>
00031 #include <math.h> // VC++ 5.0 <cmath> is not CD2 complient yet
00032   // ToDo: Update math function calls to cmath once you get a compiler that meets the standard.
00033 
00034 #include "DenseLinAlgPack_DVectorClass.hpp"
00035 #include "DenseLinAlgPack_DVectorOp.hpp"
00036 #include "DenseLinAlgPack_BLAS_Cpp.hpp"
00037 #include "DenseLinAlgPack_AssertOp.hpp"
00038 
00039 namespace {
00040 
00041 using DenseLinAlgPack::DVector;
00042 using DenseLinAlgPack::DVectorSlice;
00043 typedef DenseLinAlgPack::value_type value_type;
00044 typedef DVectorSlice::size_type size_type;
00045 typedef DVectorSlice::difference_type difference_type;
00046 
00047 //
00048 template< class T >
00049 inline
00050 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00051 //
00052 template< class T >
00053 inline
00054 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00055 
00056 // Utility for copying vector slices.  Takes care of aliasing etc. but not sizes.
00057 void i_assign(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00058   switch(vs_lhs->overlap(vs_rhs)) {
00059     case DenseLinAlgPack::SAME_MEM: return; // assignment to self so skip the copy
00060     case DenseLinAlgPack::SOME_OVERLAP:   // create a temp for the copy
00061     {
00062       DVector temp(vs_rhs);
00063       BLAS_Cpp::copy( temp.dim(), temp.raw_ptr(), 1, vs_lhs->raw_ptr(), vs_lhs->stride() );
00064       return;
00065     }
00066     default:        // no overlap so just perform the copy
00067     {
00068       BLAS_Cpp::copy( vs_rhs.dim(),vs_rhs.raw_ptr(), vs_rhs.stride()
00069         , vs_lhs->raw_ptr(), vs_lhs->stride() );
00070       return;
00071     }
00072   }
00073 }
00074 
00075 inline value_type local_prod( const value_type &v1, const value_type &v2 ) { return v1*v2; }
00076 
00077 } // end namespace
00078 
00079 // ///////////////////////
00080 // op= operations
00081 
00082 void DenseLinAlgPack::Vp_S(DVectorSlice* vs_lhs, value_type alpha) {
00083   if(vs_lhs->stride() == 1) {
00084     DVectorSlice::value_type
00085       *itr    = vs_lhs->start_ptr(),
00086       *itr_end  = itr + vs_lhs->dim();
00087     while(itr != itr_end)
00088       *itr++ += alpha;
00089   }
00090   else {
00091     DVectorSlice::iterator
00092       itr   = vs_lhs->begin(),
00093       itr_end = vs_lhs->end();
00094     while(itr != itr_end)
00095       *itr++ += alpha;
00096   }
00097 }
00098 
00099 void DenseLinAlgPack::Vt_S(DVectorSlice* vs_lhs, value_type alpha) {
00100   if( alpha == 1.0 )
00101     return;
00102   else if( alpha == 0.0 )
00103     *vs_lhs = 0.0;
00104   else
00105     BLAS_Cpp::scal( vs_lhs->dim(), alpha, vs_lhs->raw_ptr(), vs_lhs->stride() );
00106 }
00107 
00108 void DenseLinAlgPack::Vp_StV(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00109   Vp_V_assert_sizes(vs_lhs->dim(), vs_rhs.dim());
00110   BLAS_Cpp::axpy( vs_lhs->dim(), alpha, vs_rhs.raw_ptr(), vs_rhs.stride()
00111     , vs_lhs->raw_ptr(), vs_lhs->stride());
00112 }
00113 
00114 // DVectorSlice as lhs
00115 
00116 void DenseLinAlgPack::assign(DVectorSlice* vs_lhs, value_type alpha) {
00117   if(!vs_lhs->dim())
00118     throw std::length_error("DenseLinAlgPack::assign(vs_lhs,alpha): DVectorSlice must be bound and sized.");
00119   BLAS_Cpp::copy( vs_lhs->dim(), &alpha, 0, vs_lhs->raw_ptr(), vs_lhs->stride() );
00120 }
00121 void DenseLinAlgPack::assign(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00122   Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
00123   i_assign(vs_lhs,vs_rhs);
00124 }
00125 void DenseLinAlgPack::V_VpV(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
00126   VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
00127   Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs1.dim() );
00128   std::transform(vs_rhs1.begin(),vs_rhs1.end(),vs_rhs2.begin(),vs_lhs->begin(),std::plus<value_type>());
00129 }
00130 void DenseLinAlgPack::V_VmV(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
00131   VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
00132   Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs1.dim() );
00133   std::transform(vs_rhs1.begin(),vs_rhs1.end(),vs_rhs2.begin(),vs_lhs->begin(),std::minus<value_type>());
00134 }
00135 void DenseLinAlgPack::V_mV(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00136   Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
00137   (*vs_lhs) = vs_rhs;
00138   BLAS_Cpp::scal( vs_lhs->dim(), -1.0, vs_lhs->raw_ptr(), vs_lhs->stride() );
00139 }
00140 void DenseLinAlgPack::V_StV(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs)
00141 {
00142   Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
00143   (*vs_lhs) = vs_rhs;
00144   BLAS_Cpp::scal( vs_lhs->dim(), alpha, vs_lhs->raw_ptr(), vs_lhs->stride() );
00145 }
00146 
00147 // Elementwise math vector functions. VC++ 5.0 is not allowing use of ptr_fun() with overloaded
00148 // functions so I have to perform the loops straight out.
00149 // ToDo: use ptr_fun() when you get a compiler that works this out.  For now I will just use macros.
00150 #define UNARYOP_VECSLC(LHS, RHS, FUNC)                                    \
00151   DenseLinAlgPack::Vp_V_assert_sizes( (LHS)->dim(), (RHS).dim() );                      \
00152   DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;                   \
00153   for(itr_lhs = LHS->begin(), itr_rhs = RHS.begin(); itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs)     \
00154   { *itr_lhs = FUNC(*itr_rhs); }
00155 
00156 
00157 #define BINARYOP_VECSLC(LHS, RHS1, RHS2, FUNC)                                \
00158   DenseLinAlgPack::VopV_assert_sizes( (RHS1).dim(), (RHS2).dim() );                     \
00159   DenseLinAlgPack::Vp_V_assert_sizes( (LHS)->dim(), (RHS1).dim() );                     \
00160   DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;              \
00161   for(itr_lhs = LHS->begin(), itr_rhs1 = RHS1.begin(), itr_rhs2 = RHS2.begin();             \
00162     itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)                     \
00163   { *itr_lhs = FUNC(*itr_rhs1, *itr_rhs2); }
00164 
00165 #define BINARYOP_BIND1ST_VECSLC(LHS, RHS1, RHS2, FUNC)                            \
00166   DenseLinAlgPack::Vp_V_assert_sizes( (LHS)->dim(), (RHS2).dim() );                     \
00167   DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs2;                  \
00168   for(itr_lhs = LHS->begin(), itr_rhs2 = RHS2.begin();                          \
00169     itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs2)                           \
00170   { *itr_lhs = FUNC(RHS1, *itr_rhs2); }
00171 
00172 #define BINARYOP_BIND2ND_VECSLC(LHS, RHS1, RHS2, FUNC)                            \
00173   DenseLinAlgPack::Vp_V_assert_sizes( (LHS)->dim(), (RHS1).dim());                      \
00174   DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1;                  \
00175   for(itr_lhs = LHS->begin(), itr_rhs1 = RHS1.begin();                          \
00176     itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs1)                           \
00177   { *itr_lhs = FUNC(*itr_rhs1, RHS2); }
00178 
00179 void DenseLinAlgPack::abs(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00180   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::fabs)
00181 }
00182 void DenseLinAlgPack::asin(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00183   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::asin)
00184 }
00185 void DenseLinAlgPack::acos(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00186   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::acos)
00187 }
00188 void DenseLinAlgPack::atan(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00189   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::atan)
00190 }
00191 void DenseLinAlgPack::atan2(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
00192   BINARYOP_VECSLC(vs_lhs, vs_rhs1, vs_rhs2, ::atan2)
00193 }
00194 void DenseLinAlgPack::atan2(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs, value_type alpha) {
00195   BINARYOP_BIND2ND_VECSLC(vs_lhs, vs_rhs, alpha, ::atan2)
00196 }
00197 void DenseLinAlgPack::atan2(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs)
00198 {
00199   BINARYOP_BIND1ST_VECSLC(vs_lhs, alpha, vs_rhs, ::atan2)
00200 }
00201 void DenseLinAlgPack::cos(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00202   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::cos)
00203 }
00204 void DenseLinAlgPack::cosh(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00205   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::cosh)
00206 }
00207 void DenseLinAlgPack::exp(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00208   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::exp)
00209 }
00210 void DenseLinAlgPack::max(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
00211   DenseLinAlgPack::VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
00212   DenseLinAlgPack::Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs1.dim() );
00213   DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;
00214   for(itr_lhs = vs_lhs->begin(), itr_rhs1 = vs_rhs1.begin(), itr_rhs2 = vs_rhs2.begin();
00215     itr_lhs != vs_lhs->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)
00216   { *itr_lhs = my_max(*itr_rhs1, *itr_rhs2); }
00217 }
00218 void DenseLinAlgPack::max(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00219   DenseLinAlgPack::Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
00220   DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;
00221   for(itr_lhs = vs_lhs->begin(), itr_rhs = vs_rhs.begin(); itr_lhs != vs_lhs->end(); ++itr_lhs, ++itr_rhs)
00222   { *itr_lhs = my_max(alpha,*itr_rhs); }
00223 }
00224 void DenseLinAlgPack::min(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
00225   DenseLinAlgPack::VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
00226   DenseLinAlgPack::Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs1.dim() );
00227   DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;
00228   for(itr_lhs = vs_lhs->begin(), itr_rhs1 = vs_rhs1.begin(), itr_rhs2 = vs_rhs2.begin();
00229     itr_lhs != vs_lhs->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)
00230   { *itr_lhs = my_min(*itr_rhs1, *itr_rhs2); }
00231 }
00232 void DenseLinAlgPack::min(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00233   DenseLinAlgPack::Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
00234   DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;
00235   for(itr_lhs = vs_lhs->begin(), itr_rhs = vs_rhs.begin(); itr_lhs != vs_lhs->end(); ++itr_lhs, ++itr_rhs)
00236   { *itr_lhs = my_min(alpha,*itr_rhs); }
00237 }
00238 void DenseLinAlgPack::pow(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
00239   BINARYOP_VECSLC(vs_lhs, vs_rhs1, vs_rhs2, ::pow)
00240 }
00241 void DenseLinAlgPack::pow(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs, value_type alpha) {
00242   BINARYOP_BIND2ND_VECSLC(vs_lhs, vs_rhs, alpha, ::pow)
00243 }
00244 void DenseLinAlgPack::pow(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs, int n) {
00245   BINARYOP_BIND2ND_VECSLC(vs_lhs, vs_rhs, n, ::pow)
00246 }
00247 void DenseLinAlgPack::pow(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00248   BINARYOP_BIND1ST_VECSLC(vs_lhs, alpha, vs_rhs, ::pow)
00249 }
00250 void DenseLinAlgPack::prod( DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2 ) {
00251   BINARYOP_VECSLC(vs_lhs, vs_rhs1, vs_rhs2, local_prod )
00252 }
00253 void DenseLinAlgPack::sqrt(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00254   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::sqrt)
00255 }
00256 void DenseLinAlgPack::sin(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00257   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::sin)
00258 }
00259 void DenseLinAlgPack::sinh(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00260   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::sinh)
00261 }
00262 void DenseLinAlgPack::tan(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00263   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::tan)
00264 }
00265 void DenseLinAlgPack::tanh(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
00266   UNARYOP_VECSLC(vs_lhs, vs_rhs, ::tanh)
00267 }
00268 
00269 // DVector as lhs
00270 
00271 void DenseLinAlgPack::assign(DVector* v_lhs, value_type alpha)
00272 {
00273   if(!v_lhs->dim())
00274     throw std::length_error("DenseLinAlgPack::assign(v_lhs,alpha): DVector must be sized.");
00275   else
00276     BLAS_Cpp::copy( v_lhs->dim(), &alpha, 0, v_lhs->raw_ptr(), v_lhs->stride() );
00277 }
00278 void DenseLinAlgPack::assign(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00279   v_lhs->resize(vs_rhs.dim());
00280   i_assign( &(*v_lhs)(), vs_rhs );
00281 }
00282 void DenseLinAlgPack::V_VpV(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
00283 {
00284   assert_vs_sizes(vs_rhs1.dim(), vs_rhs2.dim());
00285   v_lhs->resize(vs_rhs1.dim());
00286   std::transform(vs_rhs1.begin(),vs_rhs1.end(),vs_rhs2.begin(),v_lhs->begin(),std::plus<value_type>());
00287 //  DVectorSlice::const_iterator
00288 //    v1 = vs_rhs1.begin(),
00289 //    v2 = vs_rhs2.begin();
00290 //  DVector::iterator
00291 //    v = v_lhs->begin(),
00292 //    v_end = v_lhs->end();
00293 //  while( v != v_end )
00294 //    *v++ = *v1++ + *v2++;
00295 }
00296 void DenseLinAlgPack::V_VmV(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
00297 {
00298   assert_vs_sizes(vs_rhs1.dim(), vs_rhs2.dim());
00299   v_lhs->resize(vs_rhs1.dim());
00300   std::transform(vs_rhs1.begin(),vs_rhs1.end(),vs_rhs2.begin(),v_lhs->begin(),std::minus<value_type>());
00301 }
00302 void DenseLinAlgPack::V_mV(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00303   v_lhs->resize(vs_rhs.dim());
00304   (*v_lhs) = vs_rhs;
00305   BLAS_Cpp::scal(v_lhs->dim(), -1.0, v_lhs->raw_ptr(), 1);
00306 }
00307 void DenseLinAlgPack::V_StV(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00308   v_lhs->resize(vs_rhs.dim());
00309   (*v_lhs) = vs_rhs;
00310   BLAS_Cpp::scal( v_lhs->dim(), alpha, v_lhs->raw_ptr(), 1);
00311 }
00312 
00313 void DenseLinAlgPack::rot( const value_type c, const value_type s, DVectorSlice* x, DVectorSlice* y )
00314 {
00315   assert_vs_sizes( x->dim(), y->dim() );
00316   BLAS_Cpp::rot( x->dim(), x->raw_ptr(), x->stride(), y->raw_ptr(), y->stride(), c, s );
00317 }
00318 
00319 // Elementwise math vector functions. VC++ 5.0 is not allowing use of ptr_fun() with overloaded
00320 // functions so I have to perform the loops straight out.
00321 // ToDo: use ptr_fun() when you get a compiler that works this out.  For now I will just use macros.
00322 #define UNARYOP_VEC(LHS, RHS, FUNC)                                     \
00323   LHS->resize(RHS.dim());                                       \
00324   DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;                      \
00325   for(itr_lhs = LHS->begin(), itr_rhs = RHS.begin(); itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs)     \
00326   { *itr_lhs = FUNC(*itr_rhs); }
00327 
00328 #define BINARYOP_VEC(LHS, RHS1, RHS2, FUNC)                                 \
00329   DenseLinAlgPack::assert_vs_sizes(RHS1.dim(), RHS2.dim()); LHS->resize(RHS1.dim());            \
00330   DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;               \
00331   for(itr_lhs = LHS->begin(), itr_rhs1 = RHS1.begin(), itr_rhs2 = RHS2.begin();             \
00332     itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)                     \
00333   { *itr_lhs = FUNC(*itr_rhs1, *itr_rhs2); }
00334 
00335 #define BINARYOP_BIND1ST_VEC(LHS, RHS1, RHS2, FUNC)                             \
00336   LHS->resize(RHS2.dim());                                        \
00337   DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs2;                     \
00338   for(itr_lhs = LHS->begin(), itr_rhs2 = RHS2.begin();                          \
00339     itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs2)                           \
00340   { *itr_lhs = FUNC(RHS1, *itr_rhs2); }
00341 
00342 #define BINARYOP_BIND2ND_VEC(LHS, RHS1, RHS2, FUNC)                             \
00343   LHS->resize(RHS1.dim());                                        \
00344   DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1;                     \
00345   for(itr_lhs = LHS->begin(), itr_rhs1 = RHS1.begin();                          \
00346     itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs1)                           \
00347   { *itr_lhs = FUNC(*itr_rhs1, RHS2); }
00348 
00349 void DenseLinAlgPack::abs(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00350   UNARYOP_VEC(v_lhs, vs_rhs, ::fabs)
00351 }
00352 void DenseLinAlgPack::asin(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00353   UNARYOP_VEC(v_lhs, vs_rhs, ::asin)
00354 }
00355 void DenseLinAlgPack::acos(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00356   UNARYOP_VEC(v_lhs, vs_rhs, ::acos)
00357 }
00358 void DenseLinAlgPack::atan(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00359   UNARYOP_VEC(v_lhs, vs_rhs, ::atan)
00360 }
00361 void DenseLinAlgPack::atan2(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
00362 {
00363   BINARYOP_VEC(v_lhs, vs_rhs1, vs_rhs2, ::atan2)
00364 }
00365 void DenseLinAlgPack::atan2(DVector* v_lhs, const DVectorSlice& vs_rhs, value_type alpha) {
00366   BINARYOP_BIND2ND_VEC(v_lhs, vs_rhs, alpha, ::atan2)
00367 }
00368 void DenseLinAlgPack::atan2(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00369   BINARYOP_BIND1ST_VEC(v_lhs, alpha, vs_rhs, ::atan2)
00370 }
00371 void DenseLinAlgPack::cos(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00372   UNARYOP_VEC(v_lhs, vs_rhs, ::cos)
00373 }
00374 void DenseLinAlgPack::cosh(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00375   UNARYOP_VEC(v_lhs, vs_rhs, ::cosh)
00376 }
00377 void DenseLinAlgPack::exp(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00378   UNARYOP_VEC(v_lhs, vs_rhs, ::exp)
00379 }
00380 void DenseLinAlgPack::max(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
00381 {
00382   DenseLinAlgPack::assert_vs_sizes(vs_rhs1.dim(), vs_rhs2.dim()); v_lhs->resize(vs_rhs1.dim());
00383   DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;
00384   for(itr_lhs = v_lhs->begin(), itr_rhs1 = vs_rhs1.begin(), itr_rhs2 = vs_rhs2.begin();
00385     itr_lhs != v_lhs->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)
00386   { *itr_lhs = my_max(*itr_rhs1, *itr_rhs2); }
00387 }
00388 void DenseLinAlgPack::max(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00389   v_lhs->resize(vs_rhs.dim());
00390   DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;
00391   for(itr_lhs = v_lhs->begin(), itr_rhs = vs_rhs.begin(); itr_lhs != v_lhs->end(); ++itr_lhs, ++itr_rhs)
00392   { *itr_lhs = my_max(alpha,*itr_rhs); }
00393 }
00394 void DenseLinAlgPack::min(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) 
00395 {
00396   DenseLinAlgPack::assert_vs_sizes(vs_rhs1.dim(), vs_rhs2.dim()); v_lhs->resize(vs_rhs1.dim());
00397   DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;
00398   for(itr_lhs = v_lhs->begin(), itr_rhs1 = vs_rhs1.begin(), itr_rhs2 = vs_rhs2.begin();
00399     itr_lhs != v_lhs->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)
00400   { *itr_lhs = my_min(*itr_rhs1, *itr_rhs2); }
00401 }
00402 void DenseLinAlgPack::min(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00403   v_lhs->resize(vs_rhs.dim());
00404   DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;
00405   for(itr_lhs = v_lhs->begin(), itr_rhs = vs_rhs.begin(); itr_lhs != v_lhs->end(); ++itr_lhs, ++itr_rhs)
00406   { *itr_lhs = my_min(alpha,*itr_rhs); }
00407 }
00408 void DenseLinAlgPack::pow(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
00409   BINARYOP_VEC(v_lhs, vs_rhs1, vs_rhs2, ::pow)
00410 }
00411 void DenseLinAlgPack::pow(DVector* v_lhs, const DVectorSlice& vs_rhs, value_type alpha) {
00412   BINARYOP_BIND2ND_VEC(v_lhs, vs_rhs, alpha, ::pow)
00413 }
00414 void DenseLinAlgPack::pow(DVector* v_lhs, const DVectorSlice& vs_rhs, int n) {
00415   BINARYOP_BIND2ND_VEC(v_lhs, vs_rhs, n, ::pow)
00416 }
00417 void DenseLinAlgPack::pow(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
00418   BINARYOP_BIND1ST_VEC(v_lhs, alpha, vs_rhs, ::pow)
00419 }
00420 void DenseLinAlgPack::prod( DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2 ) {
00421   BINARYOP_VEC(v_lhs, vs_rhs1, vs_rhs2, local_prod )
00422 }
00423 void DenseLinAlgPack::sqrt(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00424   UNARYOP_VEC(v_lhs, vs_rhs, ::sqrt)
00425 }
00426 void DenseLinAlgPack::sin(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00427   UNARYOP_VEC(v_lhs, vs_rhs, ::sin)
00428 }
00429 void DenseLinAlgPack::sinh(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00430   UNARYOP_VEC(v_lhs, vs_rhs, ::sinh)
00431 }
00432 void DenseLinAlgPack::tan(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00433   UNARYOP_VEC(v_lhs, vs_rhs, ::tan)
00434 }
00435 void DenseLinAlgPack::tanh(DVector* v_lhs, const DVectorSlice& vs_rhs) {
00436   UNARYOP_VEC(v_lhs, vs_rhs, ::tanh)
00437 }
00438 
00439 // //////////////////////////////
00440 // Scalar returning functions
00441 
00442 DenseLinAlgPack::value_type DenseLinAlgPack::dot(const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
00443 {
00444   VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
00445   return BLAS_Cpp::dot( vs_rhs1.dim(), vs_rhs1.raw_ptr(), vs_rhs1.stride()
00446     , vs_rhs2.raw_ptr(), vs_rhs2.stride() );
00447 }
00448 DenseLinAlgPack::value_type DenseLinAlgPack::max(const DVectorSlice& vs_rhs)
00449 { return *std::max_element(vs_rhs.begin(),vs_rhs.end()); }
00450 DenseLinAlgPack::value_type DenseLinAlgPack::min(const DVectorSlice& vs_rhs)
00451 { return *std::min_element(vs_rhs.begin(),vs_rhs.end()); }
00452 DenseLinAlgPack::value_type DenseLinAlgPack::norm_1(const DVectorSlice& vs_rhs)
00453 { return BLAS_Cpp::asum( vs_rhs.dim(), vs_rhs.raw_ptr(), vs_rhs.stride() ); }
00454 DenseLinAlgPack::value_type DenseLinAlgPack::norm_2(const DVectorSlice& vs_rhs) 
00455 { return BLAS_Cpp::nrm2( vs_rhs.dim(), vs_rhs.raw_ptr(), vs_rhs.stride() ); }
00456 DenseLinAlgPack::value_type DenseLinAlgPack::norm_inf(const DVectorSlice& vs_rhs) {
00457 //  return BLAS_Cpp::iamax( vs_rhs.dim(), vs_rhs.raw_ptr(), vs_rhs.stride() );
00458 //  For some reason iamax() is not working properly?
00459   value_type max_ele = 0.0;
00460   for(DVectorSlice::const_iterator itr = vs_rhs.begin(); itr != vs_rhs.end(); ) {
00461     value_type ele = ::fabs(*itr++);
00462     if(ele > max_ele) max_ele = ele;
00463   }
00464   return max_ele;
00465 }
00466 
00467 // /////////////////////
00468 // Misc operations
00469 
00470 void DenseLinAlgPack::swap( DVectorSlice* vs1, DVectorSlice* vs2 ) {
00471   if( vs1->overlap(*vs2) == SAME_MEM ) return;
00472   VopV_assert_sizes( vs1->dim(), vs2->dim() );
00473   BLAS_Cpp::swap( vs1->dim(), vs1->raw_ptr(), vs1->stride(), vs2->raw_ptr(), vs2->stride() );
00474 }

Generated on Tue Oct 20 12:51:46 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7