question about complex sparse linear algebra

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
c.
Reply | Threaded
Open this post in threaded view
|

question about complex sparse linear algebra

c.
Hi,

I have a question about solving a sparse linear system with
a complex system matrix.

The system I am solving is of the form

A = (J + 1i * w * M)
A * x = b

where J and M are sparse real matrices, b is a real vector and
w is a (small) real number.

J is nonsingular, M is singular, A can be proven to be invertible
(in exact arithmetics) for any value of w.

What I am experiencing is the following: if I use mldivide I get

>> A = (J + 1i * w * M);
>> x1 = A \ res;
warning: matrix singular to machine precision, rcond = 1.49212e-27

and the contents of x1 are just garbage while, if I perform LU first, I get:

>> [L, U, P, Q, R] = lu (A);
>> x2 = Q * (U \ (L \ (P * (R \ res))));

and x2 is the physically expected result.

Can anyone comment about why the the two approaches differ?
What are the underlying library methods invoked in the two cases?

Thanks,
c.


P.S. In case anyone would like to try the above examples, the matirces are available here:
http://www1.mate.polimi.it/~carlo/out.mat




_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: question about complex sparse linear algebra

siko1056
c. wrote

> Hi,
>
> I have a question about solving a sparse linear system with
> a complex system matrix.
>
> The system I am solving is of the form
>
> A = (J + 1i * w * M)
> A * x = b
>
> where J and M are sparse real matrices, b is a real vector and
> w is a (small) real number.
>
> J is nonsingular, M is singular, A can be proven to be invertible
> (in exact arithmetics) for any value of w.
>
> What I am experiencing is the following: if I use mldivide I get
>
>>> A = (J + 1i * w * M);
>>> x1 = A \ res;
> warning: matrix singular to machine precision, rcond = 1.49212e-27
>
> and the contents of x1 are just garbage while, if I perform LU first, I
> get:
>
>>> [L, U, P, Q, R] = lu (A);
>>> x2 = Q * (U \ (L \ (P * (R \ res))));
>
> and x2 is the physically expected result.
>
> Can anyone comment about why the the two approaches differ?
> What are the underlying library methods invoked in the two cases?
>
> Thanks,
> c.
>
>
> P.S. In case anyone would like to try the above examples, the matirces are
> available here:
> http://www1.mate.polimi.it/~carlo/out.mat

Maybe you already found an answer, my conclusion is, that both algorithms
use UMFPACK of SuiteSparse and that \ uses the numeric factorization [1] and
lu the symbolic factorization [2] approach.  To see where "\" and "lu" are
handled, I started a debugger (as explained in the wiki [3] + kdbg as
graphical frontend) and stepped though the code for your example for "w =
1":

For "\" I found the following trace:

1. libinterp/operators/op-scm-m.cc:

DEFBINOP (ldiv, sparse_complex_matrix, matrix)

2. liboctave/array/CSparse.cc:

ComplexMatrix
SparseComplexMatrix::solve (MatrixType& mattype, const Matrix& b,
                            octave_idx_type& err, double& rcond,
                            solve_singularity_handler sing_handler,
                            bool singular_fallback) const

ComplexMatrix
SparseComplexMatrix::fsolve (MatrixType& mattype, const Matrix& b,
                             octave_idx_type& err, double& rcond,
                             solve_singularity_handler sing_handler,
                             bool calc_cond) const

// Here UMFPACK detects the bad condition number doing a numeric
factorization.

void *
SparseComplexMatrix::factorize (octave_idx_type& err, double& rcond,
                                Matrix& Control, Matrix& Info,
                                solve_singularity_handler sing_handler,
                                bool calc_cond) const

// UMFPACK_ZNAME (numeric)


For "lu" the way goes via

1. liboctave/numeric/lu.cc:

2. liboctave/numeric/sparse-lu.cc: (directed to UMFPACK qsymbolic
factorization)

template <typename lu_type>
sparse_lu<lu_type>::sparse_lu (const lu_type& a, const Matrix& piv_thres,
bool scale)

// UMFPACK_ZNAME (qsymbolic)

Maybe as a start for further investigation,

Kai

[1]
https://github.com/sergiud/SuiteSparse/blob/master/UMFPACK/Include/umfpack_numeric.h
[2]
https://github.com/sergiud/SuiteSparse/blob/master/UMFPACK/Include/umfpack_qsymbolic.h
[3] http://wiki.octave.org/Debugging_Octave



--
Sent from: http://octave.1599824.n4.nabble.com/Octave-General-f1599825.html

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave