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