[petsc-users] PETSc-MUMPS interface, numeric and symbolic factorisation

Matthew Knepley knepley at gmail.com
Fri May 11 12:10:38 CDT 2018 

On Fri, May 11, 2018 at 1:02 PM, Y. Shidi <ys453 at cam.ac.uk> wrote:

> Thank you very much for your reply, Barry.
>
>    This is a bug in PETSc. Since you are providing a new matrix with
>> the same "state" value as the previous matrix the PC code the
>> following code
>>
> So what you mean is that every time I change the value in the matrix,
> the PETSc only determines if the nonzero pattern change but not the
> values, and if it is unchanged neither of symbolic and numeric
> happens.
>

No, that is not what Barry is saying.

PETSc looks at the matrix.
If the structure has changed, it does symbolic and numeric factorization.
If only values have changes, it does numeric factorization.

HOWEVER, you gave it a new matrix with accidentally the same state marker,
so it thought nothing had changed. We will fix this by also checking the
pointer.
For now, if you give the same Mat object back, it will do what you expect.

  Matt


> I found the following code:
>
>   if (!pc->setupcalled) {
>     ierr            = PetscInfo(pc,"Setting up PC for first
> time\n");CHKERRQ(ierr);
>     pc->flag        = DIFFERENT_NONZERO_PATTERN;
>   } else if (matstate == pc->matstate) {
>     ierr = PetscInfo(pc,"Leaving PC with identical preconditioner since
> operator is unchanged\n");CHKERRQ(ierr);
>     PetscFunctionReturn(0);
>   } else {
>     if (matnonzerostate > pc->matnonzerostate) {
>        ierr = PetscInfo(pc,"Setting up PC with different nonzero
> pattern\n");CHKERRQ(ierr);
>        pc->flag            = DIFFERENT_NONZERO_PATTERN;
>     } else {
>       ierr = PetscInfo(pc,"Setting up PC with same nonzero
> pattern\n");CHKERRQ(ierr);
>       pc->flag            = SAME_NONZERO_PATTERN;
>     }
>   }
>
> and I commend out "else if (matstate == pc->matstate){}", so it
> will do "Setting up PC with same nonzero pattern\n"; and it seems
> work in my case, only "MatFactorNumeric_MUMPS()" is calling in the
> subsequent iterations. But I am not quite sure, need some more tests.
>
> Thank you very much for your help indeed.
>
> Kind Regards,
> Shidi
>
>
>
> On 2018-05-11 16:13, Smith, Barry F. wrote:
>
>> On May 11, 2018, at 8:14 AM, Y. Shidi <ys453 at cam.ac.uk> wrote:
>>>
>>> Thank you for your reply.
>>>
>>> How are you changing the matrix? Do you remember to assemble?
>>>>
>>> I use MatCreateMPIAIJWithArrays() to create the matrix,
>>> and after that I call MatAssemblyBegin() and MatAssemblyEnd().
>>>
>>
>>   If you use MatCreateMPIAIJWithArrays() you don't need to call
>> MatAssemblyBegin() and MatAssemblyEnd().
>>
>> But I actually destroy the matrix at the end of each iteration
>>> and create the matrix at the beginning of each iteration.
>>>
>>
>>    This is a bug in PETSc. Since you are providing a new matrix with
>> the same "state" value as the previous matrix the PC code the
>> following code
>> kicks in:
>>
>>   ierr = PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
>> CHKERRQ(ierr);
>>   ierr = MatGetNonzeroState(pc->pmat,&matnonzerostate);CHKERRQ(ierr);
>>   if (!pc->setupcalled) {
>>     ierr            = PetscInfo(pc,"Setting up PC for first
>> time\n");CHKERRQ(ierr);
>>     pc->flag        = DIFFERENT_NONZERO_PATTERN;
>>   } else if (matstate == pc->matstate) {
>>     ierr = PetscInfo(pc,"Leaving PC with identical preconditioner
>> since operator is unchanged\n");CHKERRQ(ierr);
>>     PetscFunctionReturn(0);
>>
>> and it returns without refactoring.
>>
>>    We need an additional check that the matrix also remains the same.
>>
>>     We will also need a test example that reproduces the problem to
>> confirm that we have fixed it.
>>
>>   Barry
>>
>>
>>> Cheers,
>>> Shidi
>>>
>>> On 2018-05-11 12:59, Matthew Knepley wrote:
>>>
>>>> On Fri, May 11, 2018 at 7:14 AM, Y. Shidi <ys453 at cam.ac.uk> wrote:
>>>>
>>>>> Dear Matt,
>>>>> Thank you for your help last time.
>>>>> I want to get more detail about the Petsc-MUMPS factorisation;
>>>>> so I go to look the code "/src/mat/impls/aij/mpi/mumps/mumps.c".
>>>>> And I found the following functions are quite important to
>>>>> the question:
>>>>> PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS
>>>>> r,const MatFactorInfo *info);
>>>>> PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const
>>>>> MatFactorInfo *info);
>>>>> PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x);
>>>>> I print some sentence to trace when these functions are called.
>>>>> Then I test my code; the values in the matrix is changing but the
>>>>> structure stays the same. Below is the output.
>>>>> We can see that at 0th step, all the symbolic, numeric and solve
>>>>> are called; in the subsequent steps only the solve stage is called,
>>>>> the numeric step is not called.
>>>>>
>>>> How are you changing the matrix? Do you remember to assemble?
>>>>  Matt
>>>>
>>>>> Iteration 0 Step 0.0005 Time 0.0005
>>>>> [INFO]: Direct Solver setup
>>>>> MatCholeskyFactorSymbolic_MUMPS
>>>>> finish MatCholeskyFactorSymbolic_MUMPS
>>>>> MatFactorNumeric_MUMPS
>>>>> finish MatFactorNumeric_MUMPS
>>>>> MatSolve_MUMPS
>>>>> Iteration 1 Step 0.0005 Time 0.0005
>>>>> MatSolve_MUMPS
>>>>> Iteration 2 Step 0.0005 Time 0.001
>>>>> MatSolve_MUMPS
>>>>> [INFO]: End of program!!!
>>>>> I am wondering if there is any possibility to split the numeric
>>>>> and solve stage (as you mentioned using KSPSolve).
>>>>> Thank you very much indeed.
>>>>> Kind Regards,
>>>>> Shidi
>>>>> On 2018-05-04 21:10, Y. Shidi wrote:
>>>>> Thank you very much for your reply.
>>>>> That is really clear.
>>>>> Kind Regards,
>>>>> Shidi
>>>>> On 2018-05-04 21:05, Matthew Knepley wrote:
>>>>> On Fri, May 4, 2018 at 3:54 PM, Y. Shidi <ys453 at cam.ac.uk> wrote:
>>>>> Dear Matt,
>>>>> Thank you very much for your reply!
>>>>> So what you mean is that I can just do the KSPSolve() every
>>>>> iteration
>>>>> once the MUMPS is set?
>>>>> Yes.
>>>>> That means inside the KSPSolve() the numerical factorization is
>>>>> performed. If that is the case, it seems that the ksp object is
>>>>> not changed when the values in the matrix are changed.
>>>>> Yes.
>>>>> Or do I need to call both KSPSetOperators() and KSPSolve()?
>>>>> If you do SetOperators, it will redo the factorization. If you do
>>>>> not,
>>>>> it will look
>>>>> at the Mat object, determine that the structure has not changed,
>>>>> and
>>>>> just redo
>>>>> the numerical factorization.
>>>>> Thanks,
>>>>> Matt
>>>>> On 2018-05-04 14:44, Matthew Knepley wrote:
>>>>> On Fri, May 4, 2018 at 9:40 AM, Y. Shidi <ys453 at cam.ac.uk> wrote:
>>>>> Dear PETSc users,
>>>>> I am currently using MUMPS to solve linear systems directly.
>>>>> Generally, we use ICNTL(7) or ICNTL(29) to do the preprocessing
>>>>> step and then solve the system.
>>>>> In my code, the values in the matrix is changed in each iteration,
>>>>> but the structure of the matrix stays the same, which means the
>>>>> performance can be improved if symbolic factorisation is only
>>>>> performed once. Hence, it is necessary to split the symbolic
>>>>> and numeric factorisation. However, I cannot find a specific step
>>>>> (control parameter) to perform the numeric factorisation.
>>>>> I have used ICNTL(3) and ICNTL(4) to print the MUMPS information,
>>>>> it seems that the symbolic and numeric factorisation always perform
>>>>> together.
>>>>> If you use KSPSolve instead, it will automatically preserve the
>>>>> symbolic
>>>>> factorization.
>>>>> Thanks,
>>>>> Matt
>>>>> So I am wondering if anyone has an idea about it.
>>>>> Below is how I set up MUMPS solver:
>>>>> PC pc;
>>>>> PetscBool flg_mumps, flg_mumps_ch;
>>>>> flg_mumps = PETSC_FALSE;
>>>>> flg_mumps_ch = PETSC_FALSE;
>>>>> PetscOptionsGetBool(NULL, NULL, "-use_mumps_lu", &flg_mumps,
>>>>> NULL);
>>>>> PetscOptionsGetBool(NULL, NULL, "-use_mumps_ch", &flg_mumps_ch,
>>>>> NULL);
>>>>> if(flg_mumps ||flg_mumps_ch)
>>>>> {
>>>>> KSPSetType(_ksp, KSPPREONLY);
>>>>> PetscInt ival,icntl;
>>>>> PetscReal val;
>>>>> KSPGetPC(_ksp, &pc);
>>>>> /// Set preconditioner type
>>>>> if(flg_mumps)
>>>>> {
>>>>> PCSetType(pc, PCLU);
>>>>> }
>>>>> else if(flg_mumps_ch)
>>>>> {
>>>>> MatSetOption(A, MAT_SPD, PETSC_TRUE);
>>>>> PCSetType(pc, PCCHOLESKY);
>>>>> }
>>>>> PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
>>>>> PCFactorSetUpMatSolverPackage(pc);
>>>>> PCFactorGetMatrix(pc, &_F);
>>>>> icntl = 7; ival = 0;
>>>>> MatMumpsSetIcntl( _F, icntl, ival );
>>>>> MatMumpsSetIcntl(_F, 3, 6);
>>>>> MatMumpsSetIcntl(_F, 4, 2);
>>>>> }
>>>>> KSPSetUp(_ksp);
>>>>> Kind Regards,
>>>>> Shidi
>>>>> --
>>>>> What most experimenters take for granted before they begin their
>>>>> experiments is infinitely more interesting than any results to
>>>>> which
>>>>> their experiments lead.
>>>>> -- Norbert Wiener
>>>>> https://www.cse.buffalo.edu/~knepley/ [1] [1] [1]
>>>>> Links:
>>>>> ------
>>>>> [1] http://www.caam.rice.edu/~mk51/ [2] [2]
>>>>> --
>>>>> What most experimenters take for granted before they begin their
>>>>> experiments is infinitely more interesting than any results to
>>>>> which
>>>>> their experiments lead.
>>>>> -- Norbert Wiener
>>>>> https://www.cse.buffalo.edu/~knepley/ [1] [2]
>>>>> Links:
>>>>> ------
>>>>> [1] https://www.cse.buffalo.edu/~knepley/ [1]
>>>>> [2] http://www.caam.rice.edu/~mk51/ [2]
>>>>>
>>>> --
>>>> What most experimenters take for granted before they begin their
>>>> experiments is infinitely more interesting than any results to which
>>>> their experiments lead.
>>>> -- Norbert Wiener
>>>> https://www.cse.buffalo.edu/~knepley/ [2]
>>>> Links:
>>>> ------
>>>> [1] https://www.cse.buffalo.edu/~knepley/
>>>> [2] http://www.caam.rice.edu/~mk51/
>>>>
>>>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180511/39ddc4c5/attachment-0001.html>

More information about the petsc-users mailing list