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.matMaybe 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