source: NonGTP/Boost/boost/numeric/ublas/blas.hpp @ 857

Revision 857, 10.4 KB checked in by igarcia, 19 years ago (diff)
RevLine 
[857]1//
2//  Copyright (c) 2000-2002
3//  Joerg Walter, Mathias Koch
4//
5//  Permission to use, copy, modify, distribute and sell this software
6//  and its documentation for any purpose is hereby granted without fee,
7//  provided that the above copyright notice appear in all copies and
8//  that both that copyright notice and this permission notice appear
9//  in supporting documentation.  The authors make no representations
10//  about the suitability of this software for any purpose.
11//  It is provided "as is" without express or implied warranty.
12//
13//  The authors gratefully acknowledge the support of
14//  GeNeSys mbH & Co. KG in producing this work.
15//
16
17#ifndef _BOOST_UBLAS_BLAS_
18#define _BOOST_UBLAS_BLAS_
19
20#include <boost/numeric/ublas/traits.hpp>
21
22namespace boost { namespace numeric { namespace ublas {
23
24    namespace blas_1 {
25
26          /** \namespace boost::numeric::ublas::blas_1
27                  \brief wrapper functions for level 1 blas
28          */
29
30
31          /** \brief 1-Norm: \f$\sum_i |x_i|\f$
32                  \ingroup blas1
33           */
34        template<class V>
35        typename type_traits<typename V::value_type>::real_type
36        asum (const V &v) {
37            return norm_1 (v);
38        }
39          /** \brief 2-Norm: \f$\sum_i |x_i|^2\f$
40                  \ingroup blas1
41           */
42        template<class V>
43        typename type_traits<typename V::value_type>::real_type
44        nrm2 (const V &v) {
45            return norm_2 (v);
46        }
47          /** \brief element with larges absolute value: \f$\max_i |x_i|\f$
48                  \ingroup blas1
49          */                 
50        template<class V>
51        typename type_traits<typename V::value_type>::real_type
52        amax (const V &v) {
53            return norm_inf (v);
54        }
55
56          /** \brief inner product of vectors \a v1 and \a v2
57                  \ingroup blas1
58          */                 
59        template<class V1, class V2>
60        typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type
61        dot (const V1 &v1, const V2 &v2) {
62            return inner_prod (v1, v2);
63        }
64
65          /** \brief copy vector \a v2 to \a v1
66                  \ingroup blas1
67          */                 
68        template<class V1, class V2>
69        V1 &
70        copy (V1 &v1, const V2 &v2) {
71            return v1.assign (v2);
72        }
73
74          /** \brief swap vectors \a v1 and \a v2
75                  \ingroup blas1
76          */                 
77        template<class V1, class V2>
78        void swap (V1 &v1, V2 &v2) {
79            v1.swap (v2);
80        }
81
82          /** \brief scale vector \a v with scalar \a t
83                  \ingroup blas1
84          */                 
85        template<class V, class T>
86        V &
87        scal (V &v, const T &t) {
88            return v *= t;
89        }
90
91          /** \brief compute \a v1 = \a v1 + \a t * \a v2
92                  \ingroup blas1
93          */                 
94        template<class V1, class T, class V2>
95        V1 &
96        axpy (V1 &v1, const T &t, const V2 &v2) {
97            return v1.plus_assign (t * v2);
98        }
99
100          /** \brief apply plane rotation
101                  \ingroup blas1
102          */                 
103        template<class T1, class V1, class T2, class V2>
104        void
105        rot (const T1 &t1, V1 &v1, const T2 &t2, V2 &v2) {
106            typedef typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type promote_type;
107            vector<promote_type> vt (t1 * v1 + t2 * v2);
108            v2.assign (- t2 * v1 + t1 * v2);
109            v1.assign (vt);
110        }
111
112    }
113
114    namespace blas_2 {
115
116          /** \namespace boost::numeric::ublas::blas_2
117                  \brief wrapper functions for level 2 blas
118          */
119
120          /** \brief multiply vector \a v with triangular matrix \a m
121                  \ingroup blas2
122                  \todo: check that matrix is really triangular
123          */                 
124        template<class V, class M>
125        V &
126        tmv (V &v, const M &m) {
127            return v = prod (m, v);
128        }
129
130          /** \brief solve \a m \a x = \a v in place, \a m is triangular matrix
131                  \ingroup blas2
132          */                 
133        template<class V, class M, class C>
134        V &
135        tsv (V &v, const M &m, C) {
136            return v = solve (m, v, C ());
137        }
138
139          /** \brief compute \a v1 = \a t1 * \a v1 + \a t2 * (\a m * \a v2)
140                  \ingroup blas2
141          */                 
142        template<class V1, class T1, class T2, class M, class V2>
143        V1 &
144        gmv (V1 &v1, const T1 &t1, const T2 &t2, const M &m, const V2 &v2) {
145            return v1 = t1 * v1 + t2 * prod (m, v2);
146        }
147
148          /** \brief rank 1 update: \a m = \a m + \a t * (\a v1 * \a v2<sup>T</sup>)
149                  \ingroup blas2
150          */                 
151        template<class M, class T, class V1, class V2>
152        M &
153        gr (M &m, const T &t, const V1 &v1, const V2 &v2) {
154#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
155            return m += t * outer_prod (v1, v2);
156#else
157            return m = m + t * outer_prod (v1, v2);
158#endif
159        }
160
161          /** \brief symmetric rank 1 update: \a m = \a m + \a t * (\a v * \a v<sup>T</sup>)
162                  \ingroup blas2
163          */                 
164        template<class M, class T, class V>
165        M &
166        sr (M &m, const T &t, const V &v) {
167#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
168            return m += t * outer_prod (v, v);
169#else
170            return m = m + t * outer_prod (v, v);
171#endif
172        }
173          /** \brief hermitian rank 1 update: \a m = \a m + \a t * (\a v * \a v<sup>H</sup>)
174                  \ingroup blas2
175          */                 
176        template<class M, class T, class V>
177        M &
178        hr (M &m, const T &t, const V &v) {
179#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
180            return m += t * outer_prod (v, conj (v));
181#else
182            return m = m + t * outer_prod (v, conj (v));
183#endif
184        }
185
186          /** \brief symmetric rank 2 update: \a m = \a m + \a t *
187                  (\a v1 * \a v2<sup>T</sup> + \a v2 * \a v1<sup>T</sup>)
188                  \ingroup blas2
189          */                 
190        template<class M, class T, class V1, class V2>
191        M &
192        sr2 (M &m, const T &t, const V1 &v1, const V2 &v2) {
193#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
194            return m += t * (outer_prod (v1, v2) + outer_prod (v2, v1));
195#else
196            return m = m + t * (outer_prod (v1, v2) + outer_prod (v2, v1));
197#endif
198        }
199          /** \brief hermitian rank 2 update: \a m = \a m +
200                  \a t * (\a v1 * \a v2<sup>H</sup>)
201                  + \a v2 * (\a t * \a v1)<sup>H</sup>)
202                  \ingroup blas2
203          */                 
204        template<class M, class T, class V1, class V2>
205        M &
206        hr2 (M &m, const T &t, const V1 &v1, const V2 &v2) {
207#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
208            return m += t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
209#else
210            return m = m + t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
211#endif
212        }
213
214    }
215
216    namespace blas_3 {
217
218          /** \namespace boost::numeric::ublas::blas_3
219                  \brief wrapper functions for level 3 blas
220          */
221
222          /** \brief triangular matrix multiplication
223                  \ingroup blas3
224          */                 
225        template<class M1, class T, class M2, class M3>
226        M1 &
227        tmm (M1 &m1, const T &t, const M2 &m2, const M3 &m3) {
228            return m1 = t * prod (m2, m3);
229        }
230
231          /** \brief triangular solve \a m2 * \a x = \a t * \a m1 in place,
232                  \a m2 is a triangular matrix
233                  \ingroup blas3
234          */                 
235        template<class M1, class T, class M2, class C>
236        M1 &
237        tsm (M1 &m1, const T &t, const M2 &m2, C) {
238            return m1 = solve (m2, t * m1, C ());
239        }
240
241          /** \brief general matrix multiplication
242                  \ingroup blas3
243          */                 
244        template<class M1, class T1, class T2, class M2, class M3>
245        M1 &
246        gmm (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
247            return m1 = t1 * m1 + t2 * prod (m2, m3);
248        }
249
250          /** \brief symmetric rank k update: \a m1 = \a t * \a m1 +
251                  \a t2 * (\a m2 * \a m2<sup>T</sup>)
252                  \ingroup blas3
253                  \todo use opb_prod()
254          */                 
255        template<class M1, class T1, class T2, class M2>
256        M1 &
257        srk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) {
258            return m1 = t1 * m1 + t2 * prod (m2, trans (m2));
259        }
260          /** \brief hermitian rank k update: \a m1 = \a t * \a m1 +
261                  \a t2 * (\a m2 * \a m2<sup>H</sup>)
262                  \ingroup blas3
263                  \todo use opb_prod()
264          */                 
265        template<class M1, class T1, class T2, class M2>
266        M1 &
267        hrk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) {
268            return m1 = t1 * m1 + t2 * prod (m2, herm (m2));
269        }
270
271          /** \brief generalized symmetric rank k update:
272                  \a m1 = \a t1 * \a m1 + \a t2 * (\a m2 * \a m3<sup>T</sup>)
273                  + \a t2 * (\a m3 * \a m2<sup>T</sup>)
274                  \ingroup blas3
275                  \todo use opb_prod()
276          */                 
277        template<class M1, class T1, class T2, class M2, class M3>
278        M1 &
279        sr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
280            return m1 = t1 * m1 + t2 * (prod (m2, trans (m3)) + prod (m3, trans (m2)));
281        }
282          /** \brief generalized hermitian rank k update:
283                  \a m1 = \a t1 * \a m1 + \a t2 * (\a m2 * \a m3<sup>H</sup>)
284                  + (\a m3 * (\a t2 * \a m2)<sup>H</sup>)
285                  \ingroup blas3
286                  \todo use opb_prod()
287          */                 
288        template<class M1, class T1, class T2, class M2, class M3>
289        M1 &
290        hr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
291            return m1 = t1 * m1 + t2 * prod (m2, herm (m3)) + type_traits<T2>::conj (t2) * prod (m3, herm (m2));
292        }
293
294    }
295
296}}}
297
298#endif
299
300
Note: See TracBrowser for help on using the repository browser.