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:       endif

 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*/