[petsc-users] In-place MUMPS call?

Jack Poulson jack.poulson at gmail.com
Sun Mar 6 22:54:02 CST 2011


On Sun, Mar 6, 2011 at 10:50 PM, Jed Brown <jed at 59a2.org> wrote:

> On Sun, Mar 6, 2011 at 20:36, Jack Poulson <jack.poulson at gmail.com> wrote:
>
>> The documentation seems to suggest that I can choose a sparse-direct
>> solver and factor a MATMPISBAIJ matrix A out-of-place like this:
>>
>>     Mat F;
>>     MatGetFactor( A, MAT_SOLVER_MUMPS, MAT_FACTOR_CHOLESKY, &F );
>>     MatFactorInfo cholInfo;
>>     cholInfo.fill = 3.0; // arbitrary choice
>>     cholInfo.dtcol = 0.5; // arbitrary choice
>>     IS perm;
>>     ISCreateStride( comm, size, 0, 1, &perm );
>>     MatCholeskyFactorSymbolic( F, A, perm, &cholInfo );
>>     MatCholeskyFactorNumeric( F, A, &cholInfo );
>>
>> but the in-place alternative, MatCholeskyFactor
>>
>
> 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:
>
> PetscErrorCode MatCholeskyFactor_SeqSBAIJ(Mat A,IS perm,const MatFactorInfo
> *info)
> {
>   PetscErrorCode ierr;
>   Mat            C;
>
>   PetscFunctionBegin;
>   ierr = MatGetFactor(A,"petsc",MAT_FACTOR_CHOLESKY,&C);CHKERRQ(ierr);
>   ierr = MatCholeskyFactorSymbolic(C,A,perm,info);CHKERRQ(ierr);
>   ierr = MatCholeskyFactorNumeric(C,A,info);CHKERRQ(ierr);
>   A->ops->solve            = C->ops->solve;
>   A->ops->solvetranspose   = C->ops->solvetranspose;
>   ierr = MatHeaderMerge(A,C);CHKERRQ(ierr);
>   PetscFunctionReturn(0);
> }
>
> Is there a way to tell MatCholeskyFactor which approach to use (i.e.,
>> MAT_SOLVER_MUMPS),
>>
>
> No, maybe this should  be added.
>
>
>> or, failing that, can I destroy my original matrix 'A' and keep the factor
>> 'F' around for solves?
>>
>
> Do this, there is no performance penalty.
>

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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110306/de4123a8/attachment.htm>


More information about the petsc-users mailing list