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