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

Smith, Barry F. bsmith at mcs.anl.gov
Fri May 11 14:08:54 CDT 2018 


> On May 11, 2018, at 12:10 PM, Matthew Knepley <knepley at gmail.com> wrote:
> 
> 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.

  Matt, what you say below is exactly what I was 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/

More information about the petsc-users mailing list