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