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 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(:) 63 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 74 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 84 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*/ 93