New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
blas.hpp in vendors/XIOS/current/extern/boost/include/boost/numeric/ublas – NEMO

source: vendors/XIOS/current/extern/boost/include/boost/numeric/ublas/blas.hpp @ 3408

Last change on this file since 3408 was 3408, checked in by rblod, 12 years ago

importing initial XIOS vendor drop
  • Property svn:keywords set to Id
File size: 16.9 KB
Line 
1//
2//  Copyright (c) 2000-2002
3//  Joerg Walter, Mathias Koch
4//
5//  Distributed under the Boost Software License, Version 1.0. (See
6//  accompanying file LICENSE_1_0.txt or copy at
7//  http://www.boost.org/LICENSE_1_0.txt)
8//
9//  The authors gratefully acknowledge the support of
10//  GeNeSys mbH & Co. KG in producing this work.
11//
12
13#ifndef _BOOST_UBLAS_BLAS_
14#define _BOOST_UBLAS_BLAS_
15
16#include <boost/numeric/ublas/traits.hpp>
17
18namespace boost { namespace numeric { namespace ublas {
19   
20
21    /** Interface and implementation of BLAS level 1
22     * This includes functions which perform \b vector-vector operations.
23     * More information about BLAS can be found at
24     * <a href="http://en.wikipedia.org/wiki/BLAS">http://en.wikipedia.org/wiki/BLAS</a>
25     */
26    namespace blas_1 {
27
28        /** 1-Norm: \f$\sum_i |x_i|\f$ (also called \f$\mathcal{L}_1\f$ or Manhattan norm)
29    *
30    * \param v a vector or vector expression
31    * \return the 1-Norm with type of the vector's type
32    *
33    * \tparam V type of the vector (not needed by default)
34    */
35        template<class V>
36        typename type_traits<typename V::value_type>::real_type
37        asum (const V &v) {
38            return norm_1 (v);
39        }
40
41        /** 2-Norm: \f$\sum_i |x_i|^2\f$ (also called \f$\mathcal{L}_2\f$ or Euclidean norm)
42    *
43    * \param v a vector or vector expression
44    * \return the 2-Norm with type of the vector's type
45    *
46    * \tparam V type of the vector (not needed by default)
47    */
48        template<class V>
49        typename type_traits<typename V::value_type>::real_type
50        nrm2 (const V &v) {
51            return norm_2 (v);
52        }
53
54        /** Infinite-norm: \f$\max_i |x_i|\f$ (also called \f$\mathcal{L}_\infty\f$ norm)
55    *
56    * \param v a vector or vector expression
57    * \return the Infinite-Norm with type of the vector's type
58    *
59    * \tparam V type of the vector (not needed by default)
60    */
61        template<class V>
62        typename type_traits<typename V::value_type>::real_type
63        amax (const V &v) {
64            return norm_inf (v);
65        }
66
67        /** Inner product of vectors \f$v_1\f$ and \f$v_2\f$
68    *
69    * \param v1 first vector of the inner product
70    * \param v2 second vector of the inner product
71    * \return the inner product of the type of the most generic type of the 2 vectors
72    *
73    * \tparam V1 type of first vector (not needed by default)
74    * \tparam V2 type of second vector (not needed by default)
75    */
76        template<class V1, class V2>
77        typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type
78        dot (const V1 &v1, const V2 &v2) {
79            return inner_prod (v1, v2);
80        }
81
82        /** Copy vector \f$v_2\f$ to \f$v_1\f$
83    *
84    * \param v1 target vector
85    * \param v2 source vector
86    * \return a reference to the target vector
87    *
88    * \tparam V1 type of first vector (not needed by default)
89    * \tparam V2 type of second vector (not needed by default)
90    */
91        template<class V1, class V2>
92        V1 & copy (V1 &v1, const V2 &v2) 
93   {
94            return v1.assign (v2);
95        }
96
97        /** Swap vectors \f$v_1\f$ and \f$v_2\f$
98    *
99    * \param v1 first vector
100    * \param v2 second vector
101    *
102         * \tparam V1 type of first vector (not needed by default)
103    * \tparam V2 type of second vector (not needed by default)
104    */
105   template<class V1, class V2>
106        void swap (V1 &v1, V2 &v2) 
107   {
108            v1.swap (v2);
109        }
110
111        /** scale vector \f$v\f$ with scalar \f$t\f$
112    *
113    * \param v vector to be scaled
114    * \param t the scalar
115    * \return \c t*v
116    *
117    * \tparam V type of the vector (not needed by default)
118    * \tparam T type of the scalar (not needed by default)
119    */
120        template<class V, class T>
121        V & scal (V &v, const T &t) 
122   {
123            return v *= t;
124        }
125
126        /** Compute \f$v_1= v_1 +  t.v_2\f$
127    *
128    * \param v1 target and first vector
129    * \param t the scalar
130    * \param v2 second vector
131    * \return a reference to the first and target vector
132    *
133    * \tparam V1 type of the first vector (not needed by default)
134    * \tparam T type of the scalar (not needed by default)
135    * \tparam V2 type of the second vector (not needed by default)
136    */
137        template<class V1, class T, class V2>
138        V1 & axpy (V1 &v1, const T &t, const V2 &v2) 
139   {
140            return v1.plus_assign (t * v2);
141        }
142
143   /** Performs rotation of points in the plane and assign the result to the first vector
144    *
145    * Each point is defined as a pair \c v1(i) and \c v2(i), being respectively
146    * the \f$x\f$ and \f$y\f$ coordinates. The parameters \c t1 and \t2 are respectively
147    * the cosine and sine of the angle of the rotation.
148    * Results are not returned but directly written into \c v1.
149    *
150    * \param t1 cosine of the rotation
151    * \param v1 vector of \f$x\f$ values
152    * \param t2 sine of the rotation
153    * \param v2 vector of \f$y\f$ values
154    *
155    * \tparam T1 type of the cosine value (not needed by default)
156    * \tparam V1 type of the \f$x\f$ vector (not needed by default)
157    * \tparam T2 type of the sine value (not needed by default)
158    * \tparam V2 type of the \f$y\f$ vector (not needed by default)
159    */
160        template<class T1, class V1, class T2, class V2>
161        void rot (const T1 &t1, V1 &v1, const T2 &t2, V2 &v2) 
162   {
163            typedef typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type promote_type;
164            vector<promote_type> vt (t1 * v1 + t2 * v2);
165            v2.assign (- t2 * v1 + t1 * v2);
166            v1.assign (vt);
167        }
168
169    }
170
171    /** \brief Interface and implementation of BLAS level 2
172     * This includes functions which perform \b matrix-vector operations.
173     * More information about BLAS can be found at
174     * <a href="http://en.wikipedia.org/wiki/BLAS">http://en.wikipedia.org/wiki/BLAS</a>
175     */
176    namespace blas_2 {
177
178       /** \brief multiply vector \c v with triangular matrix \c m
179   *
180   * \param v a vector
181   * \param m a triangular matrix
182   * \return the result of the product
183   *
184   * \tparam V type of the vector (not needed by default)
185   * \tparam M type of the matrix (not needed by default)
186        */                 
187        template<class V, class M>
188        V & tmv (V &v, const M &m) 
189   {
190            return v = prod (m, v);
191        }
192
193        /** \brief solve \f$m.x = v\f$ in place, where \c m is a triangular matrix
194    *
195    * \param v a vector
196    * \param m a matrix
197    * \param C (this parameter is not needed)
198    * \return a result vector from the above operation
199    *
200    * \tparam V type of the vector (not needed by default)
201    * \tparam M type of the matrix (not needed by default)
202    * \tparam C n/a
203         */                 
204        template<class V, class M, class C>
205        V & tsv (V &v, const M &m, C) 
206   {
207            return v = solve (m, v, C ());
208        }
209
210        /** \brief compute \f$ v_1 = t_1.v_1 + t_2.(m.v_2)\f$, a general matrix-vector product
211    *
212    * \param v1 a vector
213    * \param t1 a scalar
214    * \param t2 another scalar
215    * \param m a matrix
216    * \param v2 another vector
217    * \return the vector \c v1 with the result from the above operation
218    *
219    * \tparam V1 type of first vector (not needed by default)
220    * \tparam T1 type of first scalar (not needed by default)
221    * \tparam T2 type of second scalar (not needed by default)
222    * \tparam M type of matrix (not needed by default)
223    * \tparam V2 type of second vector (not needed by default)
224         */                 
225        template<class V1, class T1, class T2, class M, class V2>
226        V1 & gmv (V1 &v1, const T1 &t1, const T2 &t2, const M &m, const V2 &v2) 
227   {
228            return v1 = t1 * v1 + t2 * prod (m, v2);
229        }
230
231        /** \brief Rank 1 update: \f$ m = m + t.(v_1.v_2^T)\f$
232    *
233    * \param m a matrix
234    * \param t a scalar
235    * \param v1 a vector
236    * \param v2 another vector
237    * \return a matrix with the result from the above operation
238    *
239    * \tparam M type of matrix (not needed by default)
240    * \tparam T type of scalar (not needed by default)
241    * \tparam V1 type of first vector (not needed by default)
242    * \tparam V2type of second vector (not needed by default)
243    */
244        template<class M, class T, class V1, class V2>
245        M & gr (M &m, const T &t, const V1 &v1, const V2 &v2) 
246   {
247#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
248            return m += t * outer_prod (v1, v2);
249#else
250            return m = m + t * outer_prod (v1, v2);
251#endif
252        }
253
254        /** \brief symmetric rank 1 update: \f$m = m + t.(v.v^T)\f$
255    *
256    * \param m a matrix
257    * \param t a scalar
258    * \param v a vector
259    * \return a matrix with the result from the above operation
260    *
261    * \tparam M type of matrix (not needed by default)
262    * \tparam T type of scalar (not needed by default)
263    * \tparam V type of vector (not needed by default)
264    */
265        template<class M, class T, class V>
266        M & sr (M &m, const T &t, const V &v) 
267   {
268#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
269            return m += t * outer_prod (v, v);
270#else
271            return m = m + t * outer_prod (v, v);
272#endif
273        }
274
275        /** \brief hermitian rank 1 update: \f$m = m + t.(v.v^H)\f$
276    *
277    * \param m a matrix
278    * \param t a scalar
279    * \param v a vector
280    * \return a matrix with the result from the above operation
281    *
282    * \tparam M type of matrix (not needed by default)
283    * \tparam T type of scalar (not needed by default)
284    * \tparam V type of vector (not needed by default)
285    */
286        template<class M, class T, class V>
287        M & hr (M &m, const T &t, const V &v) 
288   {
289#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
290            return m += t * outer_prod (v, conj (v));
291#else
292            return m = m + t * outer_prod (v, conj (v));
293#endif
294        }
295
296         /** \brief symmetric rank 2 update: \f$ m=m+ t.(v_1.v_2^T + v_2.v_1^T)\f$
297     *
298     * \param m a matrix
299     * \param t a scalar
300     * \param v1 a vector
301     * \param v2 another vector
302     * \return a matrix with the result from the above operation
303     *
304     * \tparam M type of matrix (not needed by default)
305     * \tparam T type of scalar (not needed by default)
306     * \tparam V1 type of first vector (not needed by default)
307     * \tparam V2type of second vector (not needed by default)
308          */                 
309        template<class M, class T, class V1, class V2>
310        M & sr2 (M &m, const T &t, const V1 &v1, const V2 &v2) 
311   {
312#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
313            return m += t * (outer_prod (v1, v2) + outer_prod (v2, v1));
314#else
315            return m = m + t * (outer_prod (v1, v2) + outer_prod (v2, v1));
316#endif
317        }
318
319        /** \brief hermitian rank 2 update: \f$m=m+t.(v_1.v_2^H) + v_2.(t.v_1)^H)\f$
320    *
321    * \param m a matrix
322    * \param t a scalar
323    * \param v1 a vector
324    * \param v2 another vector
325    * \return a matrix with the result from the above operation
326    *
327    * \tparam M type of matrix (not needed by default)
328    * \tparam T type of scalar (not needed by default)
329    * \tparam V1 type of first vector (not needed by default)
330    * \tparam V2type of second vector (not needed by default)
331         */                 
332        template<class M, class T, class V1, class V2>
333        M & hr2 (M &m, const T &t, const V1 &v1, const V2 &v2) 
334   {
335#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
336            return m += t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
337#else
338            return m = m + t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
339#endif
340        }
341
342    }
343
344    /** \brief Interface and implementation of BLAS level 3
345     * This includes functions which perform \b matrix-matrix operations.
346     * More information about BLAS can be found at
347     * <a href="http://en.wikipedia.org/wiki/BLAS">http://en.wikipedia.org/wiki/BLAS</a>
348     */
349    namespace blas_3 {
350
351        /** \brief triangular matrix multiplication \f$m_1=t.m_2.m_3\f$ where \f$m_2\f$ and \f$m_3\f$ are triangular
352    *
353    * \param m1 a matrix for storing result
354    * \param t a scalar
355    * \param m2 a triangular matrix
356    * \param m3 a triangular matrix
357    * \return the matrix \c m1
358    *
359    * \tparam M1 type of the result matrix (not needed by default)
360    * \tparam T type of the scalar (not needed by default)
361    * \tparam M2 type of the first triangular matrix (not needed by default)
362    * \tparam M3 type of the second triangular matrix (not needed by default)
363    *
364        */                 
365        template<class M1, class T, class M2, class M3>
366        M1 & tmm (M1 &m1, const T &t, const M2 &m2, const M3 &m3) 
367   {
368            return m1 = t * prod (m2, m3);
369        }
370
371        /** \brief triangular solve \f$ m_2.x = t.m_1\f$ in place, \f$m_2\f$ is a triangular matrix
372    *
373    * \param m1 a matrix
374    * \param t a scalar
375    * \param m2 a triangular matrix
376    * \param C (not used)
377    * \return the \f$m_1\f$ matrix
378    *
379    * \tparam M1 type of the first matrix (not needed by default)
380    * \tparam T type of the scalar (not needed by default)
381    * \tparam M2 type of the triangular matrix (not needed by default)
382    * \tparam C (n/a)
383         */                 
384        template<class M1, class T, class M2, class C>
385        M1 & tsm (M1 &m1, const T &t, const M2 &m2, C) 
386   {
387            return m1 = solve (m2, t * m1, C ());
388        }
389
390        /** \brief general matrix multiplication \f$m_1=t_1.m_1 + t_2.m_2.m_3\f$
391    *
392    * \param m1 first matrix
393    * \param t1 first scalar
394    * \param t2 second scalar
395    * \param m2 second matrix
396    * \param m3 third matrix
397    * \return the matrix \c m1
398    *
399    * \tparam M1 type of the first matrix (not needed by default)
400    * \tparam T1 type of the first scalar (not needed by default)
401    * \tparam T2 type of the second scalar (not needed by default)
402    * \tparam M2 type of the second matrix (not needed by default)
403    * \tparam M3 type of the third matrix (not needed by default)
404         */                 
405        template<class M1, class T1, class T2, class M2, class M3>
406        M1 & gmm (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) 
407   {
408            return m1 = t1 * m1 + t2 * prod (m2, m3);
409        }
410
411        /** \brief symmetric rank \a k update: \f$m_1=t.m_1+t_2.(m_2.m_2^T)\f$
412    *
413    * \param m1 first matrix
414    * \param t1 first scalar
415    * \param t2 second scalar
416    * \param m2 second matrix
417    * \return matrix \c m1
418    *
419    * \tparam M1 type of the first matrix (not needed by default)
420    * \tparam T1 type of the first scalar (not needed by default)
421    * \tparam T2 type of the second scalar (not needed by default)
422    * \tparam M2 type of the second matrix (not needed by default)
423    * \todo use opb_prod()
424         */                 
425        template<class M1, class T1, class T2, class M2>
426        M1 & srk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) 
427   {
428            return m1 = t1 * m1 + t2 * prod (m2, trans (m2));
429        }
430
431        /** \brief hermitian rank \a k update: \f$m_1=t.m_1+t_2.(m_2.m2^H)\f$
432    *
433    * \param m1 first matrix
434    * \param t1 first scalar
435    * \param t2 second scalar
436    * \param m2 second matrix
437    * \return matrix \c m1
438    *
439    * \tparam M1 type of the first matrix (not needed by default)
440    * \tparam T1 type of the first scalar (not needed by default)
441    * \tparam T2 type of the second scalar (not needed by default)
442    * \tparam M2 type of the second matrix (not needed by default)
443    * \todo use opb_prod()
444         */                 
445        template<class M1, class T1, class T2, class M2>
446        M1 & hrk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) 
447   {
448            return m1 = t1 * m1 + t2 * prod (m2, herm (m2));
449        }
450
451        /** \brief generalized symmetric rank \a k update: \f$m_1=t_1.m_1+t_2.(m_2.m3^T)+t_2.(m_3.m2^T)\f$
452    *
453    * \param m1 first matrix
454    * \param t1 first scalar
455    * \param t2 second scalar
456    * \param m2 second matrix
457    * \param m3 third matrix
458    * \return matrix \c m1
459    *
460    * \tparam M1 type of the first matrix (not needed by default)
461    * \tparam T1 type of the first scalar (not needed by default)
462    * \tparam T2 type of the second scalar (not needed by default)
463    * \tparam M2 type of the second matrix (not needed by default)
464    * \tparam M3 type of the third matrix (not needed by default)
465    * \todo use opb_prod()
466         */                 
467        template<class M1, class T1, class T2, class M2, class M3>
468        M1 & sr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) 
469   {
470            return m1 = t1 * m1 + t2 * (prod (m2, trans (m3)) + prod (m3, trans (m2)));
471        }
472
473        /** \brief generalized hermitian rank \a k update: * \f$m_1=t_1.m_1+t_2.(m_2.m_3^H)+(m_3.(t_2.m_2)^H)\f$
474    *
475    * \param m1 first matrix
476    * \param t1 first scalar
477    * \param t2 second scalar
478    * \param m2 second matrix
479    * \param m3 third matrix
480    * \return matrix \c m1
481    *
482    * \tparam M1 type of the first matrix (not needed by default)
483    * \tparam T1 type of the first scalar (not needed by default)
484    * \tparam T2 type of the second scalar (not needed by default)
485    * \tparam M2 type of the second matrix (not needed by default)
486    * \tparam M3 type of the third matrix (not needed by default)
487    * \todo use opb_prod()
488         */                 
489        template<class M1, class T1, class T2, class M2, class M3>
490        M1 & hr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) 
491   {
492            return m1 = 
493            t1 * m1
494          + t2 * prod (m2, herm (m3)) 
495          + type_traits<T2>::conj (t2) * prod (m3, herm (m2));
496        }
497
498    }
499
500}}}
501
502#endif
Note: See TracBrowser for help on using the repository browser.