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>
