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 const E2 expression2_type;
 133         typedef typename V::size_type size_type;
 134
 135         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
 136         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
 137         while (it1 != it1_end) {
 138             size_type index1 (it1.index1 ());
 139 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 140             typename expression1_type::const_iterator2 it2 (it1.begin ());
 141             typename expression1_type::const_iterator2 it2_end (it1.end ());
 142 #else
 143             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
 144             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
 145 #endif
 146             while (it2 != it2_end) {
 147                 v (index1) += *it2 * e2 () (it2.index2 ());
 148                 ++ it2;
 149             }
 150             ++ it1;
 151         }
 152         return v;
 153     }
 154
 155     template<class V, class E1, class E2>
 156     BOOST_UBLAS_INLINE
 157     V &
 158     axpy_prod (const matrix_expression<E1> &e1,
 159                const vector_expression<E2> &e2,
 160                V &v, packed_random_access_iterator_tag, column_major_tag) {
 161         typedef const E1 expression1_type;
 162         typedef const E2 expression2_type;
 163         typedef typename V::size_type size_type;
 164
 165         typename expression1_type::const_iterator2 it2 (e1 ().begin2 ());
 166         typename expression1_type::const_iterator2 it2_end (e1 ().end2 ());
 167         while (it2 != it2_end) {
 168             size_type index2 (it2.index2 ());
 169 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 170             typename expression1_type::const_iterator1 it1 (it2.begin ());
 171             typename expression1_type::const_iterator1 it1_end (it2.end ());
 172 #else
 173             typename expression1_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
 174             typename expression1_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
 175 #endif
 176             while (it1 != it1_end) {
 177                 v (it1.index1 ()) += *it1 * e2 () (index2);
 178                 ++ it1;
 179             }
 180             ++ it2;
 181         }
 182         return v;
 183     }
 184
 185     template<class V, class E1, class E2>
 186     BOOST_UBLAS_INLINE
 187     V &
 188     axpy_prod (const matrix_expression<E1> &e1,
 189                const vector_expression<E2> &e2,
 190                V &v, sparse_bidirectional_iterator_tag) {
 191         typedef const E1 expression1_type;
 192         typedef const E2 expression2_type;
 193         typedef typename V::size_type size_type;
 194
 195         typename expression2_type::const_iterator it (e2 ().begin ());
 196         typename expression2_type::const_iterator it_end (e2 ().end ());
 197         while (it != it_end) {
 198             v.plus_assign (column (e1 (), it.index ()) * *it);
 199             ++ it;
 200         }
 201         return v;
 202     }
 203
 204     // Dispatcher
 205     template<class V, class E1, class E2>
 206     BOOST_UBLAS_INLINE
 207     V &
 208     axpy_prod (const matrix_expression<E1> &e1,
 209                const vector_expression<E2> &e2,
 210                V &v, packed_random_access_iterator_tag) {
 211         typedef typename E1::orientation_category orientation_category;
 212         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
 213     }
 214
 215
 216   /** \brief computes <tt>v += A x</tt> or <tt>v = A x</tt> in an
 217           optimized fashion.
 218
 219           \param e1 the matrix expression \c A
 220           \param e2 the vector expression \c x
 221           \param v  the result vector \c v
 222           \param init a boolean parameter
 223
 224           <tt>axpy_prod(A, x, v, init)</tt> implements the well known
 225           axpy-product.  Setting \a init to \c true is equivalent to call
 226           <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
 227           defaults to \c true, but this may change in the future.
 228
 229           Up to now there are some specialisation for compressed
 230           matrices that give a large speed up compared to prod.
 231
 232           \ingroup blas2
 233
 234           \internal
 235
 236           template parameters:
 237           \param V type of the result vector \c v
 238           \param E1 type of a matrix expression \c A
 239           \param E2 type of a vector expression \c x
 240   */
 241     template<class V, class E1, class E2>
 242     BOOST_UBLAS_INLINE
 243     V &
 244     axpy_prod (const matrix_expression<E1> &e1,
 245                const vector_expression<E2> &e2,
 246                V &v, bool init = true) {
 247         typedef typename V::value_type value_type;
 248         typedef typename E2::const_iterator::iterator_category iterator_category;
 249
 250         if (init)
 251             v.assign (zero_vector<value_type> (e1 ().size1 ()));
 252 #if BOOST_UBLAS_TYPE_CHECK
 253         vector<value_type> cv (v);
 254         typedef typename type_traits<value_type>::real_type real_type;
 255         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
 256         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
 257 #endif
 258         axpy_prod (e1, e2, v, iterator_category ());
 259 #if BOOST_UBLAS_TYPE_CHECK
 260         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
 261 #endif
 262         return v;
 263     }
 264     template<class V, class E1, class E2>
 265     BOOST_UBLAS_INLINE
 266     V
 267     axpy_prod (const matrix_expression<E1> &e1,
 268                const vector_expression<E2> &e2) {
 269         typedef V vector_type;
 270
 271         vector_type v (e1 ().size1 ());
 272         return axpy_prod (e1, e2, v, true);
 273     }
 274
 275     template<class V, class E1, class T2, class IA2, class TA2>
 276     BOOST_UBLAS_INLINE
 277     V &
 278     axpy_prod (const vector_expression<E1> &e1,
 279                const compressed_matrix<T2, column_major, 0, IA2, TA2> &e2,
 280                V &v, column_major_tag) {
 281         typedef typename V::size_type size_type;
 282         typedef typename V::value_type value_type;
 283
 284         for (size_type j = 0; j < e2.filled1 () -1; ++ j) {
 285             size_type begin = e2.index1_data () [j];
 286             size_type end = e2.index1_data () [j + 1];
 287             value_type t (v (j));
 288             for (size_type i = begin; i < end; ++ i)
 289                 t += e2.value_data () [i] * e1 () (e2.index2_data () [i]);
 290             v (j) = t;
 291         }
 292         return v;
 293     }
 294
 295     template<class V, class E1, class T2, class IA2, class TA2>
 296     BOOST_UBLAS_INLINE
 297     V &
 298     axpy_prod (const vector_expression<E1> &e1,
 299                const compressed_matrix<T2, row_major, 0, IA2, TA2> &e2,
 300                V &v, row_major_tag) {
 301         typedef typename V::size_type size_type;
 302
 303         for (size_type i = 0; i < e2.filled1 () -1; ++ i) {
 304             size_type begin = e2.index1_data () [i];
 305             size_type end = e2.index1_data () [i + 1];
 306             for (size_type j = begin; j < end; ++ j)
 307                 v (e2.index2_data () [j]) += e2.value_data () [j] * e1 () (i);
 308         }
 309         return v;
 310     }
 311
 312     // Dispatcher
 313     template<class V, class E1, class T2, class L2, class IA2, class TA2>
 314     BOOST_UBLAS_INLINE
 315     V &
 316     axpy_prod (const vector_expression<E1> &e1,
 317                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2,
 318                V &v, bool init = true) {
 319         typedef typename V::value_type value_type;
 320         typedef typename L2::orientation_category orientation_category;
 321
 322         if (init)
 323             v.assign (zero_vector<value_type> (e2.size2 ()));
 324 #if BOOST_UBLAS_TYPE_CHECK
 325         vector<value_type> cv (v);
 326         typedef typename type_traits<value_type>::real_type real_type;
 327         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
 328         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
 329 #endif
 330         axpy_prod (e1, e2, v, orientation_category ());
 331 #if BOOST_UBLAS_TYPE_CHECK
 332         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
 333 #endif
 334         return v;
 335     }
 336     template<class V, class E1, class T2, class L2, class IA2, class TA2>
 337     BOOST_UBLAS_INLINE
 338     V
 339     axpy_prod (const vector_expression<E1> &e1,
 340                const compressed_matrix<T2, L2, 0, IA2, TA2> &e2) {
 341         typedef V vector_type;
 342
 343         vector_type v (e2.size2 ());
 344         return axpy_prod (e1, e2, v, true);
 345     }
 346
 347     template<class V, class E1, class E2>
 348     BOOST_UBLAS_INLINE
 349     V &
 350     axpy_prod (const vector_expression<E1> &e1,
 351                const matrix_expression<E2> &e2,
 352                V &v, packed_random_access_iterator_tag, column_major_tag) {
 353         typedef const E1 expression1_type;
 354         typedef const E2 expression2_type;
 355         typedef typename V::size_type size_type;
 356
 357         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
 358         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
 359         while (it2 != it2_end) {
 360             size_type index2 (it2.index2 ());
 361 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 362             typename expression2_type::const_iterator1 it1 (it2.begin ());
 363             typename expression2_type::const_iterator1 it1_end (it2.end ());
 364 #else
 365             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
 366             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
 367 #endif
 368             while (it1 != it1_end) {
 369                 v (index2) += *it1 * e1 () (it1.index1 ());
 370                 ++ it1;
 371             }
 372             ++ it2;
 373         }
 374         return v;
 375     }
 376
 377     template<class V, class E1, class E2>
 378     BOOST_UBLAS_INLINE
 379     V &
 380     axpy_prod (const vector_expression<E1> &e1,
 381                const matrix_expression<E2> &e2,
 382                V &v, packed_random_access_iterator_tag, row_major_tag) {
 383         typedef const E1 expression1_type;
 384         typedef const E2 expression2_type;
 385         typedef typename V::size_type size_type;
 386
 387         typename expression2_type::const_iterator1 it1 (e2 ().begin1 ());
 388         typename expression2_type::const_iterator1 it1_end (e2 ().end1 ());
 389         while (it1 != it1_end) {
 390             size_type index1 (it1.index1 ());
 391 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 392             typename expression2_type::const_iterator2 it2 (it1.begin ());
 393             typename expression2_type::const_iterator2 it2_end (it1.end ());
 394 #else
 395             typename expression2_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
 396             typename expression2_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
 397 #endif
 398             while (it2 != it2_end) {
 399                 v (it2.index2 ()) += *it2 * e1 () (index1);
 400                 ++ it2;
 401             }
 402             ++ it1;
 403         }
 404         return v;
 405     }
 406
 407     template<class V, class E1, class E2>
 408     BOOST_UBLAS_INLINE
 409     V &
 410     axpy_prod (const vector_expression<E1> &e1,
 411                const matrix_expression<E2> &e2,
 412                V &v, sparse_bidirectional_iterator_tag) {
 413         typedef const E1 expression1_type;
 414         typedef const E2 expression2_type;
 415         typedef typename V::size_type size_type;
 416
 417         typename expression1_type::const_iterator it (e1 ().begin ());
 418         typename expression1_type::const_iterator it_end (e1 ().end ());
 419         while (it != it_end) {
 420             v.plus_assign (*it * row (e2 (), it.index ()));
 421             ++ it;
 422         }
 423         return v;
 424     }
 425
 426     // Dispatcher
 427     template<class V, class E1, class E2>
 428     BOOST_UBLAS_INLINE
 429     V &
 430     axpy_prod (const vector_expression<E1> &e1,
 431                const matrix_expression<E2> &e2,
 432                V &v, packed_random_access_iterator_tag) {
 433         typedef typename E2::orientation_category orientation_category;
 434         return axpy_prod (e1, e2, v, packed_random_access_iterator_tag (), orientation_category ());
 435     }
 436
 437
 438   /** \brief computes <tt>v += A<sup>T</sup> x</tt> or <tt>v = A<sup>T</sup> x</tt> in an
 439           optimized fashion.
 440
 441           \param e1 the vector expression \c x
 442           \param e2 the matrix expression \c A
 443           \param v  the result vector \c v
 444           \param init a boolean parameter
 445
 446           <tt>axpy_prod(x, A, v, init)</tt> implements the well known
 447           axpy-product.  Setting \a init to \c true is equivalent to call
 448           <tt>v.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
 449           defaults to \c true, but this may change in the future.
 450
 451           Up to now there are some specialisation for compressed
 452           matrices that give a large speed up compared to prod.
 453
 454           \ingroup blas2
 455
 456           \internal
 457
 458           template parameters:
 459           \param V type of the result vector \c v
 460           \param E1 type of a vector expression \c x
 461           \param E2 type of a matrix expression \c A
 462   */
 463     template<class V, class E1, class E2>
 464     BOOST_UBLAS_INLINE
 465     V &
 466     axpy_prod (const vector_expression<E1> &e1,
 467                const matrix_expression<E2> &e2,
 468                V &v, bool init = true) {
 469         typedef typename V::value_type value_type;
 470         typedef typename E1::const_iterator::iterator_category iterator_category;
 471
 472         if (init)
 473             v.assign (zero_vector<value_type> (e2 ().size2 ()));
 474 #if BOOST_UBLAS_TYPE_CHECK
 475         vector<value_type> cv (v);
 476         typedef typename type_traits<value_type>::real_type real_type;
 477         real_type verrorbound (norm_1 (v) + norm_1 (e1) * norm_1 (e2));
 478         indexing_vector_assign<scalar_plus_assign> (cv, prod (e1, e2));
 479 #endif
 480         axpy_prod (e1, e2, v, iterator_category ());
 481 #if BOOST_UBLAS_TYPE_CHECK
 482         BOOST_UBLAS_CHECK (norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound, internal_logic ());
 483 #endif
 484         return v;
 485     }
 486     template<class V, class E1, class E2>
 487     BOOST_UBLAS_INLINE
 488     V
 489     axpy_prod (const vector_expression<E1> &e1,
 490                const matrix_expression<E2> &e2) {
 491         typedef V vector_type;
 492
 493         vector_type v (e2 ().size2 ());
 494         return axpy_prod (e1, e2, v, true);
 495     }
 496
 497     template<class M, class E1, class E2, class TRI>
 498     BOOST_UBLAS_INLINE
 499     M &
 500     axpy_prod (const matrix_expression<E1> &e1,
 501                const matrix_expression<E2> &e2,
 502                M &m, TRI,
 503                dense_proxy_tag, row_major_tag) {
 504         typedef M matrix_type;
 505         typedef const E1 expression1_type;
 506         typedef const E2 expression2_type;
 507         typedef typename M::size_type size_type;
 508         typedef typename M::value_type value_type;
 509
 510 #if BOOST_UBLAS_TYPE_CHECK
 511         matrix<value_type, row_major> cm (m);
 512         typedef typename type_traits<value_type>::real_type real_type;
 513         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 514         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
 515 #endif
 516         size_type size1 (e1 ().size1 ());
 517         size_type size2 (e1 ().size2 ());
 518         for (size_type i = 0; i < size1; ++ i)
 519             for (size_type j = 0; j < size2; ++ j)
 520                 row (m, i).plus_assign (e1 () (i, j) * row (e2 (), j));
 521 #if BOOST_UBLAS_TYPE_CHECK
 522         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 523 #endif
 524         return m;
 525     }
 526     template<class M, class E1, class E2, class TRI>
 527     BOOST_UBLAS_INLINE
 528     M &
 529     axpy_prod (const matrix_expression<E1> &e1,
 530                const matrix_expression<E2> &e2,
 531                M &m, TRI,
 532                sparse_proxy_tag, row_major_tag) {
 533         typedef M matrix_type;
 534         typedef TRI triangular_restriction;
 535         typedef const E1 expression1_type;
 536         typedef const E2 expression2_type;
 537         typedef typename M::size_type size_type;
 538         typedef typename M::value_type value_type;
 539
 540 #if BOOST_UBLAS_TYPE_CHECK
 541         matrix<value_type, row_major> cm (m);
 542         typedef typename type_traits<value_type>::real_type real_type;
 543         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 544         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
 545 #endif
 546         typename expression1_type::const_iterator1 it1 (e1 ().begin1 ());
 547         typename expression1_type::const_iterator1 it1_end (e1 ().end1 ());
 548         while (it1 != it1_end) {
 549 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 550             typename expression1_type::const_iterator2 it2 (it1.begin ());
 551             typename expression1_type::const_iterator2 it2_end (it1.end ());
 552 #else
 553             typename expression1_type::const_iterator2 it2 (boost::numeric::ublas::begin (it1, iterator1_tag ()));
 554             typename expression1_type::const_iterator2 it2_end (boost::numeric::ublas::end (it1, iterator1_tag ()));
 555 #endif
 556             while (it2 != it2_end) {
 557                 // row (m, it1.index1 ()).plus_assign (*it2 * row (e2 (), it2.index2 ()));
 558                 matrix_row<expression2_type> mr (e2 (), it2.index2 ());
 559                 typename matrix_row<expression2_type>::const_iterator itr (mr.begin ());
 560                 typename matrix_row<expression2_type>::const_iterator itr_end (mr.end ());
 561                 while (itr != itr_end) {
 562                     if (triangular_restriction::other (it1.index1 (), itr.index ()))
 563                         m (it1.index1 (), itr.index ()) += *it2 * *itr;
 564                     ++ itr;
 565                 }
 566                 ++ it2;
 567             }
 568             ++ it1;
 569         }
 570 #if BOOST_UBLAS_TYPE_CHECK
 571         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 572 #endif
 573         return m;
 574     }
 575
 576     template<class M, class E1, class E2, class TRI>
 577     BOOST_UBLAS_INLINE
 578     M &
 579     axpy_prod (const matrix_expression<E1> &e1,
 580                const matrix_expression<E2> &e2,
 581                M &m, TRI,
 582                dense_proxy_tag, column_major_tag) {
 583         typedef M matrix_type;
 584         typedef const E1 expression1_type;
 585         typedef const E2 expression2_type;
 586         typedef typename M::size_type size_type;
 587         typedef typename M::value_type value_type;
 588
 589 #if BOOST_UBLAS_TYPE_CHECK
 590         matrix<value_type, column_major> cm (m);
 591         typedef typename type_traits<value_type>::real_type real_type;
 592         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 593         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
 594 #endif
 595         size_type size1 (e2 ().size1 ());
 596         size_type size2 (e2 ().size2 ());
 597         for (size_type j = 0; j < size2; ++ j)
 598             for (size_type i = 0; i < size1; ++ i)
 599                 column (m, j).plus_assign (e2 () (i, j) * column (e1 (), i));
 600 #if BOOST_UBLAS_TYPE_CHECK
 601         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 602 #endif
 603         return m;
 604     }
 605     template<class M, class E1, class E2, class TRI>
 606     BOOST_UBLAS_INLINE
 607     M &
 608     axpy_prod (const matrix_expression<E1> &e1,
 609                const matrix_expression<E2> &e2,
 610                M &m, TRI,
 611                sparse_proxy_tag, column_major_tag) {
 612         typedef M matrix_type;
 613         typedef TRI triangular_restriction;
 614         typedef const E1 expression1_type;
 615         typedef const E2 expression2_type;
 616         typedef typename M::size_type size_type;
 617         typedef typename M::value_type value_type;
 618
 619 #if BOOST_UBLAS_TYPE_CHECK
 620         matrix<value_type, column_major> cm (m);
 621         typedef typename type_traits<value_type>::real_type real_type;
 622         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 623         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
 624 #endif
 625         typename expression2_type::const_iterator2 it2 (e2 ().begin2 ());
 626         typename expression2_type::const_iterator2 it2_end (e2 ().end2 ());
 627         while (it2 != it2_end) {
 628 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
 629             typename expression2_type::const_iterator1 it1 (it2.begin ());
 630             typename expression2_type::const_iterator1 it1_end (it2.end ());
 631 #else
 632             typename expression2_type::const_iterator1 it1 (boost::numeric::ublas::begin (it2, iterator2_tag ()));
 633             typename expression2_type::const_iterator1 it1_end (boost::numeric::ublas::end (it2, iterator2_tag ()));
 634 #endif
 635             while (it1 != it1_end) {
 636                 // column (m, it2.index2 ()).plus_assign (*it1 * column (e1 (), it1.index1 ()));
 637                 matrix_column<expression1_type> mc (e1 (), it1.index1 ());
 638                 typename matrix_column<expression1_type>::const_iterator itc (mc.begin ());
 639                 typename matrix_column<expression1_type>::const_iterator itc_end (mc.end ());
 640                 while (itc != itc_end) {
 641                     if(triangular_restriction::other (itc.index (), it2.index2 ()))
 642                        m (itc.index (), it2.index2 ()) += *it1 * *itc;
 643                     ++ itc;
 644                 }
 645                 ++ it1;
 646             }
 647             ++ it2;
 648         }
 649 #if BOOST_UBLAS_TYPE_CHECK
 650         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 651 #endif
 652         return m;
 653     }
 654
 655     // Dispatcher
 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                M &m, TRI, bool init = true) {
 662         typedef typename M::value_type value_type;
 663         typedef typename M::storage_category storage_category;
 664         typedef typename M::orientation_category orientation_category;
 665         typedef TRI triangular_restriction;
 666
 667         if (init)
 668             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
 669         return axpy_prod (e1, e2, m, triangular_restriction (), storage_category (), orientation_category ());
 670     }
 671     template<class M, class E1, class E2, class TRI>
 672     BOOST_UBLAS_INLINE
 673     M
 674     axpy_prod (const matrix_expression<E1> &e1,
 675                const matrix_expression<E2> &e2,
 676                TRI) {
 677         typedef M matrix_type;
 678         typedef TRI triangular_restriction;
 679
 680         matrix_type m (e1 ().size1 (), e2 ().size2 ());
 681         return axpy_prod (e1, e2, m, triangular_restriction (), true);
 682     }
 683
 684   /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
 685           optimized fashion.
 686
 687           \param e1 the matrix expression \c A
 688           \param e2 the matrix expression \c X
 689           \param m  the result matrix \c M
 690           \param init a boolean parameter
 691
 692           <tt>axpy_prod(A, X, M, init)</tt> implements the well known
 693           axpy-product.  Setting \a init to \c true is equivalent to call
 694           <tt>M.clear()</tt> before <tt>axpy_prod</tt>. Currently \a init
 695           defaults to \c true, but this may change in the future.
 696
 697           Up to now there are no specialisations.
 698
 699           \ingroup blas3
 700
 701           \internal
 702
 703           template parameters:
 704           \param M type of the result matrix \c M
 705           \param E1 type of a matrix expression \c A
 706           \param E2 type of a matrix expression \c X
 707   */
 708     template<class M, class E1, class E2>
 709     BOOST_UBLAS_INLINE
 710     M &
 711     axpy_prod (const matrix_expression<E1> &e1,
 712                const matrix_expression<E2> &e2,
 713                M &m, bool init = true) {
 714         typedef typename M::value_type value_type;
 715         typedef typename M::storage_category storage_category;
 716         typedef typename M::orientation_category orientation_category;
 717
 718         if (init)
 719             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
 720         return axpy_prod (e1, e2, m, full (), storage_category (), orientation_category ());
 721     }
 722     template<class M, class E1, class E2>
 723     BOOST_UBLAS_INLINE
 724     M
 725     axpy_prod (const matrix_expression<E1> &e1,
 726                const matrix_expression<E2> &e2) {
 727         typedef M matrix_type;
 728
 729         matrix_type m (e1 ().size1 (), e2 ().size2 ());
 730         return axpy_prod (e1, e2, m, full (), true);
 731     }
 732
 733
 734     template<class M, class E1, class E2>
 735     BOOST_UBLAS_INLINE
 736     M &
 737     opb_prod (const matrix_expression<E1> &e1,
 738               const matrix_expression<E2> &e2,
 739               M &m,
 740               dense_proxy_tag, row_major_tag) {
 741         typedef M matrix_type;
 742         typedef const E1 expression1_type;
 743         typedef const E2 expression2_type;
 744         typedef typename M::size_type size_type;
 745         typedef typename M::value_type value_type;
 746
 747 #if BOOST_UBLAS_TYPE_CHECK
 748         matrix<value_type, row_major> cm (m);
 749         typedef typename type_traits<value_type>::real_type real_type;
 750         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 751         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), row_major_tag ());
 752 #endif
 753         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
 754         for (size_type k = 0; k < size; ++ k) {
 755             vector<value_type> ce1 (column (e1 (), k));
 756             vector<value_type> re2 (row (e2 (), k));
 757             m.plus_assign (outer_prod (ce1, re2));
 758         }
 759 #if BOOST_UBLAS_TYPE_CHECK
 760         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 761 #endif
 762         return m;
 763     }
 764
 765     template<class M, class E1, class E2>
 766     BOOST_UBLAS_INLINE
 767     M &
 768     opb_prod (const matrix_expression<E1> &e1,
 769               const matrix_expression<E2> &e2,
 770               M &m,
 771               dense_proxy_tag, column_major_tag) {
 772         typedef M matrix_type;
 773         typedef const E1 expression1_type;
 774         typedef const E2 expression2_type;
 775         typedef typename M::size_type size_type;
 776         typedef typename M::value_type value_type;
 777
 778 #if BOOST_UBLAS_TYPE_CHECK
 779         matrix<value_type, column_major> cm (m);
 780         typedef typename type_traits<value_type>::real_type real_type;
 781         real_type merrorbound (norm_1 (m) + norm_1 (e1) * norm_1 (e2));
 782         indexing_matrix_assign<scalar_plus_assign> (cm, prod (e1, e2), column_major_tag ());
 783 #endif
 784         size_type size (BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()));
 785         for (size_type k = 0; k < size; ++ k) {
 786             vector<value_type> ce1 (column (e1 (), k));
 787             vector<value_type> re2 (row (e2 (), k));
 788             m.plus_assign (outer_prod (ce1, re2));
 789         }
 790 #if BOOST_UBLAS_TYPE_CHECK
 791         BOOST_UBLAS_CHECK (norm_1 (m - cm) <= 2 * std::numeric_limits<real_type>::epsilon () * merrorbound, internal_logic ());
 792 #endif
 793         return m;
 794     }
 795
 796     // Dispatcher
 797
 798   /** \brief computes <tt>M += A X</tt> or <tt>M = A X</tt> in an
 799           optimized fashion.
 800
 801           \param e1 the matrix expression \c A
 802           \param e2 the matrix expression \c X
 803           \param m  the result matrix \c M
 804           \param init a boolean parameter
 805
 806           <tt>opb_prod(A, X, M, init)</tt> implements the well known
 807           axpy-product. Setting \a init to \c true is equivalent to call
 808           <tt>M.clear()</tt> before <tt>opb_prod</tt>. Currently \a init
 809           defaults to \c true, but this may change in the future.
 810
 811           This function may give a speedup if \c A has less columns than
 812           rows, because the product is computed as a sum of outer
 813           products.
 814
 815           \ingroup blas3
 816
 817           \internal
 818
 819           template parameters:
 820           \param M type of the result matrix \c M
 821           \param E1 type of a matrix expression \c A
 822           \param E2 type of a matrix expression \c X
 823   */
 824     template<class M, class E1, class E2>
 825     BOOST_UBLAS_INLINE
 826     M &
 827     opb_prod (const matrix_expression<E1> &e1,
 828               const matrix_expression<E2> &e2,
 829               M &m, bool init = true) {
 830         typedef typename M::value_type value_type;
 831         typedef typename M::storage_category storage_category;
 832         typedef typename M::orientation_category orientation_category;
 833
 834         if (init)
 835             m.assign (zero_matrix<value_type> (e1 ().size1 (), e2 ().size2 ()));
 836         return opb_prod (e1, e2, m, storage_category (), orientation_category ());
 837     }
 838     template<class M, class E1, class E2>
 839     BOOST_UBLAS_INLINE
 840     M
 841     opb_prod (const matrix_expression<E1> &e1,
 842               const matrix_expression<E2> &e2) {
 843         typedef M matrix_type;
 844
 845         matrix_type m (e1 ().size1 (), e2 ().size2 ());
 846         return opb_prod (e1, e2, m, true);
 847     }
 848
 849 }}}
 850
 851 #endif