1! 2! 3! Solves the problem A x - x^3 + 1 = 0 via Picard iteration 4! 5module ex21fmodule 6 use petscsnes 7#include <petsc/finclude/petscsnes.h> 8 type userctx 9 Mat A 10 end type userctx 11end module ex21fmodule 12 13program 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_ARRAY, 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)) 49end 50 51subroutine 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(VecGetArray(f, ff, ierr)) 62 PetscCallA(VecGetArrayRead(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 6610 continue 67 PetscCallA(VecRestoreArray(f, ff, ierr)) 68 PetscCallA(VecRestoreArrayRead(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