|
|
Ublas : |
From: Riccardo Rossi (rrossi_at_[hidden])
Date: 2007-09-07 03:13:03
Hi,
Implementing a direct solver using a sparse matrix storage it is NOT a
trivial issue.
There exist a few libraries doing it and ublas has bindings to them.
Examples i am aware of (for which there are bindings) are
SuperLU (which i am using) and UMFpack. Similar libraries are PARDISO
and MUMPS
one advantage of superLU over the other is that it allows to use
compressed row storage while (as i understand) UMFpack allows only
compressed column storage
greetings
Riccardo
pinacle wrote:
> 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
>
>
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
>
>