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 |
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 |
Free forum by Nabble | Edit this page |