EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-01-18 09:43:04
0001 // Copyright (c) 2000-2011 Joerg Walter, Mathias Koch, David Bellot 0002 // 0003 // Distributed under the Boost Software License, Version 1.0. (See 0004 // accompanying file LICENSE_1_0.txt or copy at 0005 // http://www.boost.org/LICENSE_1_0.txt) 0006 // 0007 // The authors gratefully acknowledge the support of 0008 // GeNeSys mbH & Co. KG in producing this work. 0009 0010 #ifndef _BOOST_UBLAS_BLAS_ 0011 #define _BOOST_UBLAS_BLAS_ 0012 0013 #include <boost/numeric/ublas/traits.hpp> 0014 0015 namespace boost { namespace numeric { namespace ublas { 0016 0017 0018 /** Interface and implementation of BLAS level 1 0019 * This includes functions which perform \b vector-vector operations. 0020 * More information about BLAS can be found at 0021 * <a href="http://en.wikipedia.org/wiki/BLAS">http://en.wikipedia.org/wiki/BLAS</a> 0022 */ 0023 namespace blas_1 { 0024 0025 /** 1-Norm: \f$\sum_i |x_i|\f$ (also called \f$\mathcal{L}_1\f$ or Manhattan norm) 0026 * 0027 * \param v a vector or vector expression 0028 * \return the 1-Norm with type of the vector's type 0029 * 0030 * \tparam V type of the vector (not needed by default) 0031 */ 0032 template<class V> 0033 typename type_traits<typename V::value_type>::real_type 0034 asum (const V &v) { 0035 return norm_1 (v); 0036 } 0037 0038 /** 2-Norm: \f$\sum_i |x_i|^2\f$ (also called \f$\mathcal{L}_2\f$ or Euclidean norm) 0039 * 0040 * \param v a vector or vector expression 0041 * \return the 2-Norm with type of the vector's type 0042 * 0043 * \tparam V type of the vector (not needed by default) 0044 */ 0045 template<class V> 0046 typename type_traits<typename V::value_type>::real_type 0047 nrm2 (const V &v) { 0048 return norm_2 (v); 0049 } 0050 0051 /** Infinite-norm: \f$\max_i |x_i|\f$ (also called \f$\mathcal{L}_\infty\f$ norm) 0052 * 0053 * \param v a vector or vector expression 0054 * \return the Infinite-Norm with type of the vector's type 0055 * 0056 * \tparam V type of the vector (not needed by default) 0057 */ 0058 template<class V> 0059 typename type_traits<typename V::value_type>::real_type 0060 amax (const V &v) { 0061 return norm_inf (v); 0062 } 0063 0064 /** Inner product of vectors \f$v_1\f$ and \f$v_2\f$ 0065 * 0066 * \param v1 first vector of the inner product 0067 * \param v2 second vector of the inner product 0068 * \return the inner product of the type of the most generic type of the 2 vectors 0069 * 0070 * \tparam V1 type of first vector (not needed by default) 0071 * \tparam V2 type of second vector (not needed by default) 0072 */ 0073 template<class V1, class V2> 0074 typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type 0075 dot (const V1 &v1, const V2 &v2) { 0076 return inner_prod (v1, v2); 0077 } 0078 0079 /** Copy vector \f$v_2\f$ to \f$v_1\f$ 0080 * 0081 * \param v1 target vector 0082 * \param v2 source vector 0083 * \return a reference to the target vector 0084 * 0085 * \tparam V1 type of first vector (not needed by default) 0086 * \tparam V2 type of second vector (not needed by default) 0087 */ 0088 template<class V1, class V2> 0089 V1 & copy (V1 &v1, const V2 &v2) 0090 { 0091 return v1.assign (v2); 0092 } 0093 0094 /** Swap vectors \f$v_1\f$ and \f$v_2\f$ 0095 * 0096 * \param v1 first vector 0097 * \param v2 second vector 0098 * 0099 * \tparam V1 type of first vector (not needed by default) 0100 * \tparam V2 type of second vector (not needed by default) 0101 */ 0102 template<class V1, class V2> 0103 void swap (V1 &v1, V2 &v2) 0104 { 0105 v1.swap (v2); 0106 } 0107 0108 /** scale vector \f$v\f$ with scalar \f$t\f$ 0109 * 0110 * \param v vector to be scaled 0111 * \param t the scalar 0112 * \return \c t*v 0113 * 0114 * \tparam V type of the vector (not needed by default) 0115 * \tparam T type of the scalar (not needed by default) 0116 */ 0117 template<class V, class T> 0118 V & scal (V &v, const T &t) 0119 { 0120 return v *= t; 0121 } 0122 0123 /** Compute \f$v_1= v_1 + t.v_2\f$ 0124 * 0125 * \param v1 target and first vector 0126 * \param t the scalar 0127 * \param v2 second vector 0128 * \return a reference to the first and target vector 0129 * 0130 * \tparam V1 type of the first vector (not needed by default) 0131 * \tparam T type of the scalar (not needed by default) 0132 * \tparam V2 type of the second vector (not needed by default) 0133 */ 0134 template<class V1, class T, class V2> 0135 V1 & axpy (V1 &v1, const T &t, const V2 &v2) 0136 { 0137 return v1.plus_assign (t * v2); 0138 } 0139 0140 /** Performs rotation of points in the plane and assign the result to the first vector 0141 * 0142 * Each point is defined as a pair \c v1(i) and \c v2(i), being respectively 0143 * the \f$x\f$ and \f$y\f$ coordinates. The parameters \c t1 and \t2 are respectively 0144 * the cosine and sine of the angle of the rotation. 0145 * Results are not returned but directly written into \c v1. 0146 * 0147 * \param t1 cosine of the rotation 0148 * \param v1 vector of \f$x\f$ values 0149 * \param t2 sine of the rotation 0150 * \param v2 vector of \f$y\f$ values 0151 * 0152 * \tparam T1 type of the cosine value (not needed by default) 0153 * \tparam V1 type of the \f$x\f$ vector (not needed by default) 0154 * \tparam T2 type of the sine value (not needed by default) 0155 * \tparam V2 type of the \f$y\f$ vector (not needed by default) 0156 */ 0157 template<class T1, class V1, class T2, class V2> 0158 void rot (const T1 &t1, V1 &v1, const T2 &t2, V2 &v2) 0159 { 0160 typedef typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type promote_type; 0161 vector<promote_type> vt (t1 * v1 + t2 * v2); 0162 v2.assign (- t2 * v1 + t1 * v2); 0163 v1.assign (vt); 0164 } 0165 0166 } 0167 0168 /** \brief Interface and implementation of BLAS level 2 0169 * This includes functions which perform \b matrix-vector operations. 0170 * More information about BLAS can be found at 0171 * <a href="http://en.wikipedia.org/wiki/BLAS">http://en.wikipedia.org/wiki/BLAS</a> 0172 */ 0173 namespace blas_2 { 0174 0175 /** \brief multiply vector \c v with triangular matrix \c m 0176 * 0177 * \param v a vector 0178 * \param m a triangular matrix 0179 * \return the result of the product 0180 * 0181 * \tparam V type of the vector (not needed by default) 0182 * \tparam M type of the matrix (not needed by default) 0183 */ 0184 template<class V, class M> 0185 V & tmv (V &v, const M &m) 0186 { 0187 return v = prod (m, v); 0188 } 0189 0190 /** \brief solve \f$m.x = v\f$ in place, where \c m is a triangular matrix 0191 * 0192 * \param v a vector 0193 * \param m a matrix 0194 * \param C (this parameter is not needed) 0195 * \return a result vector from the above operation 0196 * 0197 * \tparam V type of the vector (not needed by default) 0198 * \tparam M type of the matrix (not needed by default) 0199 * \tparam C n/a 0200 */ 0201 template<class V, class M, class C> 0202 V & tsv (V &v, const M &m, C) 0203 { 0204 return v = solve (m, v, C ()); 0205 } 0206 0207 /** \brief compute \f$ v_1 = t_1.v_1 + t_2.(m.v_2)\f$, a general matrix-vector product 0208 * 0209 * \param v1 a vector 0210 * \param t1 a scalar 0211 * \param t2 another scalar 0212 * \param m a matrix 0213 * \param v2 another vector 0214 * \return the vector \c v1 with the result from the above operation 0215 * 0216 * \tparam V1 type of first vector (not needed by default) 0217 * \tparam T1 type of first scalar (not needed by default) 0218 * \tparam T2 type of second scalar (not needed by default) 0219 * \tparam M type of matrix (not needed by default) 0220 * \tparam V2 type of second vector (not needed by default) 0221 */ 0222 template<class V1, class T1, class T2, class M, class V2> 0223 V1 & gmv (V1 &v1, const T1 &t1, const T2 &t2, const M &m, const V2 &v2) 0224 { 0225 return v1 = t1 * v1 + t2 * prod (m, v2); 0226 } 0227 0228 /** \brief Rank 1 update: \f$ m = m + t.(v_1.v_2^T)\f$ 0229 * 0230 * \param m a matrix 0231 * \param t a scalar 0232 * \param v1 a vector 0233 * \param v2 another vector 0234 * \return a matrix with the result from the above operation 0235 * 0236 * \tparam M type of matrix (not needed by default) 0237 * \tparam T type of scalar (not needed by default) 0238 * \tparam V1 type of first vector (not needed by default) 0239 * \tparam V2type of second vector (not needed by default) 0240 */ 0241 template<class M, class T, class V1, class V2> 0242 M & gr (M &m, const T &t, const V1 &v1, const V2 &v2) 0243 { 0244 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG 0245 return m += t * outer_prod (v1, v2); 0246 #else 0247 return m = m + t * outer_prod (v1, v2); 0248 #endif 0249 } 0250 0251 /** \brief symmetric rank 1 update: \f$m = m + t.(v.v^T)\f$ 0252 * 0253 * \param m a matrix 0254 * \param t a scalar 0255 * \param v a vector 0256 * \return a matrix with the result from the above operation 0257 * 0258 * \tparam M type of matrix (not needed by default) 0259 * \tparam T type of scalar (not needed by default) 0260 * \tparam V type of vector (not needed by default) 0261 */ 0262 template<class M, class T, class V> 0263 M & sr (M &m, const T &t, const V &v) 0264 { 0265 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG 0266 return m += t * outer_prod (v, v); 0267 #else 0268 return m = m + t * outer_prod (v, v); 0269 #endif 0270 } 0271 0272 /** \brief hermitian rank 1 update: \f$m = m + t.(v.v^H)\f$ 0273 * 0274 * \param m a matrix 0275 * \param t a scalar 0276 * \param v a vector 0277 * \return a matrix with the result from the above operation 0278 * 0279 * \tparam M type of matrix (not needed by default) 0280 * \tparam T type of scalar (not needed by default) 0281 * \tparam V type of vector (not needed by default) 0282 */ 0283 template<class M, class T, class V> 0284 M & hr (M &m, const T &t, const V &v) 0285 { 0286 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG 0287 return m += t * outer_prod (v, conj (v)); 0288 #else 0289 return m = m + t * outer_prod (v, conj (v)); 0290 #endif 0291 } 0292 0293 /** \brief symmetric rank 2 update: \f$ m=m+ t.(v_1.v_2^T + v_2.v_1^T)\f$ 0294 * 0295 * \param m a matrix 0296 * \param t a scalar 0297 * \param v1 a vector 0298 * \param v2 another vector 0299 * \return a matrix with the result from the above operation 0300 * 0301 * \tparam M type of matrix (not needed by default) 0302 * \tparam T type of scalar (not needed by default) 0303 * \tparam V1 type of first vector (not needed by default) 0304 * \tparam V2type of second vector (not needed by default) 0305 */ 0306 template<class M, class T, class V1, class V2> 0307 M & sr2 (M &m, const T &t, const V1 &v1, const V2 &v2) 0308 { 0309 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG 0310 return m += t * (outer_prod (v1, v2) + outer_prod (v2, v1)); 0311 #else 0312 return m = m + t * (outer_prod (v1, v2) + outer_prod (v2, v1)); 0313 #endif 0314 } 0315 0316 /** \brief hermitian rank 2 update: \f$m=m+t.(v_1.v_2^H) + v_2.(t.v_1)^H)\f$ 0317 * 0318 * \param m a matrix 0319 * \param t a scalar 0320 * \param v1 a vector 0321 * \param v2 another vector 0322 * \return a matrix with the result from the above operation 0323 * 0324 * \tparam M type of matrix (not needed by default) 0325 * \tparam T type of scalar (not needed by default) 0326 * \tparam V1 type of first vector (not needed by default) 0327 * \tparam V2type of second vector (not needed by default) 0328 */ 0329 template<class M, class T, class V1, class V2> 0330 M & hr2 (M &m, const T &t, const V1 &v1, const V2 &v2) 0331 { 0332 #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG 0333 return m += t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1)); 0334 #else 0335 return m = m + t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1)); 0336 #endif 0337 } 0338 0339 } 0340 0341 /** \brief Interface and implementation of BLAS level 3 0342 * This includes functions which perform \b matrix-matrix operations. 0343 * More information about BLAS can be found at 0344 * <a href="http://en.wikipedia.org/wiki/BLAS">http://en.wikipedia.org/wiki/BLAS</a> 0345 */ 0346 namespace blas_3 { 0347 0348 /** \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 0349 * 0350 * \param m1 a matrix for storing result 0351 * \param t a scalar 0352 * \param m2 a triangular matrix 0353 * \param m3 a triangular matrix 0354 * \return the matrix \c m1 0355 * 0356 * \tparam M1 type of the result matrix (not needed by default) 0357 * \tparam T type of the scalar (not needed by default) 0358 * \tparam M2 type of the first triangular matrix (not needed by default) 0359 * \tparam M3 type of the second triangular matrix (not needed by default) 0360 * 0361 */ 0362 template<class M1, class T, class M2, class M3> 0363 M1 & tmm (M1 &m1, const T &t, const M2 &m2, const M3 &m3) 0364 { 0365 return m1 = t * prod (m2, m3); 0366 } 0367 0368 /** \brief triangular solve \f$ m_2.x = t.m_1\f$ in place, \f$m_2\f$ is a triangular matrix 0369 * 0370 * \param m1 a matrix 0371 * \param t a scalar 0372 * \param m2 a triangular matrix 0373 * \param C (not used) 0374 * \return the \f$m_1\f$ matrix 0375 * 0376 * \tparam M1 type of the first matrix (not needed by default) 0377 * \tparam T type of the scalar (not needed by default) 0378 * \tparam M2 type of the triangular matrix (not needed by default) 0379 * \tparam C (n/a) 0380 */ 0381 template<class M1, class T, class M2, class C> 0382 M1 & tsm (M1 &m1, const T &t, const M2 &m2, C) 0383 { 0384 return m1 = solve (m2, t * m1, C ()); 0385 } 0386 0387 /** \brief general matrix multiplication \f$m_1=t_1.m_1 + t_2.m_2.m_3\f$ 0388 * 0389 * \param m1 first matrix 0390 * \param t1 first scalar 0391 * \param t2 second scalar 0392 * \param m2 second matrix 0393 * \param m3 third matrix 0394 * \return the matrix \c m1 0395 * 0396 * \tparam M1 type of the first matrix (not needed by default) 0397 * \tparam T1 type of the first scalar (not needed by default) 0398 * \tparam T2 type of the second scalar (not needed by default) 0399 * \tparam M2 type of the second matrix (not needed by default) 0400 * \tparam M3 type of the third matrix (not needed by default) 0401 */ 0402 template<class M1, class T1, class T2, class M2, class M3> 0403 M1 & gmm (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) 0404 { 0405 return m1 = t1 * m1 + t2 * prod (m2, m3); 0406 } 0407 0408 /** \brief symmetric rank \a k update: \f$m_1=t.m_1+t_2.(m_2.m_2^T)\f$ 0409 * 0410 * \param m1 first matrix 0411 * \param t1 first scalar 0412 * \param t2 second scalar 0413 * \param m2 second matrix 0414 * \return matrix \c m1 0415 * 0416 * \tparam M1 type of the first matrix (not needed by default) 0417 * \tparam T1 type of the first scalar (not needed by default) 0418 * \tparam T2 type of the second scalar (not needed by default) 0419 * \tparam M2 type of the second matrix (not needed by default) 0420 * \todo use opb_prod() 0421 */ 0422 template<class M1, class T1, class T2, class M2> 0423 M1 & srk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) 0424 { 0425 return m1 = t1 * m1 + t2 * prod (m2, trans (m2)); 0426 } 0427 0428 /** \brief hermitian rank \a k update: \f$m_1=t.m_1+t_2.(m_2.m2^H)\f$ 0429 * 0430 * \param m1 first matrix 0431 * \param t1 first scalar 0432 * \param t2 second scalar 0433 * \param m2 second matrix 0434 * \return matrix \c m1 0435 * 0436 * \tparam M1 type of the first matrix (not needed by default) 0437 * \tparam T1 type of the first scalar (not needed by default) 0438 * \tparam T2 type of the second scalar (not needed by default) 0439 * \tparam M2 type of the second matrix (not needed by default) 0440 * \todo use opb_prod() 0441 */ 0442 template<class M1, class T1, class T2, class M2> 0443 M1 & hrk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) 0444 { 0445 return m1 = t1 * m1 + t2 * prod (m2, herm (m2)); 0446 } 0447 0448 /** \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$ 0449 * 0450 * \param m1 first matrix 0451 * \param t1 first scalar 0452 * \param t2 second scalar 0453 * \param m2 second matrix 0454 * \param m3 third matrix 0455 * \return matrix \c m1 0456 * 0457 * \tparam M1 type of the first matrix (not needed by default) 0458 * \tparam T1 type of the first scalar (not needed by default) 0459 * \tparam T2 type of the second scalar (not needed by default) 0460 * \tparam M2 type of the second matrix (not needed by default) 0461 * \tparam M3 type of the third matrix (not needed by default) 0462 * \todo use opb_prod() 0463 */ 0464 template<class M1, class T1, class T2, class M2, class M3> 0465 M1 & sr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) 0466 { 0467 return m1 = t1 * m1 + t2 * (prod (m2, trans (m3)) + prod (m3, trans (m2))); 0468 } 0469 0470 /** \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$ 0471 * 0472 * \param m1 first matrix 0473 * \param t1 first scalar 0474 * \param t2 second scalar 0475 * \param m2 second matrix 0476 * \param m3 third matrix 0477 * \return matrix \c m1 0478 * 0479 * \tparam M1 type of the first matrix (not needed by default) 0480 * \tparam T1 type of the first scalar (not needed by default) 0481 * \tparam T2 type of the second scalar (not needed by default) 0482 * \tparam M2 type of the second matrix (not needed by default) 0483 * \tparam M3 type of the third matrix (not needed by default) 0484 * \todo use opb_prod() 0485 */ 0486 template<class M1, class T1, class T2, class M2, class M3> 0487 M1 & hr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) 0488 { 0489 return m1 = 0490 t1 * m1 0491 + t2 * prod (m2, herm (m3)) 0492 + type_traits<T2>::conj (t2) * prod (m3, herm (m2)); 0493 } 0494 0495 } 0496 0497 }}} 0498 0499 #endif
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|