1*df2570feSBarry Smith! 2*df2570feSBarry Smith! 3*df2570feSBarry Smith! Solves the problem A x - x^3 + 1 = 0 via Picard iteration 4*df2570feSBarry Smith! 5*df2570feSBarry Smith module ex21fmodule 6*df2570feSBarry Smith use petscsnes 7*df2570feSBarry Smith#include <petsc/finclude/petscsnes.h> 8*df2570feSBarry Smith type userctx 9*df2570feSBarry Smith Mat A 10*df2570feSBarry Smith end type userctx 11*df2570feSBarry Smith end module ex21fmodule 12*df2570feSBarry Smith 13*df2570feSBarry Smith program main 14*df2570feSBarry Smith#include <petsc/finclude/petscsnes.h> 15*df2570feSBarry Smith use ex21fmodule 16*df2570feSBarry Smith implicit none 17*df2570feSBarry Smith SNES snes 18*df2570feSBarry Smith PetscErrorCode ierr 19*df2570feSBarry Smith Vec res,x 20*df2570feSBarry Smith type(userctx) user 21*df2570feSBarry Smith PetscScalar val 22*df2570feSBarry Smith PetscInt one,zero,two 23*df2570feSBarry Smith external FormFunction,FormJacobian 24*df2570feSBarry Smith 25*df2570feSBarry Smith call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 26*df2570feSBarry Smith if (ierr .ne. 0) then 27*df2570feSBarry Smith print*,'Unable to initialize PETSc' 28*df2570feSBarry Smith stop 29*df2570feSBarry Smith endif 30*df2570feSBarry Smith 31*df2570feSBarry Smith one = 1 32*df2570feSBarry Smith zero = 0 33*df2570feSBarry Smith two = 2 34*df2570feSBarry Smith call MatCreateSeqAIJ(PETSC_COMM_SELF,two,two,two,PETSC_NULL_INTEGER,user%A,ierr) 35*df2570feSBarry Smith val = 2.0; call MatSetValues(user%A,one,zero,one,zero,val,INSERT_VALUES,ierr) 36*df2570feSBarry Smith val = -1.0; call MatSetValues(user%A,one,zero,one,one,val,INSERT_VALUES,ierr) 37*df2570feSBarry Smith val = -1.0; call MatSetValues(user%A,one,one,one,zero,val,INSERT_VALUES,ierr) 38*df2570feSBarry Smith val = 1.0; call MatSetValues(user%A,one,one,one,one,val,INSERT_VALUES,ierr) 39*df2570feSBarry Smith call MatAssemblyBegin(user%A,MAT_FINAL_ASSEMBLY,ierr) 40*df2570feSBarry Smith call MatAssemblyEnd(user%A,MAT_FINAL_ASSEMBLY,ierr) 41*df2570feSBarry Smith 42*df2570feSBarry Smith call MatCreateVecs(user%A,x,res,ierr) 43*df2570feSBarry Smith 44*df2570feSBarry Smith call SNESCreate(PETSC_COMM_SELF,snes, ierr) 45*df2570feSBarry Smith call SNESSetPicard(snes, res, FormFunction, user%A, user%A, FormJacobian, user, ierr) 46*df2570feSBarry Smith call SNESSetFromOptions(snes,ierr) 47*df2570feSBarry Smith call SNESSolve(snes, PETSC_NULL_VEC, x, ierr) 48*df2570feSBarry Smith call VecDestroy(x,ierr) 49*df2570feSBarry Smith call VecDestroy(res,ierr) 50*df2570feSBarry Smith call MatDestroy(user%A,ierr) 51*df2570feSBarry Smith call SNESDestroy(snes,ierr) 52*df2570feSBarry Smith call PetscFinalize(ierr) 53*df2570feSBarry Smith end 54*df2570feSBarry Smith 55*df2570feSBarry Smith 56*df2570feSBarry Smith subroutine FormFunction(snes, x, f, user, ierr) 57*df2570feSBarry Smith use ex21fmodule 58*df2570feSBarry Smith SNES snes 59*df2570feSBarry Smith Vec x, f 60*df2570feSBarry Smith type(userctx) user 61*df2570feSBarry Smith PetscErrorCode ierr 62*df2570feSBarry Smith PetscInt i,n 63*df2570feSBarry Smith PetscScalar, pointer :: xx(:),ff(:) 64*df2570feSBarry Smith 65*df2570feSBarry Smith call MatMult(user%A, x, f, ierr) 66*df2570feSBarry Smith call VecGetArrayF90(f,ff,ierr) 67*df2570feSBarry Smith call VecGetArrayReadF90(x,xx,ierr) 68*df2570feSBarry Smith call VecGetLocalSize(x,n,ierr) 69*df2570feSBarry Smith do 10, i=1,n 70*df2570feSBarry Smith ff(i) = ff(i) - xx(i)*xx(i)*xx(i)*xx(i) + 1.0 71*df2570feSBarry Smith 10 continue 72*df2570feSBarry Smith call VecRestoreArrayF90(f,ff,ierr) 73*df2570feSBarry Smith call VecRestoreArrayReadF90(x,xx,ierr) 74*df2570feSBarry Smith end subroutine 75*df2570feSBarry Smith 76*df2570feSBarry Smith! The matrix is constant so no need to recompute it 77*df2570feSBarry Smith subroutine FormJacobian(snes, x, jac, jacb, user, ierr) 78*df2570feSBarry Smith use ex21fmodule 79*df2570feSBarry Smith SNES snes 80*df2570feSBarry Smith Vec x 81*df2570feSBarry Smith type(userctx) user 82*df2570feSBarry Smith Mat jac, jacb 83*df2570feSBarry Smith PetscErrorCode ierr 84*df2570feSBarry Smith end subroutine 85*df2570feSBarry Smith 86*df2570feSBarry Smith!/*TEST 87*df2570feSBarry Smith! 88*df2570feSBarry Smith! test: 89*df2570feSBarry Smith! nsize: 1 90*df2570feSBarry Smith! requires: !single 91*df2570feSBarry Smith! args: -snes_monitor -snes_converged_reason -snes_view -pc_type lu 92*df2570feSBarry Smith! 93*df2570feSBarry Smith!TEST*/ 94