1: !
  2: !
  3: !     Solves the problem A x - x^3 + 1 = 0 via Picard iteration
  4: !
  5: #include <petsc/finclude/petscsnes.h>
  6: module ex21fmodule
  7:   use petscsnes
  8:   implicit none
  9:   type userctx
 10:     Mat A
 11:   end type userctx
 12: contains
 13:   subroutine FormFunction(snes, x, f, user, ierr)
 14:     SNES snes
 15:     Vec x, f
 16:     type(userctx) user
 17:     PetscErrorCode, intent(out) :: ierr
 18:     PetscInt n
 19:     PetscScalar, pointer :: xx(:), ff(:)

 21:     PetscCallA(MatMult(user%A, x, f, ierr))
 22:     PetscCallA(VecGetArray(f, ff, ierr))
 23:     PetscCallA(VecGetArrayRead(x, xx, ierr))
 24:     PetscCallA(VecGetLocalSize(x, n, ierr))
 25:     ff(1:n) = ff(1:n) - xx(1:n)**4 + 1.0
 26:     PetscCallA(VecRestoreArray(f, ff, ierr))
 27:     PetscCallA(VecRestoreArrayRead(x, xx, ierr))
 28:   end subroutine

 30: ! The matrix is constant so no need to recompute it
 31:   subroutine FormJacobian(snes, x, jac, jacb, user, ierr)
 32:     SNES snes
 33:     Vec x
 34:     type(userctx) user
 35:     Mat jac, jacb
 36:     PetscErrorCode, intent(out) :: ierr

 38:     ierr = 0
 39:   end subroutine
 40: end module ex21fmodule

 42: program main
 43:   use petscsnes
 44:   use ex21fmodule
 45:   implicit none
 46:   SNES snes
 47:   PetscErrorCode ierr
 48:   Vec res, x
 49:   type(userctx) user
 50:   PetscScalar val

 52:   PetscCallA(PetscInitialize(ierr))

 54:   PetscCallA(MatCreateSeqAIJ(PETSC_COMM_SELF, 2_PETSC_INT_KIND, 2_PETSC_INT_KIND, 2_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, user%A, ierr))
 55:   val = 2.0
 56:   PetscCallA(MatSetValues(user%A, 1_PETSC_INT_KIND, [0_PETSC_INT_KIND], 1_PETSC_INT_KIND, [0_PETSC_INT_KIND], [val], INSERT_VALUES, ierr))
 57:   val = -1.0
 58:   PetscCallA(MatSetValues(user%A, 1_PETSC_INT_KIND, [0_PETSC_INT_KIND], 1_PETSC_INT_KIND, [1_PETSC_INT_KIND], [val], INSERT_VALUES, ierr))
 59:   val = -1.0
 60:   PetscCallA(MatSetValues(user%A, 1_PETSC_INT_KIND, [1_PETSC_INT_KIND], 1_PETSC_INT_KIND, [0_PETSC_INT_KIND], [val], INSERT_VALUES, ierr))
 61:   val = 1.0
 62:   PetscCallA(MatSetValues(user%A, 1_PETSC_INT_KIND, [1_PETSC_INT_KIND], 1_PETSC_INT_KIND, [1_PETSC_INT_KIND], [val], INSERT_VALUES, ierr))
 63:   PetscCallA(MatAssemblyBegin(user%A, MAT_FINAL_ASSEMBLY, ierr))
 64:   PetscCallA(MatAssemblyEnd(user%A, MAT_FINAL_ASSEMBLY, ierr))

 66:   PetscCallA(MatCreateVecs(user%A, x, res, ierr))

 68:   PetscCallA(SNESCreate(PETSC_COMM_SELF, snes, ierr))
 69:   PetscCallA(SNESSetPicard(snes, res, FormFunction, user%A, user%A, FormJacobian, user, ierr))
 70:   PetscCallA(SNESSetFromOptions(snes, ierr))
 71:   PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
 72:   PetscCallA(VecDestroy(x, ierr))
 73:   PetscCallA(VecDestroy(res, ierr))
 74:   PetscCallA(MatDestroy(user%A, ierr))
 75:   PetscCallA(SNESDestroy(snes, ierr))
 76:   PetscCallA(PetscFinalize(ierr))
 77: end

 79: !/*TEST
 80: !
 81: !   test:
 82: !     nsize: 1
 83: !     requires: !single
 84: !     args: -snes_monitor -snes_converged_reason -snes_view -pc_type lu
 85: !
 86: !TEST*/