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