[petsc-users] Determinant of the matrix

[Hong Zhang]

Jelena :
> Thank you very much, but I need to access this value in my program and not
> just print it, is there a way to do this? Also, I need to use Cholesky
> factorization, is there a way to avoid MUMPS so I could access the diagonal?

Google search 'determinant, mumps' leads to an email that I sent to an user:
http://lists.mcs.anl.gov/pipermail/petsc-users/2011-September/010225.html

I added mumps options ICNTL(30 - 33) to petsc-3.2.
With the latest petsc-3.2, you can get determinant with '-mat_mumps_icntl_33 1'
the result can be found using '-ksp_view'. For example,
petsc-3.2/src/ksp/ksp/examples/tutorials>./ex2 -pc_type lu
-pc_factor_mat_solver_package mumps -mat_mumps_icntl_33 1 -ksp_view
...
(RINFOG(12) RINFOG(13))*2^INFOG(34) (determinant): (0.900768,0)*(2^99)

Let us know if you need additional support - I guess you might want
MatGetMumpsRINFOG()?

Hong

>
>
> On Mon, Feb 25, 2013 at 7:13 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>
>>
>>   Jelena,
>>
>>      MUMPS has its own (complicated) data structure for storing the
>> factored matrix, you cannot easily look inside it.
>>
>>       However, fortunately MUMPS has a way of accessing the determinate of
>> the matrix directly. Look at src/mat/impls/aij/mpi/mumps.c and search for
>> the word determinant, you will find the command line argument to have it
>> computed and where it will be printed out if you call MatView() on the
>> factored matrix. You mush pass the PetscViewerFormat PETSC_VIEWER_ASCII_INFO
>> to PetscViewerSetFormat() before viewing the matrix.
>>
>>    Good luck,
>>
>>    Barry
>>
>> On Feb 25, 2013, at 6:02 PM, Jelena Slivka <slivkaje at gmail.com> wrote:
>>
>> > Hello,
>> > I am trying to find the determinant of the matrix. I have found this
>> > thread:
>> >
>> > http://lists.mcs.anl.gov/pipermail/petsc-users/2010-July/006716.html
>> >
>> > But, could you please explain how to access the data structure directly?
>> > So far I have tried the following:
>> >
>> > #include <petsc-private/matimpl.h>
>> > #include <../src/mat/impls/aij/seq/aij.h>
>> > ...
>> > // (Solving the system A*x = rhs
>> > MatSetOption(A, MAT_SPD, PETSC_TRUE);
>> > Mat F;
>> > MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F);
>> > MatFactorInfo factinfo;
>> > MatFactorInfoInitialize(&factinfo);
>> > IS perm, iscol;
>> > MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iscol);
>> > MatCholeskyFactorSymbolic(F, A, perm, &factinfo);
>> > MatCholeskyFactorNumeric(F, A, &factinfo);
>> > ISDestroy(&iscol);
>> > ISDestroy(&perm);
>> >
>> > MatSolve(F, rhs, x);
>> >
>> > // Now I would like to access the diagonal part of the matrix F, can I
>> > do something like:
>> >
>> > Mat_SeqAIJ *aa = (Mat_SeqAIJ*)F->data;
>> >
>> > Do you maybe have an example I could look at?
>> >
>> > Grateful in advance,
>> > Jelena
>>
>