1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/blas.hpp Fri Jun 15 03:10:57 2012 +0200
1.3 @@ -0,0 +1,300 @@
1.4 +//
1.5 +// Copyright (c) 2000-2002
1.6 +// Joerg Walter, Mathias Koch
1.7 +//
1.8 +// Permission to use, copy, modify, distribute and sell this software
1.9 +// and its documentation for any purpose is hereby granted without fee,
1.10 +// provided that the above copyright notice appear in all copies and
1.11 +// that both that copyright notice and this permission notice appear
1.12 +// in supporting documentation. The authors make no representations
1.13 +// about the suitability of this software for any purpose.
1.14 +// It is provided "as is" without express or implied warranty.
1.15 +//
1.16 +// The authors gratefully acknowledge the support of
1.17 +// GeNeSys mbH & Co. KG in producing this work.
1.18 +//
1.19 +
1.20 +#ifndef _BOOST_UBLAS_BLAS_
1.21 +#define _BOOST_UBLAS_BLAS_
1.22 +
1.23 +#include <boost/numeric/ublas/traits.hpp>
1.24 +
1.25 +namespace boost { namespace numeric { namespace ublas {
1.26 +
1.27 + namespace blas_1 {
1.28 +
1.29 + /** \namespace boost::numeric::ublas::blas_1
1.30 + \brief wrapper functions for level 1 blas
1.31 + */
1.32 +
1.33 +
1.34 + /** \brief 1-Norm: \f$\sum_i |x_i|\f$
1.35 + \ingroup blas1
1.36 + */
1.37 + template<class V>
1.38 + typename type_traits<typename V::value_type>::real_type
1.39 + asum (const V &v) {
1.40 + return norm_1 (v);
1.41 + }
1.42 + /** \brief 2-Norm: \f$\sum_i |x_i|^2\f$
1.43 + \ingroup blas1
1.44 + */
1.45 + template<class V>
1.46 + typename type_traits<typename V::value_type>::real_type
1.47 + nrm2 (const V &v) {
1.48 + return norm_2 (v);
1.49 + }
1.50 + /** \brief element with larges absolute value: \f$\max_i |x_i|\f$
1.51 + \ingroup blas1
1.52 + */
1.53 + template<class V>
1.54 + typename type_traits<typename V::value_type>::real_type
1.55 + amax (const V &v) {
1.56 + return norm_inf (v);
1.57 + }
1.58 +
1.59 + /** \brief inner product of vectors \a v1 and \a v2
1.60 + \ingroup blas1
1.61 + */
1.62 + template<class V1, class V2>
1.63 + typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type
1.64 + dot (const V1 &v1, const V2 &v2) {
1.65 + return inner_prod (v1, v2);
1.66 + }
1.67 +
1.68 + /** \brief copy vector \a v2 to \a v1
1.69 + \ingroup blas1
1.70 + */
1.71 + template<class V1, class V2>
1.72 + V1 &
1.73 + copy (V1 &v1, const V2 &v2) {
1.74 + return v1.assign (v2);
1.75 + }
1.76 +
1.77 + /** \brief swap vectors \a v1 and \a v2
1.78 + \ingroup blas1
1.79 + */
1.80 + template<class V1, class V2>
1.81 + void swap (V1 &v1, V2 &v2) {
1.82 + v1.swap (v2);
1.83 + }
1.84 +
1.85 + /** \brief scale vector \a v with scalar \a t
1.86 + \ingroup blas1
1.87 + */
1.88 + template<class V, class T>
1.89 + V &
1.90 + scal (V &v, const T &t) {
1.91 + return v *= t;
1.92 + }
1.93 +
1.94 + /** \brief compute \a v1 = \a v1 + \a t * \a v2
1.95 + \ingroup blas1
1.96 + */
1.97 + template<class V1, class T, class V2>
1.98 + V1 &
1.99 + axpy (V1 &v1, const T &t, const V2 &v2) {
1.100 + return v1.plus_assign (t * v2);
1.101 + }
1.102 +
1.103 + /** \brief apply plane rotation
1.104 + \ingroup blas1
1.105 + */
1.106 + template<class T1, class V1, class T2, class V2>
1.107 + void
1.108 + rot (const T1 &t1, V1 &v1, const T2 &t2, V2 &v2) {
1.109 + typedef typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type promote_type;
1.110 + vector<promote_type> vt (t1 * v1 + t2 * v2);
1.111 + v2.assign (- t2 * v1 + t1 * v2);
1.112 + v1.assign (vt);
1.113 + }
1.114 +
1.115 + }
1.116 +
1.117 + namespace blas_2 {
1.118 +
1.119 + /** \namespace boost::numeric::ublas::blas_2
1.120 + \brief wrapper functions for level 2 blas
1.121 + */
1.122 +
1.123 + /** \brief multiply vector \a v with triangular matrix \a m
1.124 + \ingroup blas2
1.125 + \todo: check that matrix is really triangular
1.126 + */
1.127 + template<class V, class M>
1.128 + V &
1.129 + tmv (V &v, const M &m) {
1.130 + return v = prod (m, v);
1.131 + }
1.132 +
1.133 + /** \brief solve \a m \a x = \a v in place, \a m is triangular matrix
1.134 + \ingroup blas2
1.135 + */
1.136 + template<class V, class M, class C>
1.137 + V &
1.138 + tsv (V &v, const M &m, C) {
1.139 + return v = solve (m, v, C ());
1.140 + }
1.141 +
1.142 + /** \brief compute \a v1 = \a t1 * \a v1 + \a t2 * (\a m * \a v2)
1.143 + \ingroup blas2
1.144 + */
1.145 + template<class V1, class T1, class T2, class M, class V2>
1.146 + V1 &
1.147 + gmv (V1 &v1, const T1 &t1, const T2 &t2, const M &m, const V2 &v2) {
1.148 + return v1 = t1 * v1 + t2 * prod (m, v2);
1.149 + }
1.150 +
1.151 + /** \brief rank 1 update: \a m = \a m + \a t * (\a v1 * \a v2<sup>T</sup>)
1.152 + \ingroup blas2
1.153 + */
1.154 + template<class M, class T, class V1, class V2>
1.155 + M &
1.156 + gr (M &m, const T &t, const V1 &v1, const V2 &v2) {
1.157 +#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1.158 + return m += t * outer_prod (v1, v2);
1.159 +#else
1.160 + return m = m + t * outer_prod (v1, v2);
1.161 +#endif
1.162 + }
1.163 +
1.164 + /** \brief symmetric rank 1 update: \a m = \a m + \a t * (\a v * \a v<sup>T</sup>)
1.165 + \ingroup blas2
1.166 + */
1.167 + template<class M, class T, class V>
1.168 + M &
1.169 + sr (M &m, const T &t, const V &v) {
1.170 +#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1.171 + return m += t * outer_prod (v, v);
1.172 +#else
1.173 + return m = m + t * outer_prod (v, v);
1.174 +#endif
1.175 + }
1.176 + /** \brief hermitian rank 1 update: \a m = \a m + \a t * (\a v * \a v<sup>H</sup>)
1.177 + \ingroup blas2
1.178 + */
1.179 + template<class M, class T, class V>
1.180 + M &
1.181 + hr (M &m, const T &t, const V &v) {
1.182 +#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1.183 + return m += t * outer_prod (v, conj (v));
1.184 +#else
1.185 + return m = m + t * outer_prod (v, conj (v));
1.186 +#endif
1.187 + }
1.188 +
1.189 + /** \brief symmetric rank 2 update: \a m = \a m + \a t *
1.190 + (\a v1 * \a v2<sup>T</sup> + \a v2 * \a v1<sup>T</sup>)
1.191 + \ingroup blas2
1.192 + */
1.193 + template<class M, class T, class V1, class V2>
1.194 + M &
1.195 + sr2 (M &m, const T &t, const V1 &v1, const V2 &v2) {
1.196 +#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1.197 + return m += t * (outer_prod (v1, v2) + outer_prod (v2, v1));
1.198 +#else
1.199 + return m = m + t * (outer_prod (v1, v2) + outer_prod (v2, v1));
1.200 +#endif
1.201 + }
1.202 + /** \brief hermitian rank 2 update: \a m = \a m +
1.203 + \a t * (\a v1 * \a v2<sup>H</sup>)
1.204 + + \a v2 * (\a t * \a v1)<sup>H</sup>)
1.205 + \ingroup blas2
1.206 + */
1.207 + template<class M, class T, class V1, class V2>
1.208 + M &
1.209 + hr2 (M &m, const T &t, const V1 &v1, const V2 &v2) {
1.210 +#ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
1.211 + return m += t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
1.212 +#else
1.213 + return m = m + t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
1.214 +#endif
1.215 + }
1.216 +
1.217 + }
1.218 +
1.219 + namespace blas_3 {
1.220 +
1.221 + /** \namespace boost::numeric::ublas::blas_3
1.222 + \brief wrapper functions for level 3 blas
1.223 + */
1.224 +
1.225 + /** \brief triangular matrix multiplication
1.226 + \ingroup blas3
1.227 + */
1.228 + template<class M1, class T, class M2, class M3>
1.229 + M1 &
1.230 + tmm (M1 &m1, const T &t, const M2 &m2, const M3 &m3) {
1.231 + return m1 = t * prod (m2, m3);
1.232 + }
1.233 +
1.234 + /** \brief triangular solve \a m2 * \a x = \a t * \a m1 in place,
1.235 + \a m2 is a triangular matrix
1.236 + \ingroup blas3
1.237 + */
1.238 + template<class M1, class T, class M2, class C>
1.239 + M1 &
1.240 + tsm (M1 &m1, const T &t, const M2 &m2, C) {
1.241 + return m1 = solve (m2, t * m1, C ());
1.242 + }
1.243 +
1.244 + /** \brief general matrix multiplication
1.245 + \ingroup blas3
1.246 + */
1.247 + template<class M1, class T1, class T2, class M2, class M3>
1.248 + M1 &
1.249 + gmm (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
1.250 + return m1 = t1 * m1 + t2 * prod (m2, m3);
1.251 + }
1.252 +
1.253 + /** \brief symmetric rank k update: \a m1 = \a t * \a m1 +
1.254 + \a t2 * (\a m2 * \a m2<sup>T</sup>)
1.255 + \ingroup blas3
1.256 + \todo use opb_prod()
1.257 + */
1.258 + template<class M1, class T1, class T2, class M2>
1.259 + M1 &
1.260 + srk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) {
1.261 + return m1 = t1 * m1 + t2 * prod (m2, trans (m2));
1.262 + }
1.263 + /** \brief hermitian rank k update: \a m1 = \a t * \a m1 +
1.264 + \a t2 * (\a m2 * \a m2<sup>H</sup>)
1.265 + \ingroup blas3
1.266 + \todo use opb_prod()
1.267 + */
1.268 + template<class M1, class T1, class T2, class M2>
1.269 + M1 &
1.270 + hrk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) {
1.271 + return m1 = t1 * m1 + t2 * prod (m2, herm (m2));
1.272 + }
1.273 +
1.274 + /** \brief generalized symmetric rank k update:
1.275 + \a m1 = \a t1 * \a m1 + \a t2 * (\a m2 * \a m3<sup>T</sup>)
1.276 + + \a t2 * (\a m3 * \a m2<sup>T</sup>)
1.277 + \ingroup blas3
1.278 + \todo use opb_prod()
1.279 + */
1.280 + template<class M1, class T1, class T2, class M2, class M3>
1.281 + M1 &
1.282 + sr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
1.283 + return m1 = t1 * m1 + t2 * (prod (m2, trans (m3)) + prod (m3, trans (m2)));
1.284 + }
1.285 + /** \brief generalized hermitian rank k update:
1.286 + \a m1 = \a t1 * \a m1 + \a t2 * (\a m2 * \a m3<sup>H</sup>)
1.287 + + (\a m3 * (\a t2 * \a m2)<sup>H</sup>)
1.288 + \ingroup blas3
1.289 + \todo use opb_prod()
1.290 + */
1.291 + template<class M1, class T1, class T2, class M2, class M3>
1.292 + M1 &
1.293 + hr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
1.294 + return m1 = t1 * m1 + t2 * prod (m2, herm (m3)) + type_traits<T2>::conj (t2) * prod (m3, herm (m2));
1.295 + }
1.296 +
1.297 + }
1.298 +
1.299 +}}}
1.300 +
1.301 +#endif
1.302 +
1.303 +