Actual source code: ex100f.F90

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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