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