On Sun, Mar 6, 2011 at 10:50 PM, Jed Brown <span dir="ltr">&lt;<a href="mailto:jed@59a2.org">jed@59a2.org</a>&gt;</span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div class="gmail_quote"><div class="im">On Sun, Mar 6, 2011 at 20:36, Jack Poulson <span dir="ltr">&lt;<a href="mailto:jack.poulson@gmail.com" target="_blank">jack.poulson@gmail.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

The documentation seems to suggest that I can choose a sparse-direct solver and factor a MATMPISBAIJ matrix A out-of-place like this:<div><br></div><div><div>    Mat F;</div><div>    MatGetFactor( A, MAT_SOLVER_MUMPS, MAT_FACTOR_CHOLESKY, &amp;F );</div>


<div>    MatFactorInfo cholInfo;</div><div>    cholInfo.fill = 3.0; // arbitrary choice</div><div>    cholInfo.dtcol = 0.5; // arbitrary choice</div><div>    IS perm;</div><div>    ISCreateStride( comm, size, 0, 1, &amp;perm );</div>


<div>    MatCholeskyFactorSymbolic( F, A, perm, &amp;cholInfo );</div><div>    MatCholeskyFactorNumeric( F, A, &amp;cholInfo );</div></div><div><br></div><div>but the in-place alternative, MatCholeskyFactor</div></blockquote>

<div><br></div></div><div>I don&#39;t know of any sparse direct solver that actually does in-place factorization. An example of what MatCholeskyFactor does to make itself look in-place:</div><div><br></div><div><div>PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat A,IS perm,const MatFactorInfo *info)</div>

<div>{</div><div>  PetscErrorCode ierr;</div><div>  Mat            C;</div><div><br></div><div>  PetscFunctionBegin;</div><div>  ierr = MatGetFactor(A,&quot;petsc&quot;,MAT_FACTOR_CHOLESKY,&amp;C);CHKERRQ(ierr);</div><div>

  ierr = MatCholeskyFactorSymbolic(C,A,perm,info);CHKERRQ(ierr);</div><div>  ierr = MatCholeskyFactorNumeric(C,A,info);CHKERRQ(ierr);</div><div>  A-&gt;ops-&gt;solve            = C-&gt;ops-&gt;solve;</div><div>  A-&gt;ops-&gt;solvetranspose   = C-&gt;ops-&gt;solvetranspose;</div>

<div>  ierr = MatHeaderMerge(A,C);CHKERRQ(ierr);</div><div>  PetscFunctionReturn(0);</div><div>}</div></div><div class="im"><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

<div>Is there a way to tell MatCholeskyFactor which approach to use (i.e., MAT_SOLVER_MUMPS), </div></blockquote><div><br></div></div><div>No, maybe this should  be added.</div><div class="im"><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

<div>or, failing that, can I destroy my original matrix &#39;A&#39; and keep the factor &#39;F&#39; around for solves?<br></div></blockquote><div><br></div></div><div>Do this, there is no performance penalty.</div></div>

</blockquote></div><br><div>Thanks Jed. I was aware that the factorization was not actually in-place, but the fact that the symbolic and numeric factorizations required the original matrix to be around made me unsure of when I was allowed to free it. This interface should be perfectly sufficient. </div>