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