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 12 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 24 25 call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 26 if (ierr .ne. 0) then 27 print*,'Unable to initialize PETSc' 28 stop 29 endif 30 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) 41 42 call MatCreateVecs(user%A,x,res,ierr) 43 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 54 55 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(:) 64 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 75 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 85 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*/ 94