[petsc-users] Problems with preconditioners, which one?

[Hong Zhang]

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