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