If you look at the eigenvalues of the matrix closely you will see it has two parts. Eigenvalues with real part between .0033 and .01 and eigenvalues between 2 and 12. I suspect this two levels of well-conditioned pieces is due to the parts of the matrix have zeros on the diagonals and the rest without zeros (this is why I asked if the problem was a Stokes problem). Perhaps a Schur complement type preconditioner done with the PCFieldSplit preconditioner is the way to go, you would eliminate in the Schur complement the rows/columns of the matrix that do not have the zeros on the diagonal.
Barry On Aug 4, 2010, at 10:10 AM, Hong Zhang wrote: > Filippo, > Sorry, your matrix is ill-conditioned, gmres+ilu does not converge in > true residual > (not preconditioned residual!): > ./ex10 -f0 $D/A_rhs_spiga -ksp_monitor_true_residual > 0 KSP preconditioned resid norm 1.011174511507e+01 true resid norm > 2.084628109126e-02 ||r(i)||/||b|| 1.000000000000e+00 > 1 KSP preconditioned resid norm 1.498441679989e+00 true resid norm > 4.919980218651e-02 ||r(i)||/||b|| 2.360123706052e+00 > 2 KSP preconditioned resid norm 1.657440211742e-01 true resid norm > 1.624190902280e-02 ||r(i)||/||b|| 7.791274113447e-01 > 3 KSP preconditioned resid norm 7.568788163764e-02 true resid norm > 1.516900048065e-02 ||r(i)||/||b|| 7.276597880572e-01 > 4 KSP preconditioned resid norm 3.158616884464e-02 true resid norm > 1.703303336172e-02 ||r(i)||/||b|| 8.170777937395e-01 > 5 KSP preconditioned resid norm 6.169977289081e-08 true resid norm > 1.629219484085e-02 ||r(i)||/||b|| 7.815396314348e-01 > > Using superlu, I get condition number > ./ex10 -f0 $D/A_rhs_spiga -pc_type lu -pc_factor_mat_solver_package > superlu -mat_superlu_conditionnumber > > Recip. condition number = 1.171784e-04 > > For such tiny matrix (27x27), condition number=1.e+4. > Using superlu_dist or mumps LU as precondition might be the only > option for your application. > > Hong > > > > On Wed, Aug 4, 2010 at 9:53 AM, Hong Zhang <hzhang at mcs.anl.gov> wrote: >> 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 >>> >>> >>
