Boost logo

Ublas :

From: Shawn D. Pautz (sdpautz_at_[hidden])
Date: 2005-10-19 10:44:21


Good point. I had tried to distill the code down; perhaps I went too
far. Here's a less distilled version that I can't compile:

template <class Matrix>
typename Matrix::value_type determinant(const Matrix& A)
    assert(A.size1() == A.size2());

    typedef typename Matrix::value_type value_type;
    typedef typename Matrix::size_type size_type;

    Matrix LU = A; // mutable type?
    ublas::permutation_matrix<size_type> pivots(LU.size1());
    if (lu_factorize(LU, pivots) != 0) // this line won't compile in
some situations
    return 0;

    value_type det = 1.;
    for (size_type i = 0; i < LU.size1(); ++i)
    det *= LU(i,i);
    for (typename ublas::permutation_matrix<size_type>::const_iterator
         pivotIter = pivots.begin();
     pivotIter != pivots.end(); ++pivotIter)
    if (*pivotIter != pivotIter.index())
        det *= -1.;

    return det;

int main()
    ublas::matrix<double> A(2,2);
    ublas::matrix<double> B(2,2);
    determinant(A); // compiles
    determinant(B); // compiles
    determinant(ublas::prod(A,B)); // won't compile
    determinant(ublas::matrix<double>(ublas::prod(A,B))); // compiles
>(ublas::prod(A,B))); // compiles
    ublas::matrix<double> tempMatrix(ublas::prod(A,B)); // compiles
    determinant(tempMatrix); // compiles

I note that only one of the lines above will not compile now; thanks for
catching the const issue. But what do I need to do to get the most
compact line above to compile? Is the problem in my local template
code? Is there a problem in the ublas conversions? More generally,
should I be able to string together various matrix operations, not just
the basic mathematical ones like operator+ but algorithms with a
matrix-like return type? Or do I need to insert explicit conversions
from time to time like above?

Thanks for assisting me to understand this.


Karl Meerbergen wrote:

> Shawn D. Pautz wrote:
>> All,
>> Here's some more information that may be relevant to this
>> discussion. The code below does not compile with gcc 3.2.3 (I've
>> suppressed the include statements; I've tried including a lot of the
>> ublas header files without success). Which lines ought to work, and
>> which lines are an error on my part?
>> int main()
>> {
>> ublas::matrix<double> A(2,2);
>> ublas::matrix<double> B(2,2);
>> lu_factorize(A); // compiles
>> lu_factorize(B); // compiles
>> lu_factorize(ublas::prod(A,B)); // won't compile
>> lu_factorize(ublas::matrix<double>(ublas::prod(A,B))); // won't
>> compile
>> lu_factorize(static_cast<ublas::matrix<double>
>> >(ublas::prod(A,B))); // won't compile
>> ublas::matrix<double> tempMatrix(ublas::prod(A,B)); // compiles
>> lu_factorize(tempMatrix); // compiles
>> }
> I am not familiar with lu_factorize(), but I guess the reason why some
> lines do not compile is that you pass a const object to a function
> that requires a reference to a non-const object. This is not possible.
> Karl

Shawn D. Pautz
Radiation Transport Department, MS 1179
Sandia National Laboratories
P.O. Box 5800
Albuquerque, NM 87185
Phone: 505-284-4291
Fax: 505-844-0092
E-Mail: sdpautz_at_[hidden]