[petsc-users] LSQR with Jacobi preconditioning

[Matthew Knepley]

On Mon, Feb 25, 2013 at 7:41 AM, Mihai Alexe <mihai at vt.edu> wrote:

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

Yes, put it in the context.

struct MyStruct {
  Mat A;
};

struct myStruct ctx;

ctx.A = A_mat;
MatShellSetContext(P_mat, &ctx);

Then use MatShellGetContext() in myDiagFunc.

   Matt


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


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