[petsc-users] LSQR with Jacobi preconditioning

[Mihai Alexe]

Matt,

Thank you. I'll try to implement some of the changes you suggested.

As for the preconditioner, I can think of two ways of going about it:

a. Define my "own" Jacobi preconditioner as in ksp/ex15.c & use
MatGetColumnNorms on A_mat to get the diagonal entries of A_mat^T * A_mat.
I don't very much like this option.

or

b. Leave P_mat as defined before and implement MATOP_GET_DIAGONAL in

PetscErrorCode myDiagFunc(Mat mat, Vec diag_elements)

But how can I access to A_mat from inside myDiagFunc? Through the shell
matrix context object? I couldn't find an example of how to do so.

Which option do you think would work best?

cheers,
Mihai


On Mon, Feb 25, 2013 at 12:32 PM, Matthew Knepley <knepley at gmail.com> wrote:

> On Mon, Feb 25, 2013 at 4:28 AM, Mihai Alexe <mihai at vt.edu> wrote:
>
>> Hello,
>>
>> I'm solving a small least-squares problem with LSQR. I want to set up
>> Jacobi preconditioning, so I do something like the following:
>>
>
> Comments on the code below:
>
> Always check all the return codes with CHKERRQ()
>
>
>> * MatCreateMPIAIJWithSplitArrays( PETSC_COMM_WORLD, *locrow, *loccol,
>> nrow,*
>> *  *ncol, onrowidx, oncolidx,*
>> *  (PetscScalar*) onvals, offrowidx, offcolidx,*
>> *  (PetscScalar*) values, &A_mat );*
>> **
>>
>
> Why would you do this instead of just using MatSetValues(). Its fragile,
> more complex, and cannot use other matrix types.
>
>
>> */* Create Linear Solver Context (and set options) */*
>> *  KSPCreate( PETSC_COMM_WORLD, &solksp );*
>> *  KSPSetType( solksp, "lsqr" );*
>>
>
> I would replace this with KSPSetFromOptions(solksp) so that you can change
> solvers from the command line using -ksp_type.
>
>
>> *  MatCreateNormal(A_mat, &P_mat);*
>>
>
> This is just a shell matrix to apply A^T A. In order to use Jacobi, you can
>
>   MatShellSetOperation(P_mat, MATOP_GET_DIAGONAL, myDiagFunc)
>
> where myDiagFunc() return the vector with a^T_i a_i for each slot.
>
>
>> *  MatSetUp(P_mat); *
>> *  KSPSetOperators( solksp, A_mat, P_mat, DIFFERENT_NONZERO_PATTERN );*
>> *
>> *
>> *  KSPGetPC( solksp, &pc);*
>> *  PCSetType( pc, PCJACOBI );*
>>
>
> Again, why not just use -pc_type jacobi.
>
>    Matt
>
>
>> *  KSPSetFromOptions( solksp );*
>> *  KSPSolve( solksp, b_vec, x_vec );*
>>
>> However, PETSc detects zeros on the diagonal of P_mat:
>>
>> [0] PCSetUp_Jacobi(): Zero detected in diagonal of matrix, using 1 at
>> those locations
>>
>> I then printed out the indices at which the zeros are detected. It turns
>> out that PETSC believes that all the diagonal elements are zero.
>> There is no zero column in my A_mat though... I confirmed this by
>> examining both A_mat and A_mat^T * A_mat in Matlab.
>>
>> What am I doing wrong?
>>
>> kind regards,
>> Mihai
>>
>
>
>
> --
> 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 --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130225/86b2b81c/attachment.html>