boost/numeric/ublas/blas.hpp

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