From: Gunter Winkler (guwi17_at_[hidden])
Date: 2007-05-30 18:01:43
Am Mittwoch, 30. Mai 2007 21:37 schrieb Albert Strasheim:
> This seems like a viable approach for small A, but I'm running into
> problems when dealing with large matrices. Code:
>
> #include <boost/numeric/ublas/matrix.hpp>
> #include <boost/numeric/ublas/matrix_proxy.hpp>
> #include <boost/numeric/ublas/vector.hpp>
.. may be
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
>
> using namespace boost::numeric::ublas;
>
> int main(int argc, char* argv[])
> {
> matrix<double> A(196608, 1024);
> vector<double> b(A.size2());
> scalar_vector<double> e(A.size1(), 1.0);
> #if 0
> noalias( A ) += outer_prod(e, b);
^^^^^^^ this prevents temporaries
> #else
> for (std::size_t i = 0; i < A.size1(); ++i)
> {
> noalias( row(A, i) ) += b;
^^^^^^^ this prevents temporaries
> }
> #endif
> return 0;
> }
> However, it seems the outer_prod doubles the memory requirements of
> this program compared to using a for-loop. Matrix A is 1536 MB, which
> can be allocated on my 32-bit Linux system without problems. However,
> I get a bad_alloc exception when trying to add the outer product.
The temporary is created because uBLAS assumes aliasing of rhs and lhs
by default (if rhs is not a plain matrix/vector). If you know for sure
that there is no aliasing, you should always use the noalias()
function.
> Indeed, for smaller sizes of A I can see the resident memory size of
> my program doubles while it does the outer product.
The only inefficiency that remains is the unnecessary computation of
1.0*b(i). This could be avoided by a new vector type 'one_vector'
similar to 'zero_vector' and a sufficiently smart compiler ...
mfg
Gunter