Actual source code: ex21f.F90

petsc-3.13.4 2020-08-01
Report Typos and Errors
  1: !
  2: !
  3: !     Solves the problem A x - x^3 + 1 = 0 via Picard iteration
  4: !
  5:       module ex21fmodule
  6:         use petscsnes
  7:  #include <petsc/finclude/petscsnes.h>
  8:         type userctx
  9:           Mat A
 10:         end type userctx
 11:       end module ex21fmodule

 13:       program main
 14:  #include <petsc/finclude/petscsnes.h>
 15:       use ex21fmodule
 16:       implicit none
 17:       SNES snes
 18:       PetscErrorCode ierr
 19:       Vec res,x
 20:       type(userctx) user
 21:       PetscScalar val
 22:       PetscInt one,zero,two
 23:       external FormFunction,FormJacobian

 25:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 26:       if (ierr .ne. 0) then
 27:         print*,'Unable to initialize PETSc'
 28:         stop
 29:       endif

 31:       one = 1
 32:       zero = 0
 33:       two = 2
 34:       call MatCreateSeqAIJ(PETSC_COMM_SELF,two,two,two,PETSC_NULL_INTEGER,user%A,ierr)
 35:       val = 2.0; call MatSetValues(user%A,one,zero,one,zero,val,INSERT_VALUES,ierr)
 36:       val = -1.0; call MatSetValues(user%A,one,zero,one,one,val,INSERT_VALUES,ierr)
 37:       val = -1.0; call MatSetValues(user%A,one,one,one,zero,val,INSERT_VALUES,ierr)
 38:       val = 1.0; call MatSetValues(user%A,one,one,one,one,val,INSERT_VALUES,ierr)
 39:       call MatAssemblyBegin(user%A,MAT_FINAL_ASSEMBLY,ierr)
 40:       call MatAssemblyEnd(user%A,MAT_FINAL_ASSEMBLY,ierr)

 42:       call MatCreateVecs(user%A,x,res,ierr)

 44:       call SNESCreate(PETSC_COMM_SELF,snes, ierr)
 45:       call SNESSetPicard(snes, res, FormFunction, user%A, user%A, FormJacobian, user, ierr)
 46:       call SNESSetFromOptions(snes,ierr)
 47:       call SNESSolve(snes, PETSC_NULL_VEC, x, ierr)
 48:       call VecDestroy(x,ierr)
 49:       call VecDestroy(res,ierr)
 50:       call MatDestroy(user%A,ierr)
 51:       call SNESDestroy(snes,ierr)
 52:       call PetscFinalize(ierr)
 53:       end


 56:       subroutine FormFunction(snes, x, f, user, ierr)
 57:       use ex21fmodule
 58:       SNES snes
 59:       Vec  x, f
 60:       type(userctx) user
 61:       PetscErrorCode ierr
 62:       PetscInt i,n
 63:       PetscScalar, pointer :: xx(:),ff(:)

 65:       call MatMult(user%A, x, f, ierr)
 66:       call VecGetArrayF90(f,ff,ierr)
 67:       call VecGetArrayReadF90(x,xx,ierr)
 68:       call VecGetLocalSize(x,n,ierr)
 69:       do 10, i=1,n
 70:          ff(i) = ff(i) - xx(i)*xx(i)*xx(i)*xx(i) + 1.0
 71:  10   continue
 72:       call VecRestoreArrayF90(f,ff,ierr)
 73:       call VecRestoreArrayReadF90(x,xx,ierr)
 74:       end subroutine

 76: !      The matrix is constant so no need to recompute it
 77:       subroutine FormJacobian(snes, x, jac, jacb, user, ierr)
 78:       use ex21fmodule
 79:       SNES snes
 80:       Vec  x
 81:       type(userctx) user
 82:       Mat  jac, jacb
 83:       PetscErrorCode ierr
 84:       end subroutine

 86: !/*TEST
 87: !
 88: !   test:
 89: !     nsize: 1
 90: !     requires: !single
 91: !     args: -snes_monitor -snes_converged_reason -snes_view -pc_type lu
 92: !
 93: !TEST*/