Re: [ublas] find2 required by Matrix?

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

Karl,

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
    determinant(static_cast<ublas::matrix<double> 
 >(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.

Shawn

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: [phone]
Fax: [phone]
E-Mail: sdpautz_at_[hidden]