[petsc-users] KSP has an extra iteration when use shell matrix
Barry Smith
bsmith at petsc.dev
Mon Feb 5 18:39:41 CST 2024
You cannot do call MatShellSetContext(jac,X,ierr) in subroutine FormJacobianShell(snes,X,jac,B,dummy,ierr) because Fortran arguments are always pass by address and the address of the X being passed in may not be valid after the routine that called FormJacobianShell() has returned. This is why the GetContext works inside this function but not later in MyMult.
Instead add a Vec Xbase member to the MatCtx type and update that inside FormJacobShell(). Like
TYPE(MatCtx),POINTER :: ctxF_pt
MatShellGetContext(jac,ctxF_ptr,ierr)
ctxF_pt%Xbase = X
The reason this works is because ctxF_pt%Xbase = X copies the value of X (the PETSc vector) to Xbase, not the address of X.
Barry
> On Feb 5, 2024, at 4:18 PM, Yi Hu <y.hu at mpie.de> wrote:
>
> Dear Barry,
>
> the code is attached.
>
> Just to let you know. When I commented out MatShellSetContext() in FormJacobianShell(), then the code seems to work, meaning that the base vector is passed to shell matrix context behind the scene.
>
> Best regards,
>
> Yi
>
> On 2/5/24 19:09, Barry Smith wrote:
>>
>> Send the entire code.
>>
>>
>>> On Feb 4, 2024, at 4:43 PM, Yi Hu <y.hu at mpie.de> <mailto:y.hu at mpie.de> wrote:
>>>
>>> Thanks for your response. You are correct. I overlooked this step.
>>>
>>> Now I am trying to correct my "shell matrix approach" for ex1f.F90 of snes solver (https://github.com/hyharry/small_petsc_test/blob/master/test_shell_jac/ex1f.F90). I realized that I need to record the base vector X in the context of shell matrix and then use this info to carry MyMult. However, the context cannot be obtained through MatShellGetContext(). Here are the critical parts of my code.
>>>
>>> INTERFACE MatCreateShell
>>> SUBROUTINE MatCreateShell(comm,mloc,nloc,m,n,ctx,mat,ierr)
>>> USE solver_context
>>> MPI_Comm :: comm
>>> PetscInt :: mloc,nloc,m,n
>>> Vec :: ctx
>>> Mat :: mat
>>> PetscErrorCode :: ierr
>>> END SUBROUTINE MatCreateShell
>>> END INTERFACE MatCreateShell
>>>
>>> INTERFACE MatShellSetContext
>>> SUBROUTINE MatShellSetContext(mat,ctx,ierr)
>>> USE solver_context
>>> Mat :: mat
>>> !TYPE(MatCtx) :: ctx
>>> Vec :: ctx
>>> PetscErrorCode :: ierr
>>> END SUBROUTINE MatShellSetContext
>>> END INTERFACE MatShellSetContext
>>>
>>> INTERFACE MatShellGetContext
>>> SUBROUTINE MatShellGetContext(mat,ctx,ierr)
>>> USE solver_context
>>> Mat :: mat
>>> Vec, Pointer :: ctx
>>> PetscErrorCode :: ierr
>>> END SUBROUTINE MatShellGetContext
>>> END INTERFACE MatShellGetContext
>>>
>>> in my FormShellJacobian() I did
>>>
>>> subroutine FormJacobianShell(snes,X,jac,B,dummy,ierr)
>>>
>>> ......
>>>
>>> call MatShellSetContext(jac,X,ierr)
>>>
>>> ......
>>>
>>> Then in MyMult() I tried to recover this context by
>>>
>>> call MatShellGetContext(J,x,ierr)
>>>
>>> call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)
>>>
>>> Then the program failed with
>>>
>>> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
>>> [0]PETSC ERROR: Null argument, when expecting valid pointer
>>> [0]PETSC ERROR: Null Pointer: Parameter # 1
>>>
>>> In MyMult, I actually defined x to be a pointer. So I am confused here.
>>>
>>> Best regards,
>>>
>>> Yi
>>>
>>> On 1/31/24 03:18, Barry Smith wrote:
>>>>
>>>> It is not running an extra KSP iteration. This "extra" matmult is normal and occurs in many of the SNESLineSearchApply_* functions, for example, https://petsc.org/release/src/snes/linesearch/impls/bt/linesearchbt.c.html#SNESLineSearchApply_BT It is used to decide if the Newton step results in sufficient decrease of the function value.
>>>>
>>>> Barry
>>>>
>>>>
>>>>
>>>>> On Jan 30, 2024, at 3:19 PM, Yi Hu <y.hu at mpie.de> <mailto:y.hu at mpie.de> wrote:
>>>>>
>>>>> Hello Barry,
>>>>>
>>>>> Thanks for your reply. The monitor options are fine. I actually meant my modification of snes tutorial ex1f.F90 does not work and has some unexpected behavior. I basically wanted to test if I can use a shell matrix as my jacobian (code is here https://github.com/hyharry/small_petsc_test/blob/master/test_shell_jac/ex1f.F90). After compile my modified version and run with these monitor options, it gives me the following,
>>>>>
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 0 SNES Function norm 6.041522986797e+00
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 6.041522986797e+00
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 5.065392549852e-16
>>>>> Linear solve converged due to CONVERGED_RTOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 1 SNES Function norm 3.512662245652e+00
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 3.512662245652e+00
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 6.230314124713e-16
>>>>> Linear solve converged due to CONVERGED_RTOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 2 SNES Function norm 8.969285922373e-01
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 8.969285922373e-01
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 0.000000000000e+00
>>>>> Linear solve converged due to CONVERGED_ATOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 3 SNES Function norm 4.863816734540e-01
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 4.863816734540e-01
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 0.000000000000e+00
>>>>> Linear solve converged due to CONVERGED_ATOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 4 SNES Function norm 3.512070785520e-01
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 3.512070785520e-01
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 0.000000000000e+00
>>>>> Linear solve converged due to CONVERGED_ATOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 5 SNES Function norm 2.769700293115e-01
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 2.769700293115e-01
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 1.104778916974e-16
>>>>> Linear solve converged due to CONVERGED_RTOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 6 SNES Function norm 2.055345318150e-01
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 2.055345318150e-01
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 0.000000000000e+00
>>>>> Linear solve converged due to CONVERGED_ATOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 7 SNES Function norm 1.267482220786e-01
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 1.267482220786e-01
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 1.498679601680e-17
>>>>> Linear solve converged due to CONVERGED_RTOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> 8 SNES Function norm 3.468150619264e-02
>>>>> ++++++++++++ in jac shell +++++++++++
>>>>> 0 KSP Residual norm 3.468150619264e-02
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> 1 KSP Residual norm 5.944160522951e-18
>>>>> Linear solve converged due to CONVERGED_RTOL iterations 1
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> ( in rhs )
>>>>> ( leave rhs )
>>>>> === start mymult ===
>>>>> === done mymult ===
>>>>> Nonlinear solve did not converge due to DIVERGED_LINE_SEARCH iterations 8
>>>>> Number of SNES iterations = 8
>>>>>
>>>>> After each "Linear solve converged due to CONVERGED_ATOL iterations", the code starts to do mymult again. So I thought it did an extra (unwanted) KSP iteration. I would like to ask if this extra iteration could be disabled, or maybe I am wrong about it.
>>>>>
>>>>> Best regards,
>>>>>
>>>>> Yi
>>>>>
>>>>> On 1/30/24 18:35, Barry Smith wrote:
>>>>>>
>>>>>> How do I see a difference? What does "hence ruin my previous converged KSP result" mean? A different answer at the end of the KSP solve?
>>>>>>
>>>>>> $ ./joe > joe.basic
>>>>>> ~/Src/petsc/src/ksp/ksp/tutorials (barry/2023-09-15/fix-log-pcmpi=) arch-fix-log-pcmpi
>>>>>> $ ./joe -ksp_monitor -ksp_converged_reason -snes_monitor > joe.monitor
>>>>>> ~/Src/petsc/src/ksp/ksp/tutorials (barry/2023-09-15/fix-log-pcmpi=) arch-fix-log-pcmpi
>>>>>> $ diff joe.basic joe.monitor
>>>>>> 0a1,36
>>>>>> > 0 SNES Function norm 6.041522986797e+00
>>>>>> > 0 KSP Residual norm 6.041522986797e+00
>>>>>> > 1 KSP Residual norm 5.065392549852e-16
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>> > 1 SNES Function norm 3.512662245652e+00
>>>>>> > 0 KSP Residual norm 3.512662245652e+00
>>>>>> > 1 KSP Residual norm 6.230314124713e-16
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>> > 2 SNES Function norm 8.969285922373e-01
>>>>>> > 0 KSP Residual norm 8.969285922373e-01
>>>>>> > 1 KSP Residual norm 0.000000000000e+00
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>> > 3 SNES Function norm 4.863816734540e-01
>>>>>> > 0 KSP Residual norm 4.863816734540e-01
>>>>>> > 1 KSP Residual norm 0.000000000000e+00
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>> > 4 SNES Function norm 3.512070785520e-01
>>>>>> > 0 KSP Residual norm 3.512070785520e-01
>>>>>> > 1 KSP Residual norm 0.000000000000e+00
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>> > 5 SNES Function norm 2.769700293115e-01
>>>>>> > 0 KSP Residual norm 2.769700293115e-01
>>>>>> > 1 KSP Residual norm 1.104778916974e-16
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>> > 6 SNES Function norm 2.055345318150e-01
>>>>>> > 0 KSP Residual norm 2.055345318150e-01
>>>>>> > 1 KSP Residual norm 1.535110861002e-17
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>> > 7 SNES Function norm 1.267482220786e-01
>>>>>> > 0 KSP Residual norm 1.267482220786e-01
>>>>>> > 1 KSP Residual norm 1.498679601680e-17
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>> > 8 SNES Function norm 3.468150619264e-02
>>>>>> > 0 KSP Residual norm 3.468150619264e-02
>>>>>> > 1 KSP Residual norm 5.944160522951e-18
>>>>>> > Linear solve converged due to CONVERGED_RTOL_NORMAL iterations 1
>>>>>>
>>>>>>
>>>>>>
>>>>>>> On Jan 30, 2024, at 11:19 AM, Yi Hu <y.hu at mpie.de> <mailto:y.hu at mpie.de> wrote:
>>>>>>>
>>>>>>> Dear PETSc team,
>>>>>>>
>>>>>>> I am still trying to sort out my previous thread https://lists.mcs.anl.gov/pipermail/petsc-users/2024-January/050079.html using a minimal working example. However, I encountered another problem. Basically I combined the basic usage of SNES solver and shell matrix and tried to make it work. The jacobian of my snes is replaced by a customized MATOP_MULT. The minimal example code can be viewed here https://github.com/hyharry/small_petsc_test/blob/master/test_shell_jac/ex1f.F90
>>>>>>>
>>>>>>> When running with -ksp_monitor -ksp_converged_reason, it shows an extra mymult step, and hence ruin my previous converged KSP result. Implement a customized converged call-back also does not help. I am wondering how to skip this extra ksp iteration. Could anyone help me on this?
>>>>>>>
>>>>>>> Thanks for your help.
>>>>>>>
>>>>>>> Best wishes,
>>>>>>> Yi
>>>>>>>
>>>>>>>
>>>>>>> -------------------------------------------------
>>>>>>> Stay up to date and follow us on LinkedIn, Twitter and YouTube.
>>>>>>>
>>>>>>> Max-Planck-Institut für Eisenforschung GmbH
>>>>>>> Max-Planck-Straße 1
>>>>>>> D-40237 Düsseldorf
>>>>>>>
>>>>>>> Handelsregister B 2533
>>>>>>> Amtsgericht Düsseldorf
>>>>>>>
>>>>>>> Geschäftsführung
>>>>>>> Prof. Dr. Gerhard Dehm
>>>>>>> Prof. Dr. Jörg Neugebauer
>>>>>>> Prof. Dr. Dierk Raabe
>>>>>>> Dr. Kai de Weldige
>>>>>>>
>>>>>>> Ust.-Id.-Nr.: DE 11 93 58 514
>>>>>>> Steuernummer: 105 5891 1000
>>>>>>>
>>>>>>>
>>>>>>> Please consider that invitations and e-mails of our institute are
>>>>>>> only valid if they end with …@mpie.de.
>>>>>>> If you are not sure of the validity please contact rco at mpie.de <mailto:rco at mpie.de>
>>>>>>>
>>>>>>> Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails
>>>>>>> aus unserem Haus nur mit der Endung …@mpie.de gültig sind.
>>>>>>> In Zweifelsfällen wenden Sie sich bitte an rco at mpie.de <mailto:rco at mpie.de>
>>>>>>> -------------------------------------------------
>>>>>>
>>>>>
>>>>>
>>>>> -------------------------------------------------
>>>>> Stay up to date and follow us on LinkedIn, Twitter and YouTube.
>>>>>
>>>>> Max-Planck-Institut für Eisenforschung GmbH
>>>>> Max-Planck-Straße 1
>>>>> D-40237 Düsseldorf
>>>>>
>>>>> Handelsregister B 2533
>>>>> Amtsgericht Düsseldorf
>>>>>
>>>>> Geschäftsführung
>>>>> Prof. Dr. Gerhard Dehm
>>>>> Prof. Dr. Jörg Neugebauer
>>>>> Prof. Dr. Dierk Raabe
>>>>> Dr. Kai de Weldige
>>>>>
>>>>> Ust.-Id.-Nr.: DE 11 93 58 514
>>>>> Steuernummer: 105 5891 1000
>>>>>
>>>>>
>>>>> Please consider that invitations and e-mails of our institute are
>>>>> only valid if they end with …@mpie.de.
>>>>> If you are not sure of the validity please contact rco at mpie.de <mailto:rco at mpie.de>
>>>>>
>>>>> Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails
>>>>> aus unserem Haus nur mit der Endung …@mpie.de gültig sind.
>>>>> In Zweifelsfällen wenden Sie sich bitte an rco at mpie.de <mailto:rco at mpie.de>
>>>>> -------------------------------------------------
>>>>
>>>
>>>
>>> -------------------------------------------------
>>> Stay up to date and follow us on LinkedIn, Twitter and YouTube.
>>>
>>> Max-Planck-Institut für Eisenforschung GmbH
>>> Max-Planck-Straße 1
>>> D-40237 Düsseldorf
>>>
>>> Handelsregister B 2533
>>> Amtsgericht Düsseldorf
>>>
>>> Geschäftsführung
>>> Prof. Dr. Gerhard Dehm
>>> Prof. Dr. Jörg Neugebauer
>>> Prof. Dr. Dierk Raabe
>>> Dr. Kai de Weldige
>>>
>>> Ust.-Id.-Nr.: DE 11 93 58 514
>>> Steuernummer: 105 5891 1000
>>>
>>>
>>> Please consider that invitations and e-mails of our institute are
>>> only valid if they end with …@mpie.de.
>>> If you are not sure of the validity please contact rco at mpie.de <mailto:rco at mpie.de>
>>>
>>> Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails
>>> aus unserem Haus nur mit der Endung …@mpie.de gültig sind.
>>> In Zweifelsfällen wenden Sie sich bitte an rco at mpie.de <mailto:rco at mpie.de>
>>> -------------------------------------------------
>>
>
>
> -------------------------------------------------
> Stay up to date and follow us on LinkedIn, Twitter and YouTube.
>
> Max-Planck-Institut für Eisenforschung GmbH
> Max-Planck-Straße 1
> D-40237 Düsseldorf
>
> Handelsregister B 2533
> Amtsgericht Düsseldorf
>
> Geschäftsführung
> Prof. Dr. Gerhard Dehm
> Prof. Dr. Jörg Neugebauer
> Prof. Dr. Dierk Raabe
> Dr. Kai de Weldige
>
> Ust.-Id.-Nr.: DE 11 93 58 514
> Steuernummer: 105 5891 1000
>
>
> Please consider that invitations and e-mails of our institute are
> only valid if they end with …@mpie.de.
> If you are not sure of the validity please contact rco at mpie.de
>
> Bitte beachten Sie, dass Einladungen zu Veranstaltungen und E-Mails
> aus unserem Haus nur mit der Endung …@mpie.de gültig sind.
> In Zweifelsfällen wenden Sie sich bitte an rco at mpie.de
> -------------------------------------------------
> <ex1f.F90>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20240205/5738878d/attachment-0001.html>
More information about the petsc-users
mailing list