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_ARRAY,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*/