Actual source code: ex100f.F90

  1: program main
  2: #include "petsc/finclude/petscksp.h"
  3:   use petscksp

  5:   PetscInt N
  6:   PetscBool draw, flg
  7:   PetscReal rnorm, rtwo
  8:   PetscScalar one, mone
  9:   Mat A
 10:   Vec b, x, r
 11:   KSP ksp
 12:   PC pc
 13:   PetscErrorCode ierr
 14:   character*80 mattype

 16:   N = 100
 17:   draw = .FALSE.
 18:   one = 1.0
 19:   mone = -1.0
 20:   rtwo = 2.0

 22:   PetscCallA(PetscInitialize(ierr))
 23:   PetscCallA(PetscPythonInitialize(PETSC_NULL_CHARACTER, PETSC_NULL_CHARACTER, ierr))

 25:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-N', N, flg, ierr))
 26:   PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-draw', draw, flg, ierr))

 28:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 29:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr))
 30:   PetscCallA(MatSetType(A, 'python', ierr))
 31:   PetscCallA(MatPythonSetType(A, 'example100.py:Laplace1D', ierr))
 32:   PetscCallA(MatPythonGetType(A, mattype, ierr))
 33:   PetscCheckA(mattype == 'example100.py:Laplace1D', PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Error')
 34:   PetscCallA(MatSetUp(A, ierr))

 36:   PetscCallA(MatCreateVecs(A, x, b, ierr))
 37:   PetscCallA(VecSet(b, one, ierr))

 39:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
 40:   PetscCallA(KSPSetType(ksp, 'python', ierr))
 41:   PetscCallA(KSPPythonSetType(ksp, 'example100.py:ConjGrad', ierr))

 43:   PetscCallA(KSPGetPC(ksp, pc, ierr))
 44:   PetscCallA(PCSetType(pc, 'python', ierr))
 45:   PetscCallA(PCPythonSetType(pc, 'example100.py:Jacobi', ierr))

 47:   PetscCallA(KSPSetOperators(ksp, A, A, ierr))
 48:   PetscCallA(KSPSetFromOptions(ksp, ierr))
 49:   PetscCallA(KSPSolve(ksp, b, x, ierr))

 51:   PetscCallA(VecDuplicate(b, r, ierr))
 52:   PetscCallA(MatMult(A, x, r, ierr))
 53:   PetscCallA(VecAYPX(r, mone, b, ierr))
 54:   PetscCallA(VecNorm(r, NORM_2, rnorm, ierr))
 55:   print *, 'error norm = ', rnorm

 57:   if (draw) then
 58:     PetscCallA(VecView(x, PETSC_VIEWER_DRAW_WORLD, ierr))
 59:     PetscCallA(PetscSleep(rtwo, ierr))
 60:   end if

 62:   PetscCallA(VecDestroy(x, ierr))
 63:   PetscCallA(VecDestroy(b, ierr))
 64:   PetscCallA(VecDestroy(r, ierr))
 65:   PetscCallA(MatDestroy(A, ierr))
 66:   PetscCallA(KSPDestroy(ksp, ierr))

 68:   PetscCallA(PetscFinalize(ierr))
 69: end

 71: !/*TEST
 72: !
 73: !    test:
 74: !      requires: petsc4py
 75: !      localrunfiles: example100.py
 76: !
 77: !TEST*/