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