1! 2! 3! Solves the problem A x - x^3 + 1 = 0 via Picard iteration 4! 5#include <petsc/finclude/petscsnes.h> 6#include <petsc/finclude/petscsnes.h> 7module ex21fmodule 8 use petscsnes 9 type userctx 10 Mat A 11 end type userctx 12end module ex21fmodule 13 14program main 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 i = 1, n 65 ff(i) = ff(i) - xx(i)*xx(i)*xx(i)*xx(i) + 1.0 66 end do 67 PetscCallA(VecRestoreArray(f, ff, ierr)) 68 PetscCallA(VecRestoreArrayRead(x, xx, ierr)) 69end subroutine 70 71! The matrix is constant so no need to recompute it 72subroutine 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 79end 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