git.sesse.net Git - casparcg/blob - dependencies/boost/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
  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