<div dir="ltr">On Mon, Feb 25, 2013 at 4:28 AM, Mihai Alexe <span dir="ltr"><<a href="mailto:mihai@vt.edu" target="_blank">mihai@vt.edu</a>></span> wrote:<br><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>Hello,</div><div><br></div><div>I'm solving a small least-squares problem with LSQR. I want to set up Jacobi preconditioning, so I do something like the following:  </div></blockquote><div><br></div><div style>Comments on the code below:</div>
<div style><br></div><div style>Always check all the return codes with CHKERRQ()</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><b> MatCreateMPIAIJWithSplitArrays( PETSC_COMM_WORLD, *locrow, *loccol, nrow,</b></div>

<div><b><span style="white-space:pre-wrap">                         </span>  *ncol, onrowidx, oncolidx,</b></div><div><b><span style="white-space:pre-wrap">                          </span>  (PetscScalar*) onvals, offrowidx, offcolidx,</b></div>
<div><b><span style="white-space:pre-wrap">                         </span>  (PetscScalar*) values, &A_mat );</b></div><div><b></b></div></blockquote><div><br></div><div style>Why would you do this instead of just using MatSetValues(). Its fragile, more complex, and cannot use other matrix types.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><b>/* Create Linear Solver Context (and set options) */</b></div><div><b>  KSPCreate( PETSC_COMM_WORLD, &solksp );</b></div>

<div><b>  KSPSetType( solksp, "lsqr" );</b></div></blockquote><div><br></div><div style>I would replace this with KSPSetFromOptions(solksp) so that you can change solvers from the command line using -ksp_type.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><b>  MatCreateNormal(A_mat, &P_mat);</b></div></blockquote><div><br></div><div style>This is just a shell matrix to apply A^T A. In order to use Jacobi, you can</div>
<div style><br></div><div style>  MatShellSetOperation(P_mat, MATOP_GET_DIAGONAL, myDiagFunc)</div><div style><br></div><div style>where myDiagFunc() return the vector with a^T_i a_i for each slot.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div><b>  MatSetUp(P_mat); </b></div><div><b>  KSPSetOperators( solksp, A_mat, P_mat, DIFFERENT_NONZERO_PATTERN );</b></div>
<div><b><br></b></div><div><b>  KSPGetPC( solksp, &pc);</b></div><div><b>  PCSetType( pc, PCJACOBI );</b></div></blockquote><div><br></div><div style>Again, why not just use -pc_type jacobi.</div><div style><br></div>
<div style>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><b>  KSPSetFromOptions( solksp );</b></div><div><b>  KSPSolve( solksp, b_vec, x_vec );</b></div>

<div><br></div><div>However, PETSc detects zeros on the diagonal of P_mat:</div><div><br></div><div><div>[0] PCSetUp_Jacobi(): Zero detected in diagonal of matrix, using 1 at those locations</div></div><div><br></div><div>

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

<div><br></div><div>What am I doing wrong?</div><div><br></div><div>kind regards,</div><div>Mihai</div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>