On Sun, Mar 6, 2011 at 10:50 PM, Jed Brown <span dir="ltr"><<a href="mailto:jed@59a2.org">jed@59a2.org</a>></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"><<a href="mailto:jack.poulson@gmail.com" target="_blank">jack.poulson@gmail.com</a>></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, &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, &perm );</div>
<div> MatCholeskyFactorSymbolic( F, A, perm, &cholInfo );</div><div> MatCholeskyFactorNumeric( F, A, &cholInfo );</div></div><div><br></div><div>but the in-place alternative, MatCholeskyFactor</div></blockquote>
<div><br></div></div><div>I don'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,"petsc",MAT_FACTOR_CHOLESKY,&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->ops->solve = C->ops->solve;</div><div> A->ops->solvetranspose = C->ops->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 'A' and keep the factor 'F' 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>