Boost logo

Ublas :

From: Gunter Winkler (guwi17_at_[hidden])
Date: 2005-01-13 10:16:20


-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On Thursday 13 January 2005 14:39, Rakesh K Sinha wrote:
> I had posted this to the Boost Archives list sometime before. Since I did
> not get any response there, I thought this list might be more appopriate
> and hence am posting it here, again.
>
> For a particular application of mine, I ned to solve the equation Ax =
> B,
> i.e. x = inv(A) * B

I would generally avoid the computation of any inverse. uBLAS currently does
not provide and probably will not provide such a function.

The way to solve linear equations is to first compute the LU factorization of
the matrix A and then solve to triangular systems. The code looks like

const int n=10;
permutation_matrix<double> P(n);
matrix<double> A(n,n);
vector<double> x(n);
vector<double> rhs(n);

// fill matrix and rhs

lu_factorize(A,P);
x = rhs;
lu_substitute(A,P,rhs);

> This fits my need perfectly.

in order to solve triangular systems you can use the solve() funtion. You can
find examples in lu.hpp. If the documentation is missing or fuzzy, please
send me a mail.

> * Are there any other caveats with this method like it works only for
> dense matrices. How does it scale to sparse matrices ?

Yes. lu_factorize() works only for dense matrices reliably. The trianglular
solves should work for any matrix type, but may fail for adapted matrices
like trianglur_adaptor<compressed_matrix<...>>. We would be very pleased if
you provide a short (regression-)test program if there are any problems.

> * Also what does that lower_tag() essentially imply ? What is its
> mathematical significance ?

This should be explained in the documentation of matrix expressions
http://cvs.sourceforge.net/viewcvs.py/*checkout*/boost/boost/libs/numeric/ublas/doc/matrix_expression.htm

There are lower (i<=j), stric_lower (i<j) and unit_lower (i<j + unit diagonal)

> * Is there any function to calculate the singularity of a square matrix

No. But you can do a LU decompositon, if it fails the determinant is zero,
otherwise the product of the diagonal elements of the LU-Matrix give the
absolute value of the determinant. The sign can be extracted from the
characteristic of the permutation. More details should be in any good book
about numerics.

short: P * A = L * U; det(A) = prod(U_ii,i=1..n) * sign(P);

(For solving sparse linear systems I prefer iterative methods.)

mfg
Gunter
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.4 (GNU/Linux)

iD8DBQFB5pDKghsXb/J0PcERAvnBAKDHAo4eGeB0JP2r+/84reC8qsb1BwCfYALY
yHPObvbxNrpaOXfhCA/9tfo=
=sW45
-----END PGP SIGNATURE-----