boost/numeric/ublas/blas.hpp

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