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>
