PCGetFactoredMatrix

[Hong Zhang]

John,

You may look at ~petsc/src/mat/examples/tests/ex74.c
in which we use
|| y - x || as an indicator for
|| A - ICC ||.
In ex74.c, x is a randomly generated vector, b=A*x,
and ICC*y = b.
If you uncomment line 318
  printf("lf: %d, error: %G\n", lf,norm2);
and run ex74, you get
lf: -1, error: 3.33036E-15
lf: 0, error: 4.44135
lf: 1, error: 4.40183
lf: 2, error: 3.13597
lf: 3, error: 2.39443
lf: 4, error: 1.79942
lf: 5, error: 1.4183
lf: 6, error: 1.11197
lf: 7, error: 0.877789
lf: 8, error: 0.750784
lf: 9, error: 0.571567

which shows the error || y - x || for ICC(lf), lf=level of fill.

Hong

On Thu, 29 Nov 2007, Barry Smith wrote:

>
>  John,
>
>    There is no immediate way to do this.
> For the SeqAIJ format, we store both the LU in a single CSR format.
> with for each row first the part of L (below the diagonal) then 1/D_i
> then the part of U for that row. You can see how the triangular solves
> are done by looking at src/mat/impls/aij/seq/aijfact.c the routine 
> MatSolve_SeqAIJ()
> Note that it is actually more complicated due to the row and column 
> permutations
> (the factored matrix is stored in the ordering of the permutations).
> For BAIJ matrix the storage is similar except it is stored by block instead 
> of point
> and the inverse of the block diagonal is stored.
>
> One could take the MatSolve_SeqAIJ() routine and modify it to do the matrix
> vector product without too much difficulty.
>
> If you decide to do this we would gladly include it in our distribution.
>
>  Barry
>
> One can ask why we don't provide this functionality in PETSc since computing
> A - LU is a reasonable thing to do if one wants to understand the convergence
> of the method. The answer is two-fold 1) time and energy and 2) though we
> like everyone to use PETSc we driven more by people who are not interested
> in the solution algorithms etc but only in getting the answer easily and 
> relatively
> efficiently.
>
>
> On Nov 29, 2007, at 12:07 PM, John R. Wicks wrote:
>
>> I would like to compute the residual A - LU, where LU is the ILU
>> factorization of A.  What is the most convenient way of doing so?
>> 
>>> -----Original Message-----
>>> From: owner-petsc-users at mcs.anl.gov
>>> [mailto:owner-petsc-users at mcs.anl.gov] On Behalf Of Matthew Knepley
>>> Sent: Thursday, November 29, 2007 12:04 PM
>>> To: petsc-users at mcs.anl.gov
>>> Subject: Re: PCGetFactoredMatrix
>>> 
>>> 
>>> It depends on the package, but the petsc stuff stores L and U
>>> in one matrix.
>>>
>>>  Matt
>>> 
>>> On Nov 29, 2007 9:03 AM, John R. Wicks <jwicks at cs.brown.edu> wrote:
>>>> The documentation for PCGetFactoredMatrix is not clear.  What does
>>>> this return for ILU(0), for example? Does it return the
>>> product LU or
>>>> the in place factorization?
>>>> 
>>>> 
>>> 
>>> 
>>> 
>>> -- 
>>> 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
>>> 
>> 
>