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: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
26: if (ierr .ne. 0) then
27: print*,'Unable to initialize PETSc'
28: stop
29: endif
31: one = 1
32: zero = 0
33: two = 2
34: call MatCreateSeqAIJ(PETSC_COMM_SELF,two,two,two,PETSC_NULL_INTEGER,user%A,ierr)
35: val = 2.0; call MatSetValues(user%A,one,zero,one,zero,val,INSERT_VALUES,ierr)
36: val = -1.0; call MatSetValues(user%A,one,zero,one,one,val,INSERT_VALUES,ierr)
37: val = -1.0; call MatSetValues(user%A,one,one,one,zero,val,INSERT_VALUES,ierr)
38: val = 1.0; call MatSetValues(user%A,one,one,one,one,val,INSERT_VALUES,ierr)
39: call MatAssemblyBegin(user%A,MAT_FINAL_ASSEMBLY,ierr)
40: call MatAssemblyEnd(user%A,MAT_FINAL_ASSEMBLY,ierr)
42: call MatCreateVecs(user%A,x,res,ierr)
44: call SNESCreate(PETSC_COMM_SELF,snes, ierr)
45: call SNESSetPicard(snes, res, FormFunction, user%A, user%A, FormJacobian, user, ierr)
46: call SNESSetFromOptions(snes,ierr)
47: call SNESSolve(snes, PETSC_NULL_VEC, x, ierr)
48: call VecDestroy(x,ierr)
49: call VecDestroy(res,ierr)
50: call MatDestroy(user%A,ierr)
51: call SNESDestroy(snes,ierr)
52: call PetscFinalize(ierr)
53: end
55: subroutine FormFunction(snes, x, f, user, ierr)
56: use ex21fmodule
57: SNES snes
58: Vec x, f
59: type(userctx) user
60: PetscErrorCode ierr
61: PetscInt i,n
62: PetscScalar, pointer :: xx(:),ff(:)
64: call MatMult(user%A, x, f, ierr)
65: call VecGetArrayF90(f,ff,ierr)
66: call VecGetArrayReadF90(x,xx,ierr)
67: call VecGetLocalSize(x,n,ierr)
68: do 10, i=1,n
69: ff(i) = ff(i) - xx(i)*xx(i)*xx(i)*xx(i) + 1.0
70: 10 continue
71: call VecRestoreArrayF90(f,ff,ierr)
72: call VecRestoreArrayReadF90(x,xx,ierr)
73: end subroutine
75: ! The matrix is constant so no need to recompute it
76: subroutine FormJacobian(snes, x, jac, jacb, user, ierr)
77: use ex21fmodule
78: SNES snes
79: Vec x
80: type(userctx) user
81: Mat jac, jacb
82: PetscErrorCode ierr
83: end subroutine
85: !/*TEST
86: !
87: ! test:
88: ! nsize: 1
89: ! requires: !single
90: ! args: -snes_monitor -snes_converged_reason -snes_view -pc_type lu
91: !
92: !TEST*/