DenseLinAlgPack: Concreate C++ Classes for Dense Blas-Compatible Linear Algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DenseLinAlgPack_DVectorOpBLAS.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <algorithm>
43 #include <functional>
44 #include <cmath>
45 
46 #include "DenseLinAlgPack_DVectorClass.hpp"
47 #include "DenseLinAlgPack_DVectorOp.hpp"
48 #include "DenseLinAlgPack_BLAS_Cpp.hpp"
49 #include "DenseLinAlgPack_AssertOp.hpp"
50 
51 namespace {
52 
53 using DenseLinAlgPack::DVector;
54 using DenseLinAlgPack::DVectorSlice;
55 typedef DenseLinAlgPack::value_type value_type;
56 typedef DVectorSlice::size_type size_type;
57 typedef DVectorSlice::difference_type difference_type;
58 
59 //
60 template< class T >
61 inline
62 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
63 //
64 template< class T >
65 inline
66 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
67 
68 // Utility for copying vector slices. Takes care of aliasing etc. but not sizes.
69 void i_assign(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
70  switch(vs_lhs->overlap(vs_rhs)) {
71  case DenseLinAlgPack::SAME_MEM: return; // assignment to self so skip the copy
72  case DenseLinAlgPack::SOME_OVERLAP: // create a temp for the copy
73  {
74  DVector temp(vs_rhs);
75  BLAS_Cpp::copy( temp.dim(), temp.raw_ptr(), 1, vs_lhs->raw_ptr(), vs_lhs->stride() );
76  return;
77  }
78  default: // no overlap so just perform the copy
79  {
80  BLAS_Cpp::copy( vs_rhs.dim(),vs_rhs.raw_ptr(), vs_rhs.stride()
81  , vs_lhs->raw_ptr(), vs_lhs->stride() );
82  return;
83  }
84  }
85 }
86 
87 inline value_type local_prod( const value_type &v1, const value_type &v2 ) { return v1*v2; }
88 
89 } // end namespace
90 
91 // ///////////////////////
92 // op= operations
93 
94 void DenseLinAlgPack::Vp_S(DVectorSlice* vs_lhs, value_type alpha) {
95  if(vs_lhs->stride() == 1) {
96  DVectorSlice::value_type
97  *itr = vs_lhs->start_ptr(),
98  *itr_end = itr + vs_lhs->dim();
99  while(itr != itr_end)
100  *itr++ += alpha;
101  }
102  else {
103  DVectorSlice::iterator
104  itr = vs_lhs->begin(),
105  itr_end = vs_lhs->end();
106  while(itr != itr_end)
107  *itr++ += alpha;
108  }
109 }
110 
111 void DenseLinAlgPack::Vt_S(DVectorSlice* vs_lhs, value_type alpha) {
112  if( alpha == 1.0 )
113  return;
114  else if( alpha == 0.0 )
115  *vs_lhs = 0.0;
116  else
117  BLAS_Cpp::scal( vs_lhs->dim(), alpha, vs_lhs->raw_ptr(), vs_lhs->stride() );
118 }
119 
120 void DenseLinAlgPack::Vp_StV(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
121  Vp_V_assert_sizes(vs_lhs->dim(), vs_rhs.dim());
122  BLAS_Cpp::axpy( vs_lhs->dim(), alpha, vs_rhs.raw_ptr(), vs_rhs.stride()
123  , vs_lhs->raw_ptr(), vs_lhs->stride());
124 }
125 
126 // DVectorSlice as lhs
127 
128 void DenseLinAlgPack::assign(DVectorSlice* vs_lhs, value_type alpha) {
129  if(!vs_lhs->dim())
130  throw std::length_error("DenseLinAlgPack::assign(vs_lhs,alpha): DVectorSlice must be bound and sized.");
131  BLAS_Cpp::copy( vs_lhs->dim(), &alpha, 0, vs_lhs->raw_ptr(), vs_lhs->stride() );
132 }
133 void DenseLinAlgPack::assign(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
134  Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
135  i_assign(vs_lhs,vs_rhs);
136 }
137 void DenseLinAlgPack::V_VpV(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
138  VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
139  Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs1.dim() );
140  std::transform(vs_rhs1.begin(),vs_rhs1.end(),vs_rhs2.begin(),vs_lhs->begin(),std::plus<value_type>());
141 }
142 void DenseLinAlgPack::V_VmV(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
143  VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
144  Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs1.dim() );
145  std::transform(vs_rhs1.begin(),vs_rhs1.end(),vs_rhs2.begin(),vs_lhs->begin(),std::minus<value_type>());
146 }
147 void DenseLinAlgPack::V_mV(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
148  Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
149  (*vs_lhs) = vs_rhs;
150  BLAS_Cpp::scal( vs_lhs->dim(), -1.0, vs_lhs->raw_ptr(), vs_lhs->stride() );
151 }
152 void DenseLinAlgPack::V_StV(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs)
153 {
154  Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
155  (*vs_lhs) = vs_rhs;
156  BLAS_Cpp::scal( vs_lhs->dim(), alpha, vs_lhs->raw_ptr(), vs_lhs->stride() );
157 }
158 
159 // Elementwise math vector functions. VC++ 5.0 is not allowing use of ptr_fun() with overloaded
160 // functions so I have to perform the loops straight out.
161 // ToDo: use ptr_fun() when you get a compiler that works this out. For now I will just use macros.
162 #define UNARYOP_VECSLC(LHS, RHS, FUNC) \
163  DenseLinAlgPack::Vp_V_assert_sizes( (LHS)->dim(), (RHS).dim() ); \
164  DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs; \
165  for(itr_lhs = LHS->begin(), itr_rhs = RHS.begin(); itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs) \
166  { *itr_lhs = FUNC(*itr_rhs); }
167 
168 
169 #define BINARYOP_VECSLC(LHS, RHS1, RHS2, FUNC) \
170  DenseLinAlgPack::VopV_assert_sizes( (RHS1).dim(), (RHS2).dim() ); \
171  DenseLinAlgPack::Vp_V_assert_sizes( (LHS)->dim(), (RHS1).dim() ); \
172  DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2; \
173  for(itr_lhs = LHS->begin(), itr_rhs1 = RHS1.begin(), itr_rhs2 = RHS2.begin(); \
174  itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2) \
175  { *itr_lhs = FUNC(*itr_rhs1, *itr_rhs2); }
176 
177 #define BINARYOP_BIND1ST_VECSLC(LHS, RHS1, RHS2, FUNC) \
178  DenseLinAlgPack::Vp_V_assert_sizes( (LHS)->dim(), (RHS2).dim() ); \
179  DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs2; \
180  for(itr_lhs = LHS->begin(), itr_rhs2 = RHS2.begin(); \
181  itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs2) \
182  { *itr_lhs = FUNC(RHS1, *itr_rhs2); }
183 
184 #define BINARYOP_BIND2ND_VECSLC(LHS, RHS1, RHS2, FUNC) \
185  DenseLinAlgPack::Vp_V_assert_sizes( (LHS)->dim(), (RHS1).dim()); \
186  DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1; \
187  for(itr_lhs = LHS->begin(), itr_rhs1 = RHS1.begin(); \
188  itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs1) \
189  { *itr_lhs = FUNC(*itr_rhs1, RHS2); }
190 
191 void DenseLinAlgPack::abs(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
192  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::fabs)
193 }
194 void DenseLinAlgPack::asin(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
195  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::asin)
196 }
197 void DenseLinAlgPack::acos(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
198  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::acos)
199 }
200 void DenseLinAlgPack::atan(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
201  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::atan)
202 }
203 void DenseLinAlgPack::atan2(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
204  BINARYOP_VECSLC(vs_lhs, vs_rhs1, vs_rhs2, ::atan2)
205 }
206 void DenseLinAlgPack::atan2(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs, value_type alpha) {
207  BINARYOP_BIND2ND_VECSLC(vs_lhs, vs_rhs, alpha, ::atan2)
208 }
209 void DenseLinAlgPack::atan2(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs)
210 {
211  BINARYOP_BIND1ST_VECSLC(vs_lhs, alpha, vs_rhs, ::atan2)
212 }
213 void DenseLinAlgPack::cos(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
214  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::cos)
215 }
216 void DenseLinAlgPack::cosh(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
217  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::cosh)
218 }
219 void DenseLinAlgPack::exp(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
220  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::exp)
221 }
222 void DenseLinAlgPack::max(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
223  DenseLinAlgPack::VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
224  DenseLinAlgPack::Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs1.dim() );
225  DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;
226  for(itr_lhs = vs_lhs->begin(), itr_rhs1 = vs_rhs1.begin(), itr_rhs2 = vs_rhs2.begin();
227  itr_lhs != vs_lhs->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)
228  { *itr_lhs = my_max(*itr_rhs1, *itr_rhs2); }
229 }
230 void DenseLinAlgPack::max(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
231  DenseLinAlgPack::Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
232  DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;
233  for(itr_lhs = vs_lhs->begin(), itr_rhs = vs_rhs.begin(); itr_lhs != vs_lhs->end(); ++itr_lhs, ++itr_rhs)
234  { *itr_lhs = my_max(alpha,*itr_rhs); }
235 }
236 void DenseLinAlgPack::min(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
237  DenseLinAlgPack::VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
238  DenseLinAlgPack::Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs1.dim() );
239  DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;
240  for(itr_lhs = vs_lhs->begin(), itr_rhs1 = vs_rhs1.begin(), itr_rhs2 = vs_rhs2.begin();
241  itr_lhs != vs_lhs->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)
242  { *itr_lhs = my_min(*itr_rhs1, *itr_rhs2); }
243 }
244 void DenseLinAlgPack::min(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
245  DenseLinAlgPack::Vp_V_assert_sizes( vs_lhs->dim(), vs_rhs.dim() );
246  DVectorSlice::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;
247  for(itr_lhs = vs_lhs->begin(), itr_rhs = vs_rhs.begin(); itr_lhs != vs_lhs->end(); ++itr_lhs, ++itr_rhs)
248  { *itr_lhs = my_min(alpha,*itr_rhs); }
249 }
250 void DenseLinAlgPack::pow(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
251  BINARYOP_VECSLC(vs_lhs, vs_rhs1, vs_rhs2, std::pow)
252 }
253 void DenseLinAlgPack::pow(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs, value_type alpha) {
254  BINARYOP_BIND2ND_VECSLC(vs_lhs, vs_rhs, alpha, std::pow)
255 }
256 void DenseLinAlgPack::pow(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs, int n) {
257  BINARYOP_BIND2ND_VECSLC(vs_lhs, vs_rhs, n, std::pow)
258 }
259 void DenseLinAlgPack::pow(DVectorSlice* vs_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
260  BINARYOP_BIND1ST_VECSLC(vs_lhs, alpha, vs_rhs, std::pow)
261 }
262 void DenseLinAlgPack::prod( DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2 ) {
263  BINARYOP_VECSLC(vs_lhs, vs_rhs1, vs_rhs2, local_prod )
264 }
265 void DenseLinAlgPack::sqrt(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
266  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::sqrt)
267 }
268 void DenseLinAlgPack::sin(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
269  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::sin)
270 }
271 void DenseLinAlgPack::sinh(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
272  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::sinh)
273 }
274 void DenseLinAlgPack::tan(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
275  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::tan)
276 }
277 void DenseLinAlgPack::tanh(DVectorSlice* vs_lhs, const DVectorSlice& vs_rhs) {
278  UNARYOP_VECSLC(vs_lhs, vs_rhs, ::tanh)
279 }
280 
281 // DVector as lhs
282 
283 void DenseLinAlgPack::assign(DVector* v_lhs, value_type alpha)
284 {
285  if(!v_lhs->dim())
286  throw std::length_error("DenseLinAlgPack::assign(v_lhs,alpha): DVector must be sized.");
287  else
288  BLAS_Cpp::copy( v_lhs->dim(), &alpha, 0, v_lhs->raw_ptr(), v_lhs->stride() );
289 }
290 void DenseLinAlgPack::assign(DVector* v_lhs, const DVectorSlice& vs_rhs) {
291  v_lhs->resize(vs_rhs.dim());
292  i_assign( &(*v_lhs)(), vs_rhs );
293 }
294 void DenseLinAlgPack::V_VpV(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
295 {
296  assert_vs_sizes(vs_rhs1.dim(), vs_rhs2.dim());
297  v_lhs->resize(vs_rhs1.dim());
298  std::transform(vs_rhs1.begin(),vs_rhs1.end(),vs_rhs2.begin(),v_lhs->begin(),std::plus<value_type>());
299 // DVectorSlice::const_iterator
300 // v1 = vs_rhs1.begin(),
301 // v2 = vs_rhs2.begin();
302 // DVector::iterator
303 // v = v_lhs->begin(),
304 // v_end = v_lhs->end();
305 // while( v != v_end )
306 // *v++ = *v1++ + *v2++;
307 }
308 void DenseLinAlgPack::V_VmV(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
309 {
310  assert_vs_sizes(vs_rhs1.dim(), vs_rhs2.dim());
311  v_lhs->resize(vs_rhs1.dim());
312  std::transform(vs_rhs1.begin(),vs_rhs1.end(),vs_rhs2.begin(),v_lhs->begin(),std::minus<value_type>());
313 }
314 void DenseLinAlgPack::V_mV(DVector* v_lhs, const DVectorSlice& vs_rhs) {
315  v_lhs->resize(vs_rhs.dim());
316  (*v_lhs) = vs_rhs;
317  BLAS_Cpp::scal(v_lhs->dim(), -1.0, v_lhs->raw_ptr(), 1);
318 }
319 void DenseLinAlgPack::V_StV(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
320  v_lhs->resize(vs_rhs.dim());
321  (*v_lhs) = vs_rhs;
322  BLAS_Cpp::scal( v_lhs->dim(), alpha, v_lhs->raw_ptr(), 1);
323 }
324 
325 void DenseLinAlgPack::rot( const value_type c, const value_type s, DVectorSlice* x, DVectorSlice* y )
326 {
327  assert_vs_sizes( x->dim(), y->dim() );
328  BLAS_Cpp::rot( x->dim(), x->raw_ptr(), x->stride(), y->raw_ptr(), y->stride(), c, s );
329 }
330 
331 // Elementwise math vector functions. VC++ 5.0 is not allowing use of ptr_fun() with overloaded
332 // functions so I have to perform the loops straight out.
333 // ToDo: use ptr_fun() when you get a compiler that works this out. For now I will just use macros.
334 #define UNARYOP_VEC(LHS, RHS, FUNC) \
335  LHS->resize(RHS.dim()); \
336  DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs; \
337  for(itr_lhs = LHS->begin(), itr_rhs = RHS.begin(); itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs) \
338  { *itr_lhs = FUNC(*itr_rhs); }
339 
340 #define BINARYOP_VEC(LHS, RHS1, RHS2, FUNC) \
341  DenseLinAlgPack::assert_vs_sizes(RHS1.dim(), RHS2.dim()); LHS->resize(RHS1.dim()); \
342  DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2; \
343  for(itr_lhs = LHS->begin(), itr_rhs1 = RHS1.begin(), itr_rhs2 = RHS2.begin(); \
344  itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2) \
345  { *itr_lhs = FUNC(*itr_rhs1, *itr_rhs2); }
346 
347 #define BINARYOP_BIND1ST_VEC(LHS, RHS1, RHS2, FUNC) \
348  LHS->resize(RHS2.dim()); \
349  DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs2; \
350  for(itr_lhs = LHS->begin(), itr_rhs2 = RHS2.begin(); \
351  itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs2) \
352  { *itr_lhs = FUNC(RHS1, *itr_rhs2); }
353 
354 #define BINARYOP_BIND2ND_VEC(LHS, RHS1, RHS2, FUNC) \
355  LHS->resize(RHS1.dim()); \
356  DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1; \
357  for(itr_lhs = LHS->begin(), itr_rhs1 = RHS1.begin(); \
358  itr_lhs != LHS->end(); ++itr_lhs, ++itr_rhs1) \
359  { *itr_lhs = FUNC(*itr_rhs1, RHS2); }
360 
361 void DenseLinAlgPack::abs(DVector* v_lhs, const DVectorSlice& vs_rhs) {
362  UNARYOP_VEC(v_lhs, vs_rhs, ::fabs)
363 }
364 void DenseLinAlgPack::asin(DVector* v_lhs, const DVectorSlice& vs_rhs) {
365  UNARYOP_VEC(v_lhs, vs_rhs, ::asin)
366 }
367 void DenseLinAlgPack::acos(DVector* v_lhs, const DVectorSlice& vs_rhs) {
368  UNARYOP_VEC(v_lhs, vs_rhs, ::acos)
369 }
370 void DenseLinAlgPack::atan(DVector* v_lhs, const DVectorSlice& vs_rhs) {
371  UNARYOP_VEC(v_lhs, vs_rhs, ::atan)
372 }
373 void DenseLinAlgPack::atan2(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
374 {
375  BINARYOP_VEC(v_lhs, vs_rhs1, vs_rhs2, ::atan2)
376 }
377 void DenseLinAlgPack::atan2(DVector* v_lhs, const DVectorSlice& vs_rhs, value_type alpha) {
378  BINARYOP_BIND2ND_VEC(v_lhs, vs_rhs, alpha, ::atan2)
379 }
380 void DenseLinAlgPack::atan2(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
381  BINARYOP_BIND1ST_VEC(v_lhs, alpha, vs_rhs, ::atan2)
382 }
383 void DenseLinAlgPack::cos(DVector* v_lhs, const DVectorSlice& vs_rhs) {
384  UNARYOP_VEC(v_lhs, vs_rhs, ::cos)
385 }
386 void DenseLinAlgPack::cosh(DVector* v_lhs, const DVectorSlice& vs_rhs) {
387  UNARYOP_VEC(v_lhs, vs_rhs, ::cosh)
388 }
389 void DenseLinAlgPack::exp(DVector* v_lhs, const DVectorSlice& vs_rhs) {
390  UNARYOP_VEC(v_lhs, vs_rhs, ::exp)
391 }
392 void DenseLinAlgPack::max(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
393 {
394  DenseLinAlgPack::assert_vs_sizes(vs_rhs1.dim(), vs_rhs2.dim()); v_lhs->resize(vs_rhs1.dim());
395  DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;
396  for(itr_lhs = v_lhs->begin(), itr_rhs1 = vs_rhs1.begin(), itr_rhs2 = vs_rhs2.begin();
397  itr_lhs != v_lhs->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)
398  { *itr_lhs = my_max(*itr_rhs1, *itr_rhs2); }
399 }
400 void DenseLinAlgPack::max(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
401  v_lhs->resize(vs_rhs.dim());
402  DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;
403  for(itr_lhs = v_lhs->begin(), itr_rhs = vs_rhs.begin(); itr_lhs != v_lhs->end(); ++itr_lhs, ++itr_rhs)
404  { *itr_lhs = my_max(alpha,*itr_rhs); }
405 }
406 void DenseLinAlgPack::min(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
407 {
408  DenseLinAlgPack::assert_vs_sizes(vs_rhs1.dim(), vs_rhs2.dim()); v_lhs->resize(vs_rhs1.dim());
409  DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs1, itr_rhs2;
410  for(itr_lhs = v_lhs->begin(), itr_rhs1 = vs_rhs1.begin(), itr_rhs2 = vs_rhs2.begin();
411  itr_lhs != v_lhs->end(); ++itr_lhs, ++itr_rhs1, ++itr_rhs2)
412  { *itr_lhs = my_min(*itr_rhs1, *itr_rhs2); }
413 }
414 void DenseLinAlgPack::min(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
415  v_lhs->resize(vs_rhs.dim());
416  DVector::iterator itr_lhs; DVectorSlice::const_iterator itr_rhs;
417  for(itr_lhs = v_lhs->begin(), itr_rhs = vs_rhs.begin(); itr_lhs != v_lhs->end(); ++itr_lhs, ++itr_rhs)
418  { *itr_lhs = my_min(alpha,*itr_rhs); }
419 }
420 void DenseLinAlgPack::pow(DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2) {
421  BINARYOP_VEC(v_lhs, vs_rhs1, vs_rhs2, std::pow)
422 }
423 void DenseLinAlgPack::pow(DVector* v_lhs, const DVectorSlice& vs_rhs, value_type alpha) {
424  BINARYOP_BIND2ND_VEC(v_lhs, vs_rhs, alpha, std::pow)
425 }
426 void DenseLinAlgPack::pow(DVector* v_lhs, const DVectorSlice& vs_rhs, int n) {
427  BINARYOP_BIND2ND_VEC(v_lhs, vs_rhs, n, std::pow)
428 }
429 void DenseLinAlgPack::pow(DVector* v_lhs, value_type alpha, const DVectorSlice& vs_rhs) {
430  BINARYOP_BIND1ST_VEC(v_lhs, alpha, vs_rhs, std::pow)
431 }
432 void DenseLinAlgPack::prod( DVector* v_lhs, const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2 ) {
433  BINARYOP_VEC(v_lhs, vs_rhs1, vs_rhs2, local_prod )
434 }
435 void DenseLinAlgPack::sqrt(DVector* v_lhs, const DVectorSlice& vs_rhs) {
436  UNARYOP_VEC(v_lhs, vs_rhs, ::sqrt)
437 }
438 void DenseLinAlgPack::sin(DVector* v_lhs, const DVectorSlice& vs_rhs) {
439  UNARYOP_VEC(v_lhs, vs_rhs, ::sin)
440 }
441 void DenseLinAlgPack::sinh(DVector* v_lhs, const DVectorSlice& vs_rhs) {
442  UNARYOP_VEC(v_lhs, vs_rhs, ::sinh)
443 }
444 void DenseLinAlgPack::tan(DVector* v_lhs, const DVectorSlice& vs_rhs) {
445  UNARYOP_VEC(v_lhs, vs_rhs, ::tan)
446 }
447 void DenseLinAlgPack::tanh(DVector* v_lhs, const DVectorSlice& vs_rhs) {
448  UNARYOP_VEC(v_lhs, vs_rhs, ::tanh)
449 }
450 
451 // //////////////////////////////
452 // Scalar returning functions
453 
454 DenseLinAlgPack::value_type DenseLinAlgPack::dot(const DVectorSlice& vs_rhs1, const DVectorSlice& vs_rhs2)
455 {
456  VopV_assert_sizes( vs_rhs1.dim(), vs_rhs2.dim() );
457  return BLAS_Cpp::dot( vs_rhs1.dim(), vs_rhs1.raw_ptr(), vs_rhs1.stride()
458  , vs_rhs2.raw_ptr(), vs_rhs2.stride() );
459 }
460 DenseLinAlgPack::value_type DenseLinAlgPack::max(const DVectorSlice& vs_rhs)
461 { return *std::max_element(vs_rhs.begin(),vs_rhs.end()); }
462 DenseLinAlgPack::value_type DenseLinAlgPack::min(const DVectorSlice& vs_rhs)
463 { return *std::min_element(vs_rhs.begin(),vs_rhs.end()); }
464 DenseLinAlgPack::value_type DenseLinAlgPack::norm_1(const DVectorSlice& vs_rhs)
465 { return BLAS_Cpp::asum( vs_rhs.dim(), vs_rhs.raw_ptr(), vs_rhs.stride() ); }
466 DenseLinAlgPack::value_type DenseLinAlgPack::norm_2(const DVectorSlice& vs_rhs)
467 { return BLAS_Cpp::nrm2( vs_rhs.dim(), vs_rhs.raw_ptr(), vs_rhs.stride() ); }
468 DenseLinAlgPack::value_type DenseLinAlgPack::norm_inf(const DVectorSlice& vs_rhs) {
469 // return BLAS_Cpp::iamax( vs_rhs.dim(), vs_rhs.raw_ptr(), vs_rhs.stride() );
470 // For some reason iamax() is not working properly?
471  value_type max_ele = 0.0;
472  for(DVectorSlice::const_iterator itr = vs_rhs.begin(); itr != vs_rhs.end(); ) {
473  value_type ele = ::fabs(*itr++);
474  if(ele > max_ele) max_ele = ele;
475  }
476  return max_ele;
477 }
478 
479 // /////////////////////
480 // Misc operations
481 
482 void DenseLinAlgPack::swap( DVectorSlice* vs1, DVectorSlice* vs2 ) {
483  if( vs1->overlap(*vs2) == SAME_MEM ) return;
484  VopV_assert_sizes( vs1->dim(), vs2->dim() );
485  BLAS_Cpp::swap( vs1->dim(), vs1->raw_ptr(), vs1->stride(), vs2->raw_ptr(), vs2->stride() );
486 }
size_type
size_t size_type
Generated on Wed Dec 2 2015 08:49:44 for DenseLinAlgPack: Concreate C++ Classes for Dense Blas-Compatible Linear Algebra by   doxygen 1.8.6