[petsc-users] LSQR with Jacobi preconditioning

[Matthew Knepley]

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/d2ee77f4/attachment.html>