Filippo, Yes, the matrix is well conditioned. Using ~petsc-dev/src/mat/examples/tests/ex78.c, I converted your matlab matrix and rhs into petsc binary format. Then run it with petsc-dev/src/ksp/ksp/examples/tutorials/ex10.c:
./ex10 -f0 $D/A_rhs_spiga -pc_type lu which gives same error on numerical factorization as yours. Then switching to superlu and mumps LU direct sovler, I get ./ex10 -f0 $D/A_rhs_spiga -pc_type lu -pc_factor_mat_solver_package superlu Number of iterations = 1 Residual norm < 1.e-12 Using petsc default sovler (gmres+ilu): ./ex10 -f0 $D/A_rhs_spiga -ksp_monitor 0 KSP Residual norm 1.011174511507e+01 1 KSP Residual norm 1.498441679989e+00 2 KSP Residual norm 1.657440211742e-01 3 KSP Residual norm 7.568788163764e-02 4 KSP Residual norm 3.158616884464e-02 5 KSP Residual norm 6.169977289081e-08 Number of iterations = 5 Residual norm 0.0162922 As you see, it converges well. For you info., I attached the binary file A_rhs_spiga here. Hong On Tue, Aug 3, 2010 at 11:09 PM, Barry Smith <bsmith at mcs.anl.gov> wrote: > > ??Matt and Hong, > ??The eigenvalues are all nonzero, so I don't see how the matrix could be > singular, all the real parts are positive so it is even a nice matrix. But a > bunch of diagonal entries are identically zero, this screws up most LU > factorization schemes that do not do pivoting. (Matlab does pivoting). > > 0.0109 + 0.0453i > 0.0109 - 0.0453i > 0.0071 + 0.0131i > 0.0071 - 0.0131i > 0.0072 + 0.0121i > 0.0072 - 0.0121i > 0.0049 + 0.0062i > 0.0049 - 0.0062i > 0.0059 + 0.0012i > 0.0059 - 0.0012i > 0.0034 + 0.0031i > 0.0034 - 0.0031i > 0.0068 + 0.0226i > 0.0068 - 0.0226i > 0.0061 + 0.0061i > 0.0061 - 0.0061i > 0.0033 + 0.0028i > 0.0033 - 0.0028i > 6.0000 > 8.0000 > 2.0000 > 8.0000 > 12.0000 > 4.0000 > 2.0000 > 4.0000 > 2.0000 > > On Aug 3, 2010, at 11:03 PM, Matthew Knepley wrote: > > On Tue, Aug 3, 2010 at 10:54 PM, Filippo Spiga > <filippo.spiga at disco.unimib.it> wrote: >> >> ?Dear Hong, >>> >>> This confirms that your Jacobian is singular, thus none of linear >>> solvers would work. >> >> So do any preconditioner not help me to solve the problem? > > There can exist no solutions when the matrix is singular, thus you have a > problem > with either: > ??a) the problem definition > ??b) the creation of your matrix in PETSc > >> >> I put some stuff here: http://tinyurl.com/fil-petsc >> - "A_LS.m" is matrix (saved by PETSc) >> - "b_LS-m" >> - the file "eigenvalues_A" contains the eigenvalues of the matrix A, >> computed by MATLAB. >> >> I used "-pc_type lu" and 1 only processor. The result is the same of my >> previous email (*). > > This shows that your matrix is singular. > >> >> Anyway if I solve the problem using MATLAB I get the right solution. The >> formulation seems correct. To be > > What does this mean? What method in MATLAB? Some methods (like CG) can > iterate on rank deficient > matrices with a compatible rhs and get a solution, but other Krylov methods > will fail. Most preconditioners > will fail. > >> >> honest, the eigenvalues don't say me nothing. But I'm a computer >> scientist, not a mathematician. I'm not able to recognize which >> preconditioner I should use or which modifications (scaling all/part of the >> rows? reformulate the system in another way?...) do to solve the problem. >> From my side, it is not possible to try all the preconditioners and also it >> is not the right way... > > Actually, I strongly disagree. Preconditioners are very problem specific, > and it is often impossible > to prove which one will work for a certain problem. There are many > well-known results along these > lines, such as the paper of Greenbaum, Strakos, and Ptak on GMRES. > Experimentation is essential. > ?? ?Matt > >> >> Once again, thanks. >> >> (*) >> [0|23:14:58]: unknown: MatLUFactorNumeric_SeqAIJ() line 668 in >> src/mat/impls/aij/seq/aijfact.c: Zero pivot row 1 value 0 tolerance >> 2.77778e-14 * rowsum 0.0277778 >> >> -- >> >> Filippo SPIGA >> >> ?Nobody will drive us out of Cantor's paradise.? >> ? ? -- David Hilbert >> > > > > -- > What most experimenters take for granted before they begin their experiments > is infinitely more interesting than any results to which their experiments > lead. > -- Norbert Wiener > > -------------- next part -------------- A non-text attachment was scrubbed... Name: A_rhs_spiga Type: application/octet-stream Size: 3024 bytes Desc: not available URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20100804/8f3b60c7/attachment.obj>
