ceph/src/boost/libs/numeric/ublas/include/boost/numeric/ublas/operation.hpp

   1 //
   2 //  Copyright (c) 2000-2002
   3 //  Joerg Walter, Mathias Koch
   4 //
   5 //  Distributed under the Boost Software License, Version 1.0. (See
   6 //  accompanying file LICENSE_1_0.txt or copy at
   7 //  http://www.boost.org/LICENSE_1_0.txt)
   8 //
   9 //  The authors gratefully acknowledge the support of
  10 //  GeNeSys mbH & Co. KG in producing this work.
  11 //
  12
  13 #ifndef _BOOST_UBLAS_OPERATION_
  14 #define _BOOST_UBLAS_OPERATION_
  15
  16 #include <boost/numeric/ublas/matrix_proxy.hpp>
  17
  18 /** \file operation.hpp
  19  *  \brief This file contains some specialized products.
  20  */
  21
  22 // axpy-based products
  23 // Alexei Novakov had a lot of ideas to improve these. Thanks.
  24 // Hendrik Kueck proposed some new kernel. Thanks again.
  25
  26 namespace boost { namespace numeric { namespace ublas {
  27
  28     template<class V, class T1, class L1, class IA1, class TA1, class E2>
  29     BOOST_UBLAS_INLINE
  30     V &
  31     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
  32                const vector_expression<E2> &e2,
  33                V &v, row_major_tag) {
  34         typedef typename V::size_type size_type;
  35         typedef typename V::value_type value_type;
  36
  37         for (size_type i = 0; i < e1.filled1 () -1; ++ i) {
  38             size_type begin = e1.index1_data () [i];
  39             size_type end = e1.index1_data () [i + 1];
  40             value_type t (v (i));
  41             for (size_type j = begin; j < end; ++ j)
  42                 t += e1.value_data () [j] * e2 () (e1.index2_data () [j]);
  43             v (i) = t;
  44         }
  45         return v;
  46     }
  47
  48     template<class V, class T1, class L1, class IA1, class TA1, class E2>
  49     BOOST_UBLAS_INLINE
  50     V &
  51     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
  52                const vector_expression<E2> &e2,
  53                V &v, column_major_tag) {
  54         typedef typename V::size_type size_type;
  55
  56         for (size_type j = 0; j < e1.filled1 () -1; ++ j) {
  57             size_type begin = e1.index1_data () [j];
  58             size_type end = e1.index1_data () [j + 1];
  59             for (size_type i = begin; i < end; ++ i)
  60                 v (e1.index2_data () [i]) += e1.value_data () [i] * e2 () (j);
  61         }
  62         return v;
  63     }
  64
  65     // Dispatcher
  66     template<class V, class T1, class L1, class IA1, class TA1, class E2>
  67     BOOST_UBLAS_INLINE
  68     V &
  69     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
  70                const vector_expression<E2> &e2,
  71                V &v, bool init = true) {
  72         typedef typename V::value_type value_type;
  73         typedef typename L1::orientation_category orientation_category;
  74
  75         if (init)
  76             v.assign (zero_vector<value_type> (e1.size1 ()));
  77 #if BOOST_UBLAS_TYPE_CHECK
  78         vector<value_type> cv (v);
  79         typedef typename type_traits<value_type>::real_type real_type;
  80         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
  81         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
  82 #endif
  83         axpy_prod (e1, e2, v, orientation_category ());
  84 #if BOOST_UBLAS_TYPE_CHECK
  85         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
  86 #endif
  87         return v;
  88     }
  89     template<class V, class T1, class L1, class IA1, class TA1, class E2>
  90     BOOST_UBLAS_INLINE
  91     V
  92     axpy_prod (const compressed_matrix<T1, L1, 0, IA1, TA1> &e1,
  93                const vector_expression<E2> &e2) {
  94         typedef V vector_type;
  95
  96         vector_type v (e1.size1 ());
  97         return axpy_prod (e1, e2, v, true);
  98     }
  99
 100     template<class V, class T1, class L1, class IA1, class TA1, class E2>
 101     BOOST_UBLAS_INLINE
 102     V &
 103     axpy_prod (const coordinate_matrix<T1, L1, 0, IA1, TA1> &e1,
 104                const vector_expression<E2> &e2,
 105                V &v, bool init = true) {
 106         typedef typename V::size_type size_type;
 107         typedef typename V::value_type value_type;
 108         typedef L1 layout_type;
 109
 110         size_type size1 = e1.size1();
 111         size_type size2 = e1.size2();
 112
 113         if (init) {
 114             noalias(v) = zero_vector<value_type>(size1);
 115         }
 116
 117         for (size_type i = 0; i < e1.nnz(); ++i) {
 118             size_type row_index = layout_type::index_M( e1.index1_data () [i], e1.index2_data () [i] );
 119             size_type col_index = layout_type::index_m( e1.index1_data () [i], e1.index2_data () [i] );
 120             v( row_index ) += e1.value_data () [i] * e2 () (col_index);
 121         }
 122         return v;
 123     }
 124
 125     template<class V, class E1, class E2>
 126     BOOST_UBLAS_INLINE
 127     V &
 128     axpy_prod (const matrix_expression<E1> &e1,
 129                const vector_expression<E2> &e2,
 130                V &v, packed_random_access_iterator_tag, row_major_tag) {
 131         typedef const E1 expression1_type;
 132         typedef typename V::size_type size_type;
 133
 134         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
 135         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
 136         while (it1 != it1_end) {
 137             size_type index1 (it1.index1 ());
 138 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 139             typename expression1_type::const_iterator2 it2 (it1.begin ());
 140             typename expression1_type::const_iterator2 it2_end (it1.end ());
 141 #else
 142             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
 143             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
 144 #endif
 145             while (it2 != it2_end) {
 146                 v (index1) += *it2 * e2 () (it2.index2 ());
 147                 ++ it2;
 148             }
 149             ++ it1;
 150         }
 151         return v;
 152     }
 153
 154     template<class V, class E1, class E2>
 155     BOOST_UBLAS_INLINE
 156     V &
 157     axpy_prod (const matrix_expression<E1> &e1,
 158                const vector_expression<E2> &e2,
 159                V &v, packed_random_access_iterator_tag, column_major_tag) {
 160         typedef const E1 expression1_type;
 161         typedef typename V::size_type size_type;
 162
 163         typename expression1_type::const_iterator2 it2 (e1 ().begin2 ());
 164         typename expression1_type::const_iterator2 it2_end (e1 ().end2 ());
 165         while (it2 != it2_end) {
 166             size_type index2 (it2.index2 ());
 167 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 168             typename expression1_type::const_iterator1 it1 (it2.begin ());
 169             typename expression1_type::const_iterator1 it1_end (it2.end ());
 170 #else
 171             typename expression1_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
 172             typename expression1_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
 173 #endif
 174             while (it1 != it1_end) {
 175                 v (it1.index1 ()) += *it1 * e2 () (index2);
 176                 ++ it1;
 177             }
 178             ++ it2;
 179         }
 180         return v;
 181     }
 182
 183     template<class V, class E1, class E2>
 184     BOOST_UBLAS_INLINE
 185     V &
 186     axpy_prod (const matrix_expression<E1> &e1,
 187                const vector_expression<E2> &e2,
 188                V &v, sparse_bidirectional_iterator_tag) {
 189         typedef const E2 expression2_type;
 190
 191         typename expression2_type::const_iterator it (e2 ().begin ());
 192         typename expression2_type::const_iterator it_end (e2 ().end ());
 193         while (it != it_end) {
 194             v.plus_assign (column (e1 (), it.index ()) * *it);
 195             ++ it;
 196         }
 197         return v;
 198     }
 199
 200     // Dispatcher
 201     template<class V, class E1, class E2>
 202     BOOST_UBLAS_INLINE
 203     V &
 204     axpy_prod (const matrix_expression<E1> &e1,
 205                const vector_expression<E2> &e2,
 206                V &v, packed_random_access_iterator_tag) {
 207         typedef typename E1::orientation_category orientation_category;
 208         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
 209     }
 210
 211
 212   /** \brief computes <tt>v += A x</tt> or <tt>v = A x</tt> in an
 213           optimized fashion.
 214
 215           \param e1 the matrix expression \c A
 216           \param e2 the vector expression \c x
 217           \param v  the result vector \c v
 218           \param init a boolean parameter
 219
 220           <tt>axpy_prod(A, x, v, init)</tt> implements the well known
 221           axpy-product.  Setting \a init to \c true is equivalent to call
 222           <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
 223           defaults to \c true, but this may change in the future.
 224
 225           Up to now there are some specialisation for compressed
 226           matrices that give a large speed up compared to prod.
 227
 228           \ingroup blas2
 229
 230           \internal
 231
 232           template parameters:
 233           \param V type of the result vector \c v
 234           \param E1 type of a matrix expression \c A
 235           \param E2 type of a vector expression \c x
 236   */
 237     template<class V, class E1, class E2>
 238     BOOST_UBLAS_INLINE
 239     V &
 240     axpy_prod (const matrix_expression<E1> &e1,
 241                const vector_expression<E2> &e2,
 242                V &v, bool init = true) {
 243         typedef typename V::value_type value_type;
 244         typedef typename E2::const_iterator::iterator_category iterator_category;
 245
 246         if (init)
 247             v.assign (zero_vector<value_type> (e1 ().size1 ()));
 248 #if BOOST_UBLAS_TYPE_CHECK
 249         vector<value_type> cv (v);
 250         typedef typename type_traits<value_type>::real_type real_type;
 251         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
 252         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
 253 #endif
 254         axpy_prod (e1, e2, v, iterator_category ());
 255 #if BOOST_UBLAS_TYPE_CHECK
 256         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
 257 #endif
 258         return v;
 259     }
 260     template<class V, class E1, class E2>
 261     BOOST_UBLAS_INLINE
 262     V
 263     axpy_prod (const matrix_expression<E1> &e1,
 264                const vector_expression<E2> &e2) {
 265         typedef V vector_type;
 266
 267         vector_type v (e1 ().size1 ());
 268         return axpy_prod (e1, e2, v, true);
 269     }
 270
 271     template<class V, class E1, class T2, class IA2, class TA2>
 272     BOOST_UBLAS_INLINE
 273     V &
 274     axpy_prod (const vector_expression<E1> &e1,
 275                const compressed_matrix<T2, column_major, 0, IA2, TA2> &e2,
 276                V &v, column_major_tag) {
 277         typedef typename V::size_type size_type;
 278         typedef typename V::value_type value_type;
 279
 280         for (size_type j = 0; j < e2.filled1 () -1; ++ j) {
 281             size_type begin = e2.index1_data () [j];
 282             size_type end = e2.index1_data () [j + 1];
 283             value_type t (v (j));
 284             for (size_type i = begin; i < end; ++ i)
 285                 t += e2.value_data () [i] * e1 () (e2.index2_data () [i]);
 286             v (j) = t;
 287         }
 288         return v;
 289     }
 290
 291     template<class V, class E1, class T2, class IA2, class TA2>
 292     BOOST_UBLAS_INLINE
 293     V &
 294     axpy_prod (const vector_expression<E1> &e1,
 295                const compressed_matrix<T2, row_major, 0, IA2, TA2> &e2,
 296                V &v, row_major_tag) {
 297         typedef typename V::size_type size_type;
 298
 299         for (size_type i = 0; i < e2.filled1 () -1; ++ i) {
 300             size_type begin = e2.index1_data () [i];
 301             size_type end = e2.index1_data () [i + 1];
 302             for (size_type j = begin; j < end; ++ j)
 303                 v (e2.index2_data () [j]) += e2.value_data () [j] * e1 () (i);
 304         }
 305         return v;
 306     }
 307
 308     // Dispatcher
 309     template<class V, class E1, class T2, class L2, class IA2, class TA2>
 310     BOOST_UBLAS_INLINE
 311     V &
 312     axpy_prod (const vector_expression<E1> &e1,
 313                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2,
 314                V &v, bool init = true) {
 315         typedef typename V::value_type value_type;
 316         typedef typename L2::orientation_category orientation_category;
 317
 318         if (init)
 319             v.assign (zero_vector<value_type> (e2.size2 ()));
 320 #if BOOST_UBLAS_TYPE_CHECK
 321         vector<value_type> cv (v);
 322         typedef typename type_traits<value_type>::real_type real_type;
 323         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
 324         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
 325 #endif
 326         axpy_prod (e1, e2, v, orientation_category ());
 327 #if BOOST_UBLAS_TYPE_CHECK
 328         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
 329 #endif
 330         return v;
 331     }
 332     template<class V, class E1, class T2, class L2, class IA2, class TA2>
 333     BOOST_UBLAS_INLINE
 334     V
 335     axpy_prod (const vector_expression<E1> &e1,
 336                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2) {
 337         typedef V vector_type;
 338
 339         vector_type v (e2.size2 ());
 340         return axpy_prod (e1, e2, v, true);
 341     }
 342
 343     template<class V, class E1, class E2>
 344     BOOST_UBLAS_INLINE
 345     V &
 346     axpy_prod (const vector_expression<E1> &e1,
 347                const matrix_expression<E2> &e2,
 348                V &v, packed_random_access_iterator_tag, column_major_tag) {
 349         typedef const E2 expression2_type;
 350         typedef typename V::size_type size_type;
 351
 352         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
 353         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
 354         while (it2 != it2_end) {
 355             size_type index2 (it2.index2 ());
 356 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 357             typename expression2_type::const_iterator1 it1 (it2.begin ());
 358             typename expression2_type::const_iterator1 it1_end (it2.end ());
 359 #else
 360             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
 361             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
 362 #endif
 363             while (it1 != it1_end) {
 364                 v (index2) += *it1 * e1 () (it1.index1 ());
 365                 ++ it1;
 366             }
 367             ++ it2;
 368         }
 369         return v;
 370     }
 371
 372     template<class V, class E1, class E2>
 373     BOOST_UBLAS_INLINE
 374     V &
 375     axpy_prod (const vector_expression<E1> &e1,
 376                const matrix_expression<E2> &e2,
 377                V &v, packed_random_access_iterator_tag, row_major_tag) {
 378         typedef const E2 expression2_type;
 379         typedef typename V::size_type size_type;
 380
 381         typename expression2_type::const_iterator1 it1 (e2 ().begin1 ());
 382         typename expression2_type::const_iterator1 it1_end (e2 ().end1 ());
 383         while (it1 != it1_end) {
 384             size_type index1 (it1.index1 ());
 385 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 386             typename expression2_type::const_iterator2 it2 (it1.begin ());
 387             typename expression2_type::const_iterator2 it2_end (it1.end ());
 388 #else
 389             typename expression2_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
 390             typename expression2_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
 391 #endif
 392             while (it2 != it2_end) {
 393                 v (it2.index2 ()) += *it2 * e1 () (index1);
 394                 ++ it2;
 395             }
 396             ++ it1;
 397         }
 398         return v;
 399     }
 400
 401     template<class V, class E1, class E2>
 402     BOOST_UBLAS_INLINE
 403     V &
 404     axpy_prod (const vector_expression<E1> &e1,
 405                const matrix_expression<E2> &e2,
 406                V &v, sparse_bidirectional_iterator_tag) {
 407         typedef const E1 expression1_type;
 408
 409         typename expression1_type::const_iterator it (e1 ().begin ());
 410         typename expression1_type::const_iterator it_end (e1 ().end ());
 411         while (it != it_end) {
 412             v.plus_assign (*it * row (e2 (), it.index ()));
 413             ++ it;
 414         }
 415         return v;
 416     }
 417
 418     // Dispatcher
 419     template<class V, class E1, class E2>
 420     BOOST_UBLAS_INLINE
 421     V &
 422     axpy_prod (const vector_expression<E1> &e1,
 423                const matrix_expression<E2> &e2,
 424                V &v, packed_random_access_iterator_tag) {
 425         typedef typename E2::orientation_category orientation_category;
 426         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
 427     }
 428
 429
 430   /** \brief computes <tt>v += A<sup>T</sup> x</tt> or <tt>v = A<sup>T</sup> x</tt> in an
 431           optimized fashion.
 432
 433           \param e1 the vector expression \c x
 434           \param e2 the matrix expression \c A
 435           \param v  the result vector \c v
 436           \param init a boolean parameter
 437
 438           <tt>axpy_prod(x, A, v, init)</tt> implements the well known
 439           axpy-product.  Setting \a init to \c true is equivalent to call
 440           <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
 441           defaults to \c true, but this may change in the future.
 442
 443           Up to now there are some specialisation for compressed
 444           matrices that give a large speed up compared to prod.
 445
 446           \ingroup blas2
 447
 448           \internal
 449
 450           template parameters:
 451           \param V type of the result vector \c v
 452           \param E1 type of a vector expression \c x
 453           \param E2 type of a matrix expression \c A
 454   */
 455     template<class V, class E1, class E2>
 456     BOOST_UBLAS_INLINE
 457     V &
 458     axpy_prod (const vector_expression<E1> &e1,
 459                const matrix_expression<E2> &e2,
 460                V &v, bool init = true) {
 461         typedef typename V::value_type value_type;
 462         typedef typename E1::const_iterator::iterator_category iterator_category;
 463
 464         if (init)
 465             v.assign (zero_vector<value_type> (e2 ().size2 ()));
 466 #if BOOST_UBLAS_TYPE_CHECK
 467         vector<value_type> cv (v);
 468         typedef typename type_traits<value_type>::real_type real_type;
 469         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
 470         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
 471 #endif
 472         axpy_prod (e1, e2, v, iterator_category ());
 473 #if BOOST_UBLAS_TYPE_CHECK
 474         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
 475 #endif
 476         return v;
 477     }
 478     template<class V, class E1, class E2>
 479     BOOST_UBLAS_INLINE
 480     V
 481     axpy_prod (const vector_expression<E1> &e1,
 482                const matrix_expression<E2> &e2) {
 483         typedef V vector_type;
 484
 485         vector_type v (e2 ().size2 ());
 486         return axpy_prod (e1, e2, v, true);
 487     }
 488
 489     template<class M, class E1, class E2, class TRI>
 490     BOOST_UBLAS_INLINE
 491     M &
 492     axpy_prod (const matrix_expression<E1> &e1,
 493                const matrix_expression<E2> &e2,
 494                M &m, TRI,
 495                dense_proxy_tag, row_major_tag) {
 496
 497         typedef typename M::size_type size_type;
 498
 499 #if BOOST_UBLAS_TYPE_CHECK
 500         typedef typename M::value_type value_type;
 501         matrix<value_type, row_major> cm (m);
 502         typedef typename type_traits<value_type>::real_type real_type;
 503         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 504         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
 505 #endif
 506         size_type size1 (e1 ().size1 ());
 507         size_type size2 (e1 ().size2 ());
 508         for (size_type i = 0; i < size1; ++ i)
 509             for (size_type j = 0; j < size2; ++ j)
 510                 row (m, i).plus_assign (e1 () (i, j) * row (e2 (), j));
 511 #if BOOST_UBLAS_TYPE_CHECK
 512         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 513 #endif
 514         return m;
 515     }
 516     template<class M, class E1, class E2, class TRI>
 517     BOOST_UBLAS_INLINE
 518     M &
 519     axpy_prod (const matrix_expression<E1> &e1,
 520                const matrix_expression<E2> &e2,
 521                M &m, TRI,
 522                sparse_proxy_tag, row_major_tag) {
 523
 524         typedef TRI triangular_restriction;
 525         typedef const E1 expression1_type;
 526         typedef const E2 expression2_type;
 527
 528 #if BOOST_UBLAS_TYPE_CHECK
 529         typedef typename M::value_type value_type;
 530         matrix<value_type, row_major> cm (m);
 531         typedef typename type_traits<value_type>::real_type real_type;
 532         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 533         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
 534 #endif
 535         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
 536         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
 537         while (it1 != it1_end) {
 538 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 539             typename expression1_type::const_iterator2 it2 (it1.begin ());
 540             typename expression1_type::const_iterator2 it2_end (it1.end ());
 541 #else
 542             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
 543             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
 544 #endif
 545             while (it2 != it2_end) {
 546                 // row (m, it1.index1 ()).plus_assign (*it2 * row (e2 (), it2.index2 ()));
 547                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
 548                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
 549                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
 550                 while (itr != itr_end) {
 551                     if (triangular_restriction::other (it1.index1 (), itr.index ()))
 552                         m (it1.index1 (), itr.index ()) += *it2 * *itr;
 553                     ++ itr;
 554                 }
 555                 ++ it2;
 556             }
 557             ++ it1;
 558         }
 559 #if BOOST_UBLAS_TYPE_CHECK
 560         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 561 #endif
 562         return m;
 563     }
 564
 565     template<class M, class E1, class E2, class TRI>
 566     BOOST_UBLAS_INLINE
 567     M &
 568     axpy_prod (const matrix_expression<E1> &e1,
 569                const matrix_expression<E2> &e2,
 570                M &m, TRI,
 571                dense_proxy_tag, column_major_tag) {
 572         typedef typename M::size_type size_type;
 573
 574 #if BOOST_UBLAS_TYPE_CHECK
 575         typedef typename M::value_type value_type;
 576         matrix<value_type, column_major> cm (m);
 577         typedef typename type_traits<value_type>::real_type real_type;
 578         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 579         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
 580 #endif
 581         size_type size1 (e2 ().size1 ());
 582         size_type size2 (e2 ().size2 ());
 583         for (size_type j = 0; j < size2; ++ j)
 584             for (size_type i = 0; i < size1; ++ i)
 585                 column (m, j).plus_assign (e2 () (i, j) * column (e1 (), i));
 586 #if BOOST_UBLAS_TYPE_CHECK
 587         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 588 #endif
 589         return m;
 590     }
 591     template<class M, class E1, class E2, class TRI>
 592     BOOST_UBLAS_INLINE
 593     M &
 594     axpy_prod (const matrix_expression<E1> &e1,
 595                const matrix_expression<E2> &e2,
 596                M &m, TRI,
 597                sparse_proxy_tag, column_major_tag) {
 598         typedef TRI triangular_restriction;
 599         typedef const E1 expression1_type;
 600         typedef const E2 expression2_type;
 601
 602
 603 #if BOOST_UBLAS_TYPE_CHECK
 604         typedef typename M::value_type value_type;
 605         matrix<value_type, column_major> cm (m);
 606         typedef typename type_traits<value_type>::real_type real_type;
 607         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 608         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
 609 #endif
 610         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
 611         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
 612         while (it2 != it2_end) {
 613 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 614             typename expression2_type::const_iterator1 it1 (it2.begin ());
 615             typename expression2_type::const_iterator1 it1_end (it2.end ());
 616 #else
 617             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
 618             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
 619 #endif
 620             while (it1 != it1_end) {
 621                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
 622                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
 623                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
 624                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
 625                 while (itc != itc_end) {
 626                     if(triangular_restriction::other (itc.index (), it2.index2 ()))
 627                        m (itc.index (), it2.index2 ()) += *it1 * *itc;
 628                     ++ itc;
 629                 }
 630                 ++ it1;
 631             }
 632             ++ it2;
 633         }
 634 #if BOOST_UBLAS_TYPE_CHECK
 635         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 636 #endif
 637         return m;
 638     }
 639
 640     // Dispatcher
 641     template<class M, class E1, class E2, class TRI>
 642     BOOST_UBLAS_INLINE
 643     M &
 644     axpy_prod (const matrix_expression<E1> &e1,
 645                const matrix_expression<E2> &e2,
 646                M &m, TRI, bool init = true) {
 647         typedef typename M::value_type value_type;
 648         typedef typename M::storage_category storage_category;
 649         typedef typename M::orientation_category orientation_category;
 650         typedef TRI triangular_restriction;
 651
 652         if (init)
 653             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
 654         return axpy_prod (e1, e2, m, triangular_restriction (), storage_category (), orientation_category ());
 655     }
 656     template<class M, class E1, class E2, class TRI>
 657     BOOST_UBLAS_INLINE
 658     M
 659     axpy_prod (const matrix_expression<E1> &e1,
 660                const matrix_expression<E2> &e2,
 661                TRI) {
 662         typedef M matrix_type;
 663         typedef TRI triangular_restriction;
 664
 665         matrix_type m (e1 ().size1 (), e2 ().size2 ());
 666         return axpy_prod (e1, e2, m, triangular_restriction (), true);
 667     }
 668
 669   /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
 670           optimized fashion.
 671
 672           \param e1 the matrix expression \c A
 673           \param e2 the matrix expression \c X
 674           \param m  the result matrix \c M
 675           \param init a boolean parameter
 676
 677           <tt>axpy_prod(A, X, M, init)</tt> implements the well known
 678           axpy-product.  Setting \a init to \c true is equivalent to call
 679           <tt>M.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
 680           defaults to \c true, but this may change in the future.
 681
 682           Up to now there are no specialisations.
 683
 684           \ingroup blas3
 685
 686           \internal
 687
 688           template parameters:
 689           \param M type of the result matrix \c M
 690           \param E1 type of a matrix expression \c A
 691           \param E2 type of a matrix expression \c X
 692   */
 693     template<class M, class E1, class E2>
 694     BOOST_UBLAS_INLINE
 695     M &
 696     axpy_prod (const matrix_expression<E1> &e1,
 697                const matrix_expression<E2> &e2,
 698                M &m, bool init = true) {
 699         typedef typename M::value_type value_type;
 700         typedef typename M::storage_category storage_category;
 701         typedef typename M::orientation_category orientation_category;
 702
 703         if (init)
 704             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
 705         return axpy_prod (e1, e2, m, full (), storage_category (), orientation_category ());
 706     }
 707     template<class M, class E1, class E2>
 708     BOOST_UBLAS_INLINE
 709     M
 710     axpy_prod (const matrix_expression<E1> &e1,
 711                const matrix_expression<E2> &e2) {
 712         typedef M matrix_type;
 713
 714         matrix_type m (e1 ().size1 (), e2 ().size2 ());
 715         return axpy_prod (e1, e2, m, full (), true);
 716     }
 717
 718
 719     template<class M, class E1, class E2>
 720     BOOST_UBLAS_INLINE
 721     M &
 722     opb_prod (const matrix_expression<E1> &e1,
 723               const matrix_expression<E2> &e2,
 724               M &m,
 725               dense_proxy_tag, row_major_tag) {
 726         typedef typename M::size_type size_type;
 727         typedef typename M::value_type value_type;
 728
 729 #if BOOST_UBLAS_TYPE_CHECK
 730         matrix<value_type, row_major> cm (m);
 731         typedef typename type_traits<value_type>::real_type real_type;
 732         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 733         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
 734 #endif
 735         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
 736         for (size_type k = 0; k < size; ++ k) {
 737             vector<value_type> ce1 (column (e1 (), k));
 738             vector<value_type> re2 (row (e2 (), k));
 739             m.plus_assign (outer_prod (ce1, re2));
 740         }
 741 #if BOOST_UBLAS_TYPE_CHECK
 742         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 743 #endif
 744         return m;
 745     }
 746
 747     template<class M, class E1, class E2>
 748     BOOST_UBLAS_INLINE
 749     M &
 750     opb_prod (const matrix_expression<E1> &e1,
 751               const matrix_expression<E2> &e2,
 752               M &m,
 753               dense_proxy_tag, column_major_tag) {
 754         typedef typename M::size_type size_type;
 755         typedef typename M::value_type value_type;
 756
 757 #if BOOST_UBLAS_TYPE_CHECK
 758         matrix<value_type, column_major> cm (m);
 759         typedef typename type_traits<value_type>::real_type real_type;
 760         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 761         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
 762 #endif
 763         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
 764         for (size_type k = 0; k < size; ++ k) {
 765             vector<value_type> ce1 (column (e1 (), k));
 766             vector<value_type> re2 (row (e2 (), k));
 767             m.plus_assign (outer_prod (ce1, re2));
 768         }
 769 #if BOOST_UBLAS_TYPE_CHECK
 770         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 771 #endif
 772         return m;
 773     }
 774
 775     // Dispatcher
 776
 777   /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
 778           optimized fashion.
 779
 780           \param e1 the matrix expression \c A
 781           \param e2 the matrix expression \c X
 782           \param m  the result matrix \c M
 783           \param init a boolean parameter
 784
 785           <tt>opb_prod(A, X, M, init)</tt> implements the well known
 786           axpy-product. Setting \a init to \c true is equivalent to call
 787           <tt>M.clear()</tt> before <tt>opb_prod</tt>. Currently \a init
 788           defaults to \c true, but this may change in the future.
 789
 790           This function may give a speedup if \c A has less columns than
 791           rows, because the product is computed as a sum of outer
 792           products.
 793
 794           \ingroup blas3
 795
 796           \internal
 797
 798           template parameters:
 799           \param M type of the result matrix \c M
 800           \param E1 type of a matrix expression \c A
 801           \param E2 type of a matrix expression \c X
 802   */
 803     template<class M, class E1, class E2>
 804     BOOST_UBLAS_INLINE
 805     M &
 806     opb_prod (const matrix_expression<E1> &e1,
 807               const matrix_expression<E2> &e2,
 808               M &m, bool init = true) {
 809         typedef typename M::value_type value_type;
 810         typedef typename M::storage_category storage_category;
 811         typedef typename M::orientation_category orientation_category;
 812
 813         if (init)
 814             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
 815         return opb_prod (e1, e2, m, storage_category (), orientation_category ());
 816     }
 817     template<class M, class E1, class E2>
 818     BOOST_UBLAS_INLINE
 819     M
 820     opb_prod (const matrix_expression<E1> &e1,
 821               const matrix_expression<E2> &e2) {
 822         typedef M matrix_type;
 823
 824         matrix_type m (e1 ().size1 (), e2 ().size2 ());
 825         return opb_prod (e1, e2, m, true);
 826     }
 827
 828 }}}
 829
 830 #endif