[ublas] lu factorize cannot be compiled on sparse complex matrix?

From: pinacle (pinacle2000_at_[hidden])
Date: 2007-09-06 19:13:38

Dear all,

I've been using UBLAS for solving problems like Ax=b, where A is a complex
dense matrix, and x and b are complex dense vectors. Today, I tried to
change A to a sparse matrix, but I got compiling errors as:
Error	10	error C2784: 'matrix_binary_scalar2_traits<E1,const
T2,boost::numeric::ublas::scalar_divides<E1::value_type,T2>>::result_type
boost::numeric::ublas::operator /(const
boost::numeric::ublas::matrix_expression<E> &,const T2 &)' : could not
deduce template argument for 'const
boost::numeric::ublas::matrix_expression<E> &' from 'std::complex<double>'
D:\boost\boost\numeric\ublas\lu.hpp	147

But if I use real sparse A, no error appears. So the conclusion now is that
lu_factorize works for real matrices no matter dense or sparse, and works
for complex dense matrix, but doesn't work for complex sparse matrix.

I do need a lu_factorize (and lu_substitute of course) working for complex
sparse matrices. How could I solve the problem? BTW, I am using boost
1.34.0.

Here is a sample code for testing:

#include <complex>
#include <iostream>

#include <boost/numeric/ublas/storage_sparse.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>

using namespace boost::numeric::ublas;
using namespace std;

int main(int argc, char *argv[])
{
    const size_t size = 5;

    // error:
    //compressed_matrix<complex<double> > A(size, size);
    //compressed_matrix<complex<double> > ALU(size, size);
    // good:
    compressed_matrix<double> A(size, size);
    compressed_matrix<double> ALU(size, size);

    permutation_matrix<> pm(size);

    for (size_t i=0; i<A.size1(); ++i)	{
        if (i>=1) { 
            A.push_back(i,i-1,-2.0); 
        }
        A.push_back(i,i,4);
        if (i+1<A.size2()) { 
            A.push_back(i,i+1,-1.0); 
        }
    }
    ALU = A;
    lu_factorize(ALU,pm); // LU factorized

    cout << "A = " << A << endl;
    cout << "ALU = " << ALU << endl;

    return EXIT_SUCCESS;
};

Thanks,
Feng Gao