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

[Hong Zhang]

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>