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