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 asum(const V & v)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 nrm2(const V & v)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 amax(const V & v)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 dot(const V1 & v1,const V2 & v2)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> copy(V1 & v1,const V2 & 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> swap(V1 & v1,V2 & 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> scal(V & v,const T & 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> axpy(V1 & v1,const T & t,const V2 & 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> rot(const T1 & t1,V1 & v1,const T2 & t2,V2 & 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> tmv(V & v,const M & 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> tsv(V & v,const M & m,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> gmv(V1 & v1,const T1 & t1,const T2 & t2,const M & m,const V2 & 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> gr(M & m,const T & t,const V1 & v1,const V2 & 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> sr(M & m,const T & t,const V & 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> hr(M & m,const T & t,const V & 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> sr2(M & m,const T & t,const V1 & v1,const V2 & 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> hr2(M & m,const T & t,const V1 & v1,const V2 & 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> tmm(M1 & m1,const T & t,const M2 & m2,const M3 & 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> tsm(M1 & m1,const T & t,const M2 & m2,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> gmm(M1 & m1,const T1 & t1,const T2 & t2,const M2 & m2,const M3 & 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> srk(M1 & m1,const T1 & t1,const T2 & t2,const M2 & 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> hrk(M1 & m1,const T1 & t1,const T2 & t2,const M2 & 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> sr2k(M1 & m1,const T1 & t1,const T2 & t2,const M2 & m2,const M3 & 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> hr2k(M1 & m1,const T1 & t1,const T2 & t2,const M2 & m2,const M3 & 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 500