Symaptic: os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/operation.hpp@260cb5ec6c19

     1 //

     2 //  Copyright (c) 2000-2002

     3 //  Joerg Walter, Mathias Koch

     4 //

     5 //  Permission to use, copy, modify, distribute and sell this software

     6 //  and its documentation for any purpose is hereby granted without fee,

     7 //  provided that the above copyright notice appear in all copies and

     8 //  that both that copyright notice and this permission notice appear

     9 //  in supporting documentation.  The authors make no representations

    10 //  about the suitability of this software for any purpose.

    11 //  It is provided "as is" without express or implied warranty.

    12 //

    13 //  The authors gratefully acknowledge the support of

    14 //  GeNeSys mbH & Co. KG in producing this work.

    15 //

    17 #ifndef _BOOST_UBLAS_OPERATION_

    18 #define _BOOST_UBLAS_OPERATION_

    20 #include <boost/numeric/ublas/matrix_proxy.hpp>

    22 /** \file operation.hpp

    23  *  \brief This file contains some specialized products.

    24  */

    26 // axpy-based products

    27 // Alexei Novakov had a lot of ideas to improve these. Thanks.

    28 // Hendrik Kueck proposed some new kernel. Thanks again.

    30 namespace boost { namespace numeric { namespace ublas {

    32     template<class V, class T1, class IA1, class TA1, class E2>

    33     BOOST_UBLAS_INLINE

    34     V &

    35     axpy_prod (const compressed_matrix<T1, row_major, 0, IA1, TA1> &e1,

    36                const vector_expression<E2> &e2,

    37                V &v, row_major_tag) {

    38         typedef typename V::size_type size_type;

    39         typedef typename V::value_type value_type;

    41         for (size_type i = 0; i < e1.filled1 () -1; ++ i) {

    42             size_type begin = e1.index1_data () [i];

    43             size_type end = e1.index1_data () [i + 1];

    44             value_type t (v (i));

    45             for (size_type j = begin; j < end; ++ j)

    46                 t += e1.value_data () [j] * e2 () (e1.index2_data () [j]);

    47             v (i) = t;

    48         }

    49         return v;

    50     }

    52     template<class V, class T1, class IA1, class TA1, class E2>

    53     BOOST_UBLAS_INLINE

    54     V &

    55     axpy_prod (const compressed_matrix<T1, column_major, 0, IA1, TA1> &e1,

    56                const vector_expression<E2> &e2,

    57                V &v, column_major_tag) {

    58         typedef typename V::size_type size_type;

    60         for (size_type j = 0; j < e1.filled1 () -1; ++ j) {

    61             size_type begin = e1.index1_data () [j];

    62             size_type end = e1.index1_data () [j + 1];

    63             for (size_type i = begin; i < end; ++ i)

    64                 v (e1.index2_data () [i]) += e1.value_data () [i] * e2 () (j);

    65         }

    66         return v;

    67     }

    69     // Dispatcher

    70     template<class V, class T1, class L1, class IA1, class TA1, class E2>

    71     BOOST_UBLAS_INLINE

    72     V &

    73     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,

    74                const vector_expression<E2> &e2,

    75                V &v, bool init = true) {

    76         typedef typename V::value_type value_type;

    77         typedef typename L1::orientation_category orientation_category;

    79         if (init)

    80             v.assign (zero_vector<value_type> (e1.size1 ()));

    81 #if BOOST_UBLAS_TYPE_CHECK

    82         vector<value_type> cv (v);

    83         typedef typename type_traits<value_type>::real_type real_type;

    84         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));

    85         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));

    86 #endif

    87         axpy_prod (e1, e2, v, orientation_category ());

    88 #if BOOST_UBLAS_TYPE_CHECK

    89         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());

    90 #endif

    91         return v;

    92     }

    93     template<class V, class T1, class L1, class IA1, class TA1, class E2>

    94     BOOST_UBLAS_INLINE

    95     V

    96     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,

    97                const vector_expression<E2> &e2) {

    98         typedef V vector_type;

   100         vector_type v (e1.size1 ());

   101         return axpy_prod (e1, e2, v, true);

   102     }

   104     template<class V, class T1, class L1, class IA1, class TA1, class E2>

   105     BOOST_UBLAS_INLINE

   106     V &

   107     axpy_prod (const coordinate_matrix<T1, L1, 0, IA1, TA1> &e1,

   108                const vector_expression<E2> &e2,

   109                V &v, bool init = true) {

   110         typedef typename V::size_type size_type;

   111         typedef typename V::value_type value_type;

   112         typedef L1 layout_type;

   114         size_type size1 = e1.size1();

   115         size_type size2 = e1.size2();

   117         if (init) {

   118             noalias(v) = zero_vector<value_type>(size1);

   119         }

   121         for (size_type i = 0; i < e1.nnz(); ++i) {

   122             size_type row_index = layout_type::element1( e1.index1_data () [i], size1, e1.index2_data () [i], size2 );

   123             size_type col_index = layout_type::element2( e1.index1_data () [i], size1, e1.index2_data () [i], size2 );

   124             v( row_index ) += e1.value_data () [i] * e2 () (col_index);

   125         }

   126         return v;

   127     }

   129     template<class V, class E1, class E2>

   130     BOOST_UBLAS_INLINE

   131     V &

   132     axpy_prod (const matrix_expression<E1> &e1,

   133                const vector_expression<E2> &e2,

   134                V &v, packed_random_access_iterator_tag, row_major_tag) {

   135         typedef const E1 expression1_type;

   136         typedef const E2 expression2_type;

   137         typedef typename V::size_type size_type;

   139         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());

   140         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());

   141         while (it1 != it1_end) {

   142             size_type index1 (it1.index1 ());

   143 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   144             typename expression1_type::const_iterator2 it2 (it1.begin ());

   145             typename expression1_type::const_iterator2 it2_end (it1.end ());

   146 #else

   147             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));

   148             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));

   149 #endif

   150             while (it2 != it2_end) {

   151                 v (index1) += *it2 * e2 () (it2.index2 ());

   152                 ++ it2;

   153             }

   154             ++ it1;

   155         }

   156         return v;

   157     }

   159     template<class V, class E1, class E2>

   160     BOOST_UBLAS_INLINE

   161     V &

   162     axpy_prod (const matrix_expression<E1> &e1,

   163                const vector_expression<E2> &e2,

   164                V &v, packed_random_access_iterator_tag, column_major_tag) {

   165         typedef const E1 expression1_type;

   166         typedef const E2 expression2_type;

   167         typedef typename V::size_type size_type;

   169         typename expression1_type::const_iterator2 it2 (e1 ().begin2 ());

   170         typename expression1_type::const_iterator2 it2_end (e1 ().end2 ());

   171         while (it2 != it2_end) {

   172             size_type index2 (it2.index2 ());

   173 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   174             typename expression1_type::const_iterator1 it1 (it2.begin ());

   175             typename expression1_type::const_iterator1 it1_end (it2.end ());

   176 #else

   177             typename expression1_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));

   178             typename expression1_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));

   179 #endif

   180             while (it1 != it1_end) {

   181                 v (it1.index1 ()) += *it1 * e2 () (index2);

   182                 ++ it1;

   183             }

   184             ++ it2;

   185         }

   186         return v;

   187     }

   189     template<class V, class E1, class E2>

   190     BOOST_UBLAS_INLINE

   191     V &

   192     axpy_prod (const matrix_expression<E1> &e1,

   193                const vector_expression<E2> &e2,

   194                V &v, sparse_bidirectional_iterator_tag) {

   195         typedef const E1 expression1_type;

   196         typedef const E2 expression2_type;

   197         typedef typename V::size_type size_type;

   199         typename expression2_type::const_iterator it (e2 ().begin ());

   200         typename expression2_type::const_iterator it_end (e2 ().end ());

   201         while (it != it_end) {

   202             v.plus_assign (column (e1 (), it.index ()) * *it);

   203             ++ it;

   204         }

   205         return v;

   206     }

   208     // Dispatcher

   209     template<class V, class E1, class E2>

   210     BOOST_UBLAS_INLINE

   211     V &

   212     axpy_prod (const matrix_expression<E1> &e1,

   213                const vector_expression<E2> &e2,

   214                V &v, packed_random_access_iterator_tag) {

   215         typedef typename E1::orientation_category orientation_category;

   216         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());

   217     }

   220   /** \brief computes <tt>v += A x</tt> or <tt>v = A x</tt> in an

   221           optimized fashion.

   223           \param e1 the matrix expression \c A

   224           \param e2 the vector expression \c x

   225           \param v  the result vector \c v

   226           \param init a boolean parameter

   228           <tt>axpy_prod(A, x, v, init)</tt> implements the well known

   229           axpy-product.  Setting \a init to \c true is equivalent to call

   230           <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init

   231           defaults to \c true, but this may change in the future.

   233           Up to now there are some specialisation for compressed

   234           matrices that give a large speed up compared to prod.

   236           \ingroup blas2

   238           \internal

   240           template parameters:

   241           \param V type of the result vector \c v

   242           \param E1 type of a matrix expression \c A

   243           \param E2 type of a vector expression \c x

   244   */

   245     template<class V, class E1, class E2>

   246     BOOST_UBLAS_INLINE

   247     V &

   248     axpy_prod (const matrix_expression<E1> &e1,

   249                const vector_expression<E2> &e2,

   250                V &v, bool init = true) {

   251         typedef typename V::value_type value_type;

   252         typedef typename E2::const_iterator::iterator_category iterator_category;

   254         if (init)

   255             v.assign (zero_vector<value_type> (e1 ().size1 ()));

   256 #if BOOST_UBLAS_TYPE_CHECK

   257         vector<value_type> cv (v);

   258         typedef typename type_traits<value_type>::real_type real_type;

   259         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));

   260         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));

   261 #endif

   262         axpy_prod (e1, e2, v, iterator_category ());

   263 #if BOOST_UBLAS_TYPE_CHECK

   264         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());

   265 #endif

   266         return v;

   267     }

   268     template<class V, class E1, class E2>

   269     BOOST_UBLAS_INLINE

   270     V

   271     axpy_prod (const matrix_expression<E1> &e1,

   272                const vector_expression<E2> &e2) {

   273         typedef V vector_type;

   275         vector_type v (e1 ().size1 ());

   276         return axpy_prod (e1, e2, v, true);

   277     }

   279     template<class V, class E1, class T2, class IA2, class TA2>

   280     BOOST_UBLAS_INLINE

   281     V &

   282     axpy_prod (const vector_expression<E1> &e1,

   283                const compressed_matrix<T2, column_major, 0, IA2, TA2> &e2,

   284                V &v, column_major_tag) {

   285         typedef typename V::size_type size_type;

   286         typedef typename V::value_type value_type;

   288         for (size_type j = 0; j < e2.filled1 () -1; ++ j) {

   289             size_type begin = e2.index1_data () [j];

   290             size_type end = e2.index1_data () [j + 1];

   291             value_type t (v (j));

   292             for (size_type i = begin; i < end; ++ i)

   293                 t += e2.value_data () [i] * e1 () (e2.index2_data () [i]);

   294             v (j) = t;

   295         }

   296         return v;

   297     }

   299     template<class V, class E1, class T2, class IA2, class TA2>

   300     BOOST_UBLAS_INLINE

   301     V &

   302     axpy_prod (const vector_expression<E1> &e1,

   303                const compressed_matrix<T2, row_major, 0, IA2, TA2> &e2,

   304                V &v, row_major_tag) {

   305         typedef typename V::size_type size_type;

   307         for (size_type i = 0; i < e2.filled1 () -1; ++ i) {

   308             size_type begin = e2.index1_data () [i];

   309             size_type end = e2.index1_data () [i + 1];

   310             for (size_type j = begin; j < end; ++ j)

   311                 v (e2.index2_data () [j]) += e2.value_data () [j] * e1 () (i);

   312         }

   313         return v;

   314     }

   316     // Dispatcher

   317     template<class V, class E1, class T2, class L2, class IA2, class TA2>

   318     BOOST_UBLAS_INLINE

   319     V &

   320     axpy_prod (const vector_expression<E1> &e1,

   321                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2,

   322                V &v, bool init = true) {

   323         typedef typename V::value_type value_type;

   324         typedef typename L2::orientation_category orientation_category;

   326         if (init)

   327             v.assign (zero_vector<value_type> (e2.size2 ()));

   328 #if BOOST_UBLAS_TYPE_CHECK

   329         vector<value_type> cv (v);

   330         typedef typename type_traits<value_type>::real_type real_type;

   331         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));

   332         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));

   333 #endif

   334         axpy_prod (e1, e2, v, orientation_category ());

   335 #if BOOST_UBLAS_TYPE_CHECK

   336         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());

   337 #endif

   338         return v;

   339     }

   340     template<class V, class E1, class T2, class L2, class IA2, class TA2>

   341     BOOST_UBLAS_INLINE

   342     V

   343     axpy_prod (const vector_expression<E1> &e1,

   344                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2) {

   345         typedef V vector_type;

   347         vector_type v (e2.size2 ());

   348         return axpy_prod (e1, e2, v, true);

   349     }

   351     template<class V, class E1, class E2>

   352     BOOST_UBLAS_INLINE

   353     V &

   354     axpy_prod (const vector_expression<E1> &e1,

   355                const matrix_expression<E2> &e2,

   356                V &v, packed_random_access_iterator_tag, column_major_tag) {

   357         typedef const E1 expression1_type;

   358         typedef const E2 expression2_type;

   359         typedef typename V::size_type size_type;

   361         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());

   362         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());

   363         while (it2 != it2_end) {

   364             size_type index2 (it2.index2 ());

   365 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   366             typename expression2_type::const_iterator1 it1 (it2.begin ());

   367             typename expression2_type::const_iterator1 it1_end (it2.end ());

   368 #else

   369             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));

   370             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));

   371 #endif

   372             while (it1 != it1_end) {

   373                 v (index2) += *it1 * e1 () (it1.index1 ());

   374                 ++ it1;

   375             }

   376             ++ it2;

   377         }

   378         return v;

   379     }

   381     template<class V, class E1, class E2>

   382     BOOST_UBLAS_INLINE

   383     V &

   384     axpy_prod (const vector_expression<E1> &e1,

   385                const matrix_expression<E2> &e2,

   386                V &v, packed_random_access_iterator_tag, row_major_tag) {

   387         typedef const E1 expression1_type;

   388         typedef const E2 expression2_type;

   389         typedef typename V::size_type size_type;

   391         typename expression2_type::const_iterator1 it1 (e2 ().begin1 ());

   392         typename expression2_type::const_iterator1 it1_end (e2 ().end1 ());

   393         while (it1 != it1_end) {

   394             size_type index1 (it1.index1 ());

   395 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   396             typename expression2_type::const_iterator2 it2 (it1.begin ());

   397             typename expression2_type::const_iterator2 it2_end (it1.end ());

   398 #else

   399             typename expression2_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));

   400             typename expression2_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));

   401 #endif

   402             while (it2 != it2_end) {

   403                 v (it2.index2 ()) += *it2 * e1 () (index1);

   404                 ++ it2;

   405             }

   406             ++ it1;

   407         }

   408         return v;

   409     }

   411     template<class V, class E1, class E2>

   412     BOOST_UBLAS_INLINE

   413     V &

   414     axpy_prod (const vector_expression<E1> &e1,

   415                const matrix_expression<E2> &e2,

   416                V &v, sparse_bidirectional_iterator_tag) {

   417         typedef const E1 expression1_type;

   418         typedef const E2 expression2_type;

   419         typedef typename V::size_type size_type;

   421         typename expression1_type::const_iterator it (e1 ().begin ());

   422         typename expression1_type::const_iterator it_end (e1 ().end ());

   423         while (it != it_end) {

   424             v.plus_assign (*it * row (e2 (), it.index ()));

   425             ++ it;

   426         }

   427         return v;

   428     }

   430     // Dispatcher

   431     template<class V, class E1, class E2>

   432     BOOST_UBLAS_INLINE

   433     V &

   434     axpy_prod (const vector_expression<E1> &e1,

   435                const matrix_expression<E2> &e2,

   436                V &v, packed_random_access_iterator_tag) {

   437         typedef typename E2::orientation_category orientation_category;

   438         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());

   439     }

   442   /** \brief computes <tt>v += A<sup>T</sup> x</tt> or <tt>v = A<sup>T</sup> x</tt> in an

   443           optimized fashion.

   445           \param e1 the vector expression \c x

   446           \param e2 the matrix expression \c A

   447           \param v  the result vector \c v

   448           \param init a boolean parameter

   450           <tt>axpy_prod(x, A, v, init)</tt> implements the well known

   451           axpy-product.  Setting \a init to \c true is equivalent to call

   452           <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init

   453           defaults to \c true, but this may change in the future.

   455           Up to now there are some specialisation for compressed

   456           matrices that give a large speed up compared to prod.

   458           \ingroup blas2

   460           \internal

   462           template parameters:

   463           \param V type of the result vector \c v

   464           \param E1 type of a vector expression \c x

   465           \param E2 type of a matrix expression \c A

   466   */

   467     template<class V, class E1, class E2>

   468     BOOST_UBLAS_INLINE

   469     V &

   470     axpy_prod (const vector_expression<E1> &e1,

   471                const matrix_expression<E2> &e2,

   472                V &v, bool init = true) {

   473         typedef typename V::value_type value_type;

   474         typedef typename E1::const_iterator::iterator_category iterator_category;

   476         if (init)

   477             v.assign (zero_vector<value_type> (e2 ().size2 ()));

   478 #if BOOST_UBLAS_TYPE_CHECK

   479         vector<value_type> cv (v);

   480         typedef typename type_traits<value_type>::real_type real_type;

   481         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));

   482         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));

   483 #endif

   484         axpy_prod (e1, e2, v, iterator_category ());

   485 #if BOOST_UBLAS_TYPE_CHECK

   486         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());

   487 #endif

   488         return v;

   489     }

   490     template<class V, class E1, class E2>

   491     BOOST_UBLAS_INLINE

   492     V

   493     axpy_prod (const vector_expression<E1> &e1,

   494                const matrix_expression<E2> &e2) {

   495         typedef V vector_type;

   497         vector_type v (e2 ().size2 ());

   498         return axpy_prod (e1, e2, v, true);

   499     }

   501     template<class M, class E1, class E2, class TRI>

   502     BOOST_UBLAS_INLINE

   503     M &

   504     axpy_prod (const matrix_expression<E1> &e1,

   505                const matrix_expression<E2> &e2,

   506                M &m, TRI,

   507                dense_proxy_tag, row_major_tag) {

   508         typedef M matrix_type;

   509         typedef const E1 expression1_type;

   510         typedef const E2 expression2_type;

   511         typedef typename M::size_type size_type;

   512         typedef typename M::value_type value_type;

   514 #if BOOST_UBLAS_TYPE_CHECK

   515         matrix<value_type, row_major> cm (m);

   516         typedef typename type_traits<value_type>::real_type real_type;

   517         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));

   518         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());

   519 #endif

   520         size_type size1 (e1 ().size1 ());

   521         size_type size2 (e1 ().size2 ());

   522         for (size_type i = 0; i < size1; ++ i)

   523             for (size_type j = 0; j < size2; ++ j)

   524                 row (m, i).plus_assign (e1 () (i, j) * row (e2 (), j));

   525 #if BOOST_UBLAS_TYPE_CHECK

   526         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());

   527 #endif

   528         return m;

   529     }

   530     template<class M, class E1, class E2, class TRI>

   531     BOOST_UBLAS_INLINE

   532     M &

   533     axpy_prod (const matrix_expression<E1> &e1,

   534                const matrix_expression<E2> &e2,

   535                M &m, TRI,

   536                sparse_proxy_tag, row_major_tag) {

   537         typedef M matrix_type;

   538         typedef TRI triangular_restriction;

   539         typedef const E1 expression1_type;

   540         typedef const E2 expression2_type;

   541         typedef typename M::size_type size_type;

   542         typedef typename M::value_type value_type;

   544 #if BOOST_UBLAS_TYPE_CHECK

   545         matrix<value_type, row_major> cm (m);

   546         typedef typename type_traits<value_type>::real_type real_type;

   547         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));

   548         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());

   549 #endif

   550         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());

   551         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());

   552         while (it1 != it1_end) {

   553 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   554             typename expression1_type::const_iterator2 it2 (it1.begin ());

   555             typename expression1_type::const_iterator2 it2_end (it1.end ());

   556 #else

   557             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));

   558             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));

   559 #endif

   560             while (it2 != it2_end) {

   561                 // row (m, it1.index1 ()).plus_assign (*it2 * row (e2 (), it2.index2 ()));

   562                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());

   563                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());

   564                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());

   565                 while (itr != itr_end) {

   566                     if (triangular_restriction::other (it1.index1 (), itr.index ()))

   567                         m (it1.index1 (), itr.index ()) += *it2 * *itr;

   568                     ++ itr;

   569                 }

   570                 ++ it2;

   571             }

   572             ++ it1;

   573         }

   574 #if BOOST_UBLAS_TYPE_CHECK

   575         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());

   576 #endif

   577         return m;

   578     }

   580     template<class M, class E1, class E2, class TRI>

   581     BOOST_UBLAS_INLINE

   582     M &

   583     axpy_prod (const matrix_expression<E1> &e1,

   584                const matrix_expression<E2> &e2,

   585                M &m, TRI,

   586                dense_proxy_tag, column_major_tag) {

   587         typedef M matrix_type;

   588         typedef const E1 expression1_type;

   589         typedef const E2 expression2_type;

   590         typedef typename M::size_type size_type;

   591         typedef typename M::value_type value_type;

   593 #if BOOST_UBLAS_TYPE_CHECK

   594         matrix<value_type, column_major> cm (m);

   595         typedef typename type_traits<value_type>::real_type real_type;

   596         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));

   597         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());

   598 #endif

   599         size_type size1 (e2 ().size1 ());

   600         size_type size2 (e2 ().size2 ());

   601         for (size_type j = 0; j < size2; ++ j)

   602             for (size_type i = 0; i < size1; ++ i)

   603                 column (m, j).plus_assign (e2 () (i, j) * column (e1 (), i));

   604 #if BOOST_UBLAS_TYPE_CHECK

   605         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());

   606 #endif

   607         return m;

   608     }

   609     template<class M, class E1, class E2, class TRI>

   610     BOOST_UBLAS_INLINE

   611     M &

   612     axpy_prod (const matrix_expression<E1> &e1,

   613                const matrix_expression<E2> &e2,

   614                M &m, TRI,

   615                sparse_proxy_tag, column_major_tag) {

   616         typedef M matrix_type;

   617         typedef TRI triangular_restriction;

   618         typedef const E1 expression1_type;

   619         typedef const E2 expression2_type;

   620         typedef typename M::size_type size_type;

   621         typedef typename M::value_type value_type;

   623 #if BOOST_UBLAS_TYPE_CHECK

   624         matrix<value_type, column_major> cm (m);

   625         typedef typename type_traits<value_type>::real_type real_type;

   626         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));

   627         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());

   628 #endif

   629         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());

   630         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());

   631         while (it2 != it2_end) {

   632 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION

   633             typename expression2_type::const_iterator1 it1 (it2.begin ());

   634             typename expression2_type::const_iterator1 it1_end (it2.end ());

   635 #else

   636             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));

   637             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));

   638 #endif

   639             while (it1 != it1_end) {

   640                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));

   641                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());

   642                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());

   643                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());

   644                 while (itc != itc_end) {

   645                     if (triangular_restriction::functor_type ().other (itc.index (), it2.index2 ()))

   646                         m (itc.index (), it2.index2 ()) += *it1 * *itc;

   647                     ++ itc;

   648                 }

   649                 ++ it1;

   650             }

   651             ++ it2;

   652         }

   653 #if BOOST_UBLAS_TYPE_CHECK

   654         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());

   655 #endif

   656         return m;

   657     }

   659     // Dispatcher

   660     template<class M, class E1, class E2, class TRI>

   661     BOOST_UBLAS_INLINE

   662     M &

   663     axpy_prod (const matrix_expression<E1> &e1,

   664                const matrix_expression<E2> &e2,

   665                M &m, TRI, bool init = true) {

   666         typedef typename M::value_type value_type;

   667         typedef typename M::storage_category storage_category;

   668         typedef typename M::orientation_category orientation_category;

   669         typedef TRI triangular_restriction;

   671         if (init)

   672             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));

   673         return axpy_prod (e1, e2, m, triangular_restriction (), storage_category (), orientation_category ());

   674     }

   675     template<class M, class E1, class E2, class TRI>

   676     BOOST_UBLAS_INLINE

   677     M

   678     axpy_prod (const matrix_expression<E1> &e1,

   679                const matrix_expression<E2> &e2,

   680                TRI) {

   681         typedef M matrix_type;

   682         typedef TRI triangular_restriction;

   684         matrix_type m (e1 ().size1 (), e2 ().size2 ());

   685         return axpy_prod (e1, e2, m, triangular_restriction (), true);

   686     }

   688   /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an

   689           optimized fashion.

   691           \param e1 the matrix expression \c A

   692           \param e2 the matrix expression \c X

   693           \param m  the result matrix \c M

   694           \param init a boolean parameter

   696           <tt>axpy_prod(A, X, M, init)</tt> implements the well known

   697           axpy-product.  Setting \a init to \c true is equivalent to call

   698           <tt>M.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init

   699           defaults to \c true, but this may change in the future.

   701           Up to now there are no specialisations.

   703           \ingroup blas3

   705           \internal

   707           template parameters:

   708           \param M type of the result matrix \c M

   709           \param E1 type of a matrix expression \c A

   710           \param E2 type of a matrix expression \c X

   711   */

   712     template<class M, class E1, class E2>

   713     BOOST_UBLAS_INLINE

   714     M &

   715     axpy_prod (const matrix_expression<E1> &e1,

   716                const matrix_expression<E2> &e2,

   717                M &m, bool init = true) {

   718         typedef typename M::value_type value_type;

   719         typedef typename M::storage_category storage_category;

   720         typedef typename M::orientation_category orientation_category;

   722         if (init)

   723             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));

   724         return axpy_prod (e1, e2, m, full (), storage_category (), orientation_category ());

   725     }

   726     template<class M, class E1, class E2>

   727     BOOST_UBLAS_INLINE

   728     M

   729     axpy_prod (const matrix_expression<E1> &e1,

   730                const matrix_expression<E2> &e2) {

   731         typedef M matrix_type;

   733         matrix_type m (e1 ().size1 (), e2 ().size2 ());

   734         return axpy_prod (e1, e2, m, full (), true);

   735     }

   738     template<class M, class E1, class E2>

   739     BOOST_UBLAS_INLINE

   740     M &

   741     opb_prod (const matrix_expression<E1> &e1,

   742               const matrix_expression<E2> &e2,

   743               M &m,

   744               dense_proxy_tag, row_major_tag) {

   745         typedef M matrix_type;

   746         typedef const E1 expression1_type;

   747         typedef const E2 expression2_type;

   748         typedef typename M::size_type size_type;

   749         typedef typename M::value_type value_type;

   751 #if BOOST_UBLAS_TYPE_CHECK

   752         matrix<value_type, row_major> cm (m);

   753         typedef typename type_traits<value_type>::real_type real_type;

   754         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));

   755         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());

   756 #endif

   757         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));

   758         for (size_type k = 0; k < size; ++ k) {

   759             vector<value_type> ce1 (column (e1 (), k));

   760             vector<value_type> re2 (row (e2 (), k));

   761             m.plus_assign (outer_prod (ce1, re2));

   762         }

   763 #if BOOST_UBLAS_TYPE_CHECK

   764         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());

   765 #endif

   766         return m;

   767     }

   769     template<class M, class E1, class E2>

   770     BOOST_UBLAS_INLINE

   771     M &

   772     opb_prod (const matrix_expression<E1> &e1,

   773               const matrix_expression<E2> &e2,

   774               M &m,

   775               dense_proxy_tag, column_major_tag) {

   776         typedef M matrix_type;

   777         typedef const E1 expression1_type;

   778         typedef const E2 expression2_type;

   779         typedef typename M::size_type size_type;

   780         typedef typename M::value_type value_type;

   782 #if BOOST_UBLAS_TYPE_CHECK

   783         matrix<value_type, column_major> cm (m);

   784         typedef typename type_traits<value_type>::real_type real_type;

   785         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));

   786         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());

   787 #endif

   788         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));

   789         for (size_type k = 0; k < size; ++ k) {

   790             vector<value_type> ce1 (column (e1 (), k));

   791             vector<value_type> re2 (row (e2 (), k));

   792             m.plus_assign (outer_prod (ce1, re2));

   793         }

   794 #if BOOST_UBLAS_TYPE_CHECK

   795         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());

   796 #endif

   797         return m;

   798     }

   800     // Dispatcher

   802   /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an

   803           optimized fashion.

   805           \param e1 the matrix expression \c A

   806           \param e2 the matrix expression \c X

   807           \param m  the result matrix \c M

   808           \param init a boolean parameter

   810           <tt>opb_prod(A, X, M, init)</tt> implements the well known

   811           axpy-product. Setting \a init to \c true is equivalent to call

   812           <tt>M.clear()</tt> before <tt>opb_prod</tt>. Currently \a init

   813           defaults to \c true, but this may change in the future.

   815           This function may give a speedup if \c A has less columns than

   816           rows, because the product is computed as a sum of outer

   817           products.

   819           \ingroup blas3

   821           \internal

   823           template parameters:

   824           \param M type of the result matrix \c M

   825           \param E1 type of a matrix expression \c A

   826           \param E2 type of a matrix expression \c X

   827   */

   828     template<class M, class E1, class E2>

   829     BOOST_UBLAS_INLINE

   830     M &

   831     opb_prod (const matrix_expression<E1> &e1,

   832               const matrix_expression<E2> &e2,

   833               M &m, bool init = true) {

   834         typedef typename M::value_type value_type;

   835         typedef typename M::storage_category storage_category;

   836         typedef typename M::orientation_category orientation_category;

   838         if (init)

   839             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));

   840         return opb_prod (e1, e2, m, storage_category (), orientation_category ());

   841     }

   842     template<class M, class E1, class E2>

   843     BOOST_UBLAS_INLINE

   844     M

   845     opb_prod (const matrix_expression<E1> &e1,

   846               const matrix_expression<E2> &e2) {

   847         typedef M matrix_type;

   849         matrix_type m (e1 ().size1 (), e2 ().size2 ());

   850         return opb_prod (e1, e2, m, true);

   851     }

   853 }}}

   855 #endif

author	sl
	Tue, 10 Jun 2014 14:32:02 +0200
changeset 1	260cb5ec6c19
permissions	-rw-r--r--