sl@0: //
sl@0: //  Copyright (c) 2000-2002
sl@0: //  Joerg Walter, Mathias Koch
sl@0: //
sl@0: //  Permission to use, copy, modify, distribute and sell this software
sl@0: //  and its documentation for any purpose is hereby granted without fee,
sl@0: //  provided that the above copyright notice appear in all copies and
sl@0: //  that both that copyright notice and this permission notice appear
sl@0: //  in supporting documentation.  The authors make no representations
sl@0: //  about the suitability of this software for any purpose.
sl@0: //  It is provided "as is" without express or implied warranty.
sl@0: //
sl@0: //  The authors gratefully acknowledge the support of
sl@0: //  GeNeSys mbH & Co. KG in producing this work.
sl@0: //
sl@0: 
sl@0: #ifndef _BOOST_UBLAS_BLAS_
sl@0: #define _BOOST_UBLAS_BLAS_
sl@0: 
sl@0: #include <boost/numeric/ublas/traits.hpp>
sl@0: 
sl@0: namespace boost { namespace numeric { namespace ublas {
sl@0: 
sl@0:     namespace blas_1 {
sl@0: 
sl@0:           /** \namespace boost::numeric::ublas::blas_1
sl@0:                   \brief wrapper functions for level 1 blas
sl@0:           */
sl@0: 
sl@0: 
sl@0:           /** \brief 1-Norm: \f$\sum_i |x_i|\f$
sl@0:                   \ingroup blas1
sl@0:            */
sl@0:         template<class V>
sl@0:         typename type_traits<typename V::value_type>::real_type
sl@0:         asum (const V &v) {
sl@0:             return norm_1 (v);
sl@0:         }
sl@0:           /** \brief 2-Norm: \f$\sum_i |x_i|^2\f$
sl@0:                   \ingroup blas1
sl@0:            */
sl@0:         template<class V>
sl@0:         typename type_traits<typename V::value_type>::real_type
sl@0:         nrm2 (const V &v) {
sl@0:             return norm_2 (v);
sl@0:         }
sl@0:           /** \brief element with larges absolute value: \f$\max_i |x_i|\f$
sl@0:                   \ingroup blas1
sl@0:           */                 
sl@0:         template<class V>
sl@0:         typename type_traits<typename V::value_type>::real_type
sl@0:         amax (const V &v) {
sl@0:             return norm_inf (v);
sl@0:         }
sl@0: 
sl@0:           /** \brief inner product of vectors \a v1 and \a v2
sl@0:                   \ingroup blas1
sl@0:           */                 
sl@0:         template<class V1, class V2>
sl@0:         typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type
sl@0:         dot (const V1 &v1, const V2 &v2) {
sl@0:             return inner_prod (v1, v2);
sl@0:         }
sl@0: 
sl@0:           /** \brief copy vector \a v2 to \a v1
sl@0:                   \ingroup blas1
sl@0:           */                 
sl@0:         template<class V1, class V2>
sl@0:         V1 &
sl@0:         copy (V1 &v1, const V2 &v2) {
sl@0:             return v1.assign (v2);
sl@0:         }
sl@0: 
sl@0:           /** \brief swap vectors \a v1 and \a v2
sl@0:                   \ingroup blas1
sl@0:           */                 
sl@0:         template<class V1, class V2>
sl@0:         void swap (V1 &v1, V2 &v2) {
sl@0:             v1.swap (v2);
sl@0:         }
sl@0: 
sl@0:           /** \brief scale vector \a v with scalar \a t
sl@0:                   \ingroup blas1
sl@0:           */                 
sl@0:         template<class V, class T>
sl@0:         V &
sl@0:         scal (V &v, const T &t) {
sl@0:             return v *= t;
sl@0:         }
sl@0: 
sl@0:           /** \brief compute \a v1 = \a v1 + \a t * \a v2
sl@0:                   \ingroup blas1
sl@0:           */                 
sl@0:         template<class V1, class T, class V2>
sl@0:         V1 &
sl@0:         axpy (V1 &v1, const T &t, const V2 &v2) {
sl@0:             return v1.plus_assign (t * v2);
sl@0:         }
sl@0: 
sl@0:           /** \brief apply plane rotation
sl@0:                   \ingroup blas1
sl@0:           */                 
sl@0:         template<class T1, class V1, class T2, class V2>
sl@0:         void
sl@0:         rot (const T1 &t1, V1 &v1, const T2 &t2, V2 &v2) {
sl@0:             typedef typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type promote_type;
sl@0:             vector<promote_type> vt (t1 * v1 + t2 * v2);
sl@0:             v2.assign (- t2 * v1 + t1 * v2);
sl@0:             v1.assign (vt);
sl@0:         }
sl@0: 
sl@0:     }
sl@0: 
sl@0:     namespace blas_2 {
sl@0: 
sl@0:           /** \namespace boost::numeric::ublas::blas_2
sl@0:                   \brief wrapper functions for level 2 blas
sl@0:           */
sl@0: 
sl@0:           /** \brief multiply vector \a v with triangular matrix \a m
sl@0:                   \ingroup blas2
sl@0:                   \todo: check that matrix is really triangular
sl@0:           */                 
sl@0:         template<class V, class M>
sl@0:         V &
sl@0:         tmv (V &v, const M &m) {
sl@0:             return v = prod (m, v);
sl@0:         }
sl@0: 
sl@0:           /** \brief solve \a m \a x = \a v in place, \a m is triangular matrix
sl@0:                   \ingroup blas2
sl@0:           */                 
sl@0:         template<class V, class M, class C>
sl@0:         V &
sl@0:         tsv (V &v, const M &m, C) {
sl@0:             return v = solve (m, v, C ());
sl@0:         }
sl@0: 
sl@0:           /** \brief compute \a v1 = \a t1 * \a v1 + \a t2 * (\a m * \a v2)
sl@0:                   \ingroup blas2
sl@0:           */                 
sl@0:         template<class V1, class T1, class T2, class M, class V2>
sl@0:         V1 &
sl@0:         gmv (V1 &v1, const T1 &t1, const T2 &t2, const M &m, const V2 &v2) {
sl@0:             return v1 = t1 * v1 + t2 * prod (m, v2);
sl@0:         }
sl@0: 
sl@0:           /** \brief rank 1 update: \a m = \a m + \a t * (\a v1 * \a v2<sup>T</sup>)
sl@0:                   \ingroup blas2
sl@0:           */                 
sl@0:         template<class M, class T, class V1, class V2>
sl@0:         M &
sl@0:         gr (M &m, const T &t, const V1 &v1, const V2 &v2) {
sl@0: #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
sl@0:             return m += t * outer_prod (v1, v2);
sl@0: #else
sl@0:             return m = m + t * outer_prod (v1, v2);
sl@0: #endif
sl@0:         }
sl@0: 
sl@0:           /** \brief symmetric rank 1 update: \a m = \a m + \a t * (\a v * \a v<sup>T</sup>)
sl@0:                   \ingroup blas2
sl@0:           */                 
sl@0:         template<class M, class T, class V>
sl@0:         M &
sl@0:         sr (M &m, const T &t, const V &v) {
sl@0: #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
sl@0:             return m += t * outer_prod (v, v);
sl@0: #else
sl@0:             return m = m + t * outer_prod (v, v);
sl@0: #endif
sl@0:         }
sl@0:           /** \brief hermitian rank 1 update: \a m = \a m + \a t * (\a v * \a v<sup>H</sup>)
sl@0:                   \ingroup blas2
sl@0:           */                 
sl@0:         template<class M, class T, class V>
sl@0:         M &
sl@0:         hr (M &m, const T &t, const V &v) {
sl@0: #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
sl@0:             return m += t * outer_prod (v, conj (v));
sl@0: #else
sl@0:             return m = m + t * outer_prod (v, conj (v));
sl@0: #endif
sl@0:         }
sl@0: 
sl@0:           /** \brief symmetric rank 2 update: \a m = \a m + \a t * 
sl@0:                   (\a v1 * \a v2<sup>T</sup> + \a v2 * \a v1<sup>T</sup>) 
sl@0:                   \ingroup blas2
sl@0:           */                 
sl@0:         template<class M, class T, class V1, class V2>
sl@0:         M &
sl@0:         sr2 (M &m, const T &t, const V1 &v1, const V2 &v2) {
sl@0: #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
sl@0:             return m += t * (outer_prod (v1, v2) + outer_prod (v2, v1));
sl@0: #else
sl@0:             return m = m + t * (outer_prod (v1, v2) + outer_prod (v2, v1));
sl@0: #endif
sl@0:         }
sl@0:           /** \brief hermitian rank 2 update: \a m = \a m + 
sl@0:                   \a t * (\a v1 * \a v2<sup>H</sup>)
sl@0:                   + \a v2 * (\a t * \a v1)<sup>H</sup>) 
sl@0:                   \ingroup blas2
sl@0:           */                 
sl@0:         template<class M, class T, class V1, class V2>
sl@0:         M &
sl@0:         hr2 (M &m, const T &t, const V1 &v1, const V2 &v2) {
sl@0: #ifndef BOOST_UBLAS_SIMPLE_ET_DEBUG
sl@0:             return m += t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
sl@0: #else
sl@0:             return m = m + t * outer_prod (v1, conj (v2)) + type_traits<T>::conj (t) * outer_prod (v2, conj (v1));
sl@0: #endif
sl@0:         }
sl@0: 
sl@0:     }
sl@0: 
sl@0:     namespace blas_3 {
sl@0: 
sl@0:           /** \namespace boost::numeric::ublas::blas_3
sl@0:                   \brief wrapper functions for level 3 blas
sl@0:           */
sl@0: 
sl@0:           /** \brief triangular matrix multiplication
sl@0:                   \ingroup blas3
sl@0:           */                 
sl@0:         template<class M1, class T, class M2, class M3>
sl@0:         M1 &
sl@0:         tmm (M1 &m1, const T &t, const M2 &m2, const M3 &m3) {
sl@0:             return m1 = t * prod (m2, m3);
sl@0:         }
sl@0: 
sl@0:           /** \brief triangular solve \a m2 * \a x = \a t * \a m1 in place,
sl@0:                   \a m2 is a triangular matrix
sl@0:                   \ingroup blas3
sl@0:           */                 
sl@0:         template<class M1, class T, class M2, class C>
sl@0:         M1 &
sl@0:         tsm (M1 &m1, const T &t, const M2 &m2, C) {
sl@0:             return m1 = solve (m2, t * m1, C ());
sl@0:         }
sl@0: 
sl@0:           /** \brief general matrix multiplication
sl@0:                   \ingroup blas3
sl@0:           */                 
sl@0:         template<class M1, class T1, class T2, class M2, class M3>
sl@0:         M1 &
sl@0:         gmm (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
sl@0:             return m1 = t1 * m1 + t2 * prod (m2, m3);
sl@0:         }
sl@0: 
sl@0:           /** \brief symmetric rank k update: \a m1 = \a t * \a m1 + 
sl@0:                   \a t2 * (\a m2 * \a m2<sup>T</sup>)
sl@0:                   \ingroup blas3
sl@0:                   \todo use opb_prod()
sl@0:           */                 
sl@0:         template<class M1, class T1, class T2, class M2>
sl@0:         M1 &
sl@0:         srk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) {
sl@0:             return m1 = t1 * m1 + t2 * prod (m2, trans (m2));
sl@0:         }
sl@0:           /** \brief hermitian rank k update: \a m1 = \a t * \a m1 + 
sl@0:                   \a t2 * (\a m2 * \a m2<sup>H</sup>)
sl@0:                   \ingroup blas3
sl@0:                   \todo use opb_prod()
sl@0:           */                 
sl@0:         template<class M1, class T1, class T2, class M2>
sl@0:         M1 &
sl@0:         hrk (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2) {
sl@0:             return m1 = t1 * m1 + t2 * prod (m2, herm (m2));
sl@0:         }
sl@0: 
sl@0:           /** \brief generalized symmetric rank k update:
sl@0:                   \a m1 = \a t1 * \a m1 + \a t2 * (\a m2 * \a m3<sup>T</sup>)
sl@0:                   + \a t2 * (\a m3 * \a m2<sup>T</sup>)
sl@0:                   \ingroup blas3
sl@0:                   \todo use opb_prod()
sl@0:           */                 
sl@0:         template<class M1, class T1, class T2, class M2, class M3>
sl@0:         M1 &
sl@0:         sr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
sl@0:             return m1 = t1 * m1 + t2 * (prod (m2, trans (m3)) + prod (m3, trans (m2)));
sl@0:         }
sl@0:           /** \brief generalized hermitian rank k update:
sl@0:                   \a m1 = \a t1 * \a m1 + \a t2 * (\a m2 * \a m3<sup>H</sup>)
sl@0:                   + (\a m3 * (\a t2 * \a m2)<sup>H</sup>)
sl@0:                   \ingroup blas3
sl@0:                   \todo use opb_prod()
sl@0:           */                 
sl@0:         template<class M1, class T1, class T2, class M2, class M3>
sl@0:         M1 &
sl@0:         hr2k (M1 &m1, const T1 &t1, const T2 &t2, const M2 &m2, const M3 &m3) {
sl@0:             return m1 = t1 * m1 + t2 * prod (m2, herm (m3)) + type_traits<T2>::conj (t2) * prod (m3, herm (m2));
sl@0:         }
sl@0: 
sl@0:     }
sl@0: 
sl@0: }}}
sl@0: 
sl@0: #endif
sl@0: 
sl@0: