<div dir="ltr">On Mon, Feb 25, 2013 at 7:41 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">
Matt,<div><br></div><div>Thank you. I'll try to implement some of the changes you suggested.</div><div><br></div><div>As for the preconditioner, I can think of two ways of going about it:</div><div><br></div><div>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.</div>
<div><br></div><div>or</div><div><br></div><div>b. Leave P_mat as defined before and implement MATOP_GET_DIAGONAL in</div><div><br></div><div>PetscErrorCode myDiagFunc(Mat mat, Vec diag_elements)</div><div><br></div><div>
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.</div></blockquote><div><br></div><div style>Yes, put it in the context.</div>
<div style><br></div><div style>struct MyStruct {</div><div style> Mat A;</div><div style>};</div><div style><br></div><div style>struct myStruct ctx;</div><div style><br></div><div style>ctx.A = A_mat;</div><div style>MatShellSetContext(P_mat, &ctx);</div>
<div style><br></div><div style>Then use MatShellGetContext() in myDiagFunc.</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>Which option do you think would work best?</div><div>
<br></div><div>cheers,</div><div>Mihai</div><div><br><br><div class="gmail_quote">On Mon, Feb 25, 2013 at 12:32 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>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><div class="gmail_extra"><div class="gmail_quote"><div><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><div>Comments on the code below:</div>
<div><br></div><div>Always check all the return codes with CHKERRQ()</div><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><div>Why would you do this instead of just using MatSetValues(). Its fragile, more complex, and cannot use other matrix types.</div>
<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><div>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>This is just a shell matrix to apply A^T A. In order to use Jacobi, you can</div>
<div><br></div><div> MatShellSetOperation(P_mat, MATOP_GET_DIAGONAL, myDiagFunc)</div><div><br></div><div>where myDiagFunc() return the vector with a^T_i a_i for each slot.</div><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><div>Again, why not just use -pc_type jacobi.</div><div><br></div>
<div> Matt</div><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></div><span><font color="#888888"><br><br clear="all"><span class="HOEnZb"><font color="#888888"><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
</font></span></font></span></div></div>
</blockquote></div><br></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>