Actual source code: ex21f.F90
petsc-3.13.4 2020-08-01
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
56: subroutine FormFunction(snes, x, f, user, ierr)
57: use ex21fmodule
58: SNES snes
59: Vec x, f
60: type(userctx) user
61: PetscErrorCode ierr
62: PetscInt i,n
63: PetscScalar, pointer :: xx(:),ff(:)
65: call MatMult(user%A, x, f, ierr)
66: call VecGetArrayF90(f,ff,ierr)
67: call VecGetArrayReadF90(x,xx,ierr)
68: call VecGetLocalSize(x,n,ierr)
69: do 10, i=1,n
70: ff(i) = ff(i) - xx(i)*xx(i)*xx(i)*xx(i) + 1.0
71: 10 continue
72: call VecRestoreArrayF90(f,ff,ierr)
73: call VecRestoreArrayReadF90(x,xx,ierr)
74: end subroutine
76: ! The matrix is constant so no need to recompute it
77: subroutine FormJacobian(snes, x, jac, jacb, user, ierr)
78: use ex21fmodule
79: SNES snes
80: Vec x
81: type(userctx) user
82: Mat jac, jacb
83: PetscErrorCode ierr
84: end subroutine
86: !/*TEST
87: !
88: ! test:
89: ! nsize: 1
90: ! requires: !single
91: ! args: -snes_monitor -snes_converged_reason -snes_view -pc_type lu
92: !
93: !TEST*/