Actual source code: ex21f.F90

  1: !
  2: !
  3: !     Solves the problem A x - x^3 + 1 = 0 via Picard iteration
  4: !
  5:       module ex21fmodule
  6:         use petscsnes
  7: #include <petsc/finclude/petscsnes.h>
  8:         type userctx
  9:           Mat A
 10:         end type userctx
 11:       end module ex21fmodule

 13:       program main
 14: #include <petsc/finclude/petscsnes.h>
 15:       use ex21fmodule
 16:       implicit none
 17:       SNES snes
 18:       PetscErrorCode ierr
 19:       Vec res,x
 20:       type(userctx) user
 21:       PetscScalar val
 22:       PetscInt one,zero,two
 23:       external FormFunction,FormJacobian

 25:       PetscCallA(PetscInitialize(ierr))

 27:       one = 1
 28:       zero = 0
 29:       two = 2
 30:       PetscCallA(MatCreateSeqAIJ(PETSC_COMM_SELF,two,two,two,PETSC_NULL_INTEGER,user%A,ierr))
 31:       val = 2.0; PetscCallA(MatSetValues(user%A,one,zero,one,zero,val,INSERT_VALUES,ierr))
 32:       val = -1.0; PetscCallA(MatSetValues(user%A,one,zero,one,one,val,INSERT_VALUES,ierr))
 33:       val = -1.0; PetscCallA(MatSetValues(user%A,one,one,one,zero,val,INSERT_VALUES,ierr))
 34:       val = 1.0; PetscCallA(MatSetValues(user%A,one,one,one,one,val,INSERT_VALUES,ierr))
 35:       PetscCallA(MatAssemblyBegin(user%A,MAT_FINAL_ASSEMBLY,ierr))
 36:       PetscCallA(MatAssemblyEnd(user%A,MAT_FINAL_ASSEMBLY,ierr))

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

 40:       PetscCallA(SNESCreate(PETSC_COMM_SELF,snes, ierr))
 41:       PetscCallA(SNESSetPicard(snes, res, FormFunction, user%A, user%A, FormJacobian, user, ierr))
 42:       PetscCallA(SNESSetFromOptions(snes,ierr))
 43:       PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
 44:       PetscCallA(VecDestroy(x,ierr))
 45:       PetscCallA(VecDestroy(res,ierr))
 46:       PetscCallA(MatDestroy(user%A,ierr))
 47:       PetscCallA(SNESDestroy(snes,ierr))
 48:       PetscCallA(PetscFinalize(ierr))
 49:       end

 51:       subroutine FormFunction(snes, x, f, user, ierr)
 52:       use ex21fmodule
 53:       SNES snes
 54:       Vec  x, f
 55:       type(userctx) user
 56:       PetscErrorCode ierr
 57:       PetscInt i,n
 58:       PetscScalar, pointer :: xx(:),ff(:)

 60:       PetscCallA(MatMult(user%A, x, f, ierr))
 61:       PetscCallA(VecGetArrayF90(f,ff,ierr))
 62:       PetscCallA(VecGetArrayReadF90(x,xx,ierr))
 63:       PetscCallA(VecGetLocalSize(x,n,ierr))
 64:       do 10, i=1,n
 65:          ff(i) = ff(i) - xx(i)*xx(i)*xx(i)*xx(i) + 1.0
 66:  10   continue
 67:       PetscCallA(VecRestoreArrayF90(f,ff,ierr))
 68:       PetscCallA(VecRestoreArrayReadF90(x,xx,ierr))
 69:       end subroutine

 71: !      The matrix is constant so no need to recompute it
 72:       subroutine FormJacobian(snes, x, jac, jacb, user, ierr)
 73:       use ex21fmodule
 74:       SNES snes
 75:       Vec  x
 76:       type(userctx) user
 77:       Mat  jac, jacb
 78:       PetscErrorCode ierr
 79:       end subroutine

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