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 ierr
 18:     PetscInt i, 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:     do i = 1, n
 26:       ff(i) = ff(i) - xx(i)*xx(i)*xx(i)*xx(i) + 1.0
 27:     end do
 28:     PetscCallA(VecRestoreArray(f, ff, ierr))
 29:     PetscCallA(VecRestoreArrayRead(x, xx, ierr))
 30:   end subroutine

 32: !      The matrix is constant so no need to recompute it
 33:   subroutine FormJacobian(snes, x, jac, jacb, user, ierr)
 34:     SNES snes
 35:     Vec x
 36:     type(userctx) user
 37:     Mat jac, jacb
 38:     PetscErrorCode ierr
 39:   end subroutine
 40: end module ex21fmodule

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

 52:   PetscCallA(PetscInitialize(ierr))

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

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

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

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