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 <[email protected]> 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 <[email protected]> <mailto:[email protected]> >>> 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 <[email protected]> <mailto:[email protected]> >>>>> 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 <[email protected]> >>>>>>> <mailto:[email protected]> 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 [email protected] >>>>>>> <mailto:[email protected]> >>>>>>> >>>>>>> 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 [email protected] >>>>>>> <mailto:[email protected]> >>>>>>> ------------------------------------------------- >>>>>> >>>>> >>>>> >>>>> ------------------------------------------------- >>>>> 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 [email protected] >>>>> <mailto:[email protected]> >>>>> >>>>> 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 [email protected] >>>>> <mailto:[email protected]> >>>>> ------------------------------------------------- >>>> >>> >>> >>> ------------------------------------------------- >>> 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 [email protected] >>> <mailto:[email protected]> >>> >>> 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 [email protected] >>> <mailto:[email protected]> >>> ------------------------------------------------- >> > > > ------------------------------------------------- > 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 [email protected] > > 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 [email protected] > ------------------------------------------------- > <ex1f.F90>
