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