|
|
Ublas : |
From: pinacle (pinacle2000_at_[hidden])
Date: 2007-09-07 12:25:51
Hi Karl,
I understand UBLAS may be slower since it uses a lot c++ features. For
accessing the i-th (or i+1) element in an array A, A(i) for UBLAS vector
will be slower than A[i] for C array. But UBLAS is quite straight forward
for me to implement general calculations comparing with using arrays with
pointers. If the binding of UMFPACK could improve the performance to a good
extent, I'll give it a try.
BTW, I've solved my compile problem by changing the LINE 147 of lu.cpp from
project (mci, range (i + 1, size1)) *= value_type (1) / m (i, i);
to
project (mci, range (i + 1, size1)) *= value_type (1);
project (mci, range (i + 1, size1)) /= m (i, i);
Hope the developer could give a fix on this file. In fact, I feel the lu.cpp
is not well checked or maintained. I saw several discussions about it
before, but it still remains unchanged.
Cheers!
Feng
-----Original Message-----
From: ublas-bounces_at_[hidden] [mailto:ublas-bounces_at_[hidden]]
On Behalf Of Karl Meerbergen
Sent: Thursday, September 06, 2007 11:33 PM
To: ublas mailing list
Subject: Re: [ublas] lu factorize cannot be compiled on sparse complex
matrix?
Dear Feng Gao,
It surprises me that many people are trying to use ublas routines for
solving linear systems, while excellent external software is available:
LAPACK for dense matrices, UMFPACK, MUMPS for sparse.
In your example, you are trying to solve a 5x5 system of a sparse matrix. I
guess this is just a small test example to make sure the code works. Suppose
you want to solve a 20,000x20,0000 sparse system. Which performance do you
expect to obtain with lu_factorize? For a sparse matrix, this is, in
general, going to be a disaster.
I suggest to spend some time to develop a link with external software rather
than (sorry for being so rude) wasting your time using lu_factorize. There
are C++ bindings to UMFPACK in boost.sandbox. I am currently working on
MUMPS bindings, which I hope to make available soon. You will have to spend
some time installing packages and compiling and linking, maybe a day or half
a day. It may be a frustrating and hard job, but once you have done this,
your code will have a good performance.
Have a nice day,
Karl
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_t
> ype 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
>