source: XIOS/dev/dev_olga/src/extern/boost/include/boost/numeric/ublas/blas.hpp @ 1022

Last change on this file since 1022 was 1022, checked in by mhnguyen, 7 years ago
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.