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 PetscCallA(PetscInitialize(ierr)) 26 27 one = 1 28 zero = 0 29 two = 2 30 PetscCallA(MatCreateSeqAIJ(PETSC_COMM_SELF,two,two,two,PETSC_NULL_INTEGER,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)) 37 38 PetscCallA(MatCreateVecs(user%A,x,res,ierr)) 39 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 50 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(:) 59 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 70 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 80 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*/ 89