Actual source code: ex45f.F90

      1: program main
      2: #include <petsc/finclude/petscdmda.h>
      3: #include <petsc/finclude/petscksp.h>
      4:   use petscdm
      5:   use petscdmda
      6:   use petscksp
      7:   implicit none

      9:   PetscInt is, js, iw, jw
     10:   PetscInt one, three
     11:   PetscErrorCode ierr
     12:   KSP ksp
     13:   DM dm
     14:   external ComputeRHS, ComputeMatrix, ComputeInitialGuess

     16:   one = 1
     17:   three = 3

     19:   PetscCallA(PetscInitialize(ierr))
     20:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
     21:   PetscCallA(DMDACreate2D(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, three, three, PETSC_DECIDE, PETSC_DECIDE, one, one, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr))
     22:   PetscCallA(DMSetFromOptions(dm, ierr))
     23:   PetscCallA(DMSetUp(dm, ierr))
     24:   PetscCallA(KSPSetDM(ksp, dm, ierr))
     25:   PetscCallA(KSPSetComputeInitialGuess(ksp, ComputeInitialGuess, 0, ierr))
     26:   PetscCallA(KSPSetComputeRHS(ksp, ComputeRHS, 0, ierr))
     27:   PetscCallA(KSPSetComputeOperators(ksp, ComputeMatrix, 0, ierr))
     28:   PetscCallA(DMDAGetCorners(dm, is, js, PETSC_NULL_INTEGER, iw, jw, PETSC_NULL_INTEGER, ierr))
     29:   PetscCallA(KSPSetFromOptions(ksp, ierr))
     30:   PetscCallA(KSPSetUp(ksp, ierr))
     31:   PetscCallA(KSPSolve(ksp, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))
     32:   PetscCallA(KSPDestroy(ksp, ierr))
     33:   PetscCallA(DMDestroy(dm, ierr))
     34:   PetscCallA(PetscFinalize(ierr))
     35: end

     37: subroutine ComputeInitialGuess(ksp, b, ctx, ierr)
     38:   use petsckspdef
     39:   use petscvec
     40:   implicit none
     41:   PetscErrorCode ierr
     42:   KSP ksp
     43:   PetscInt ctx(*)
     44:   Vec b
     45:   PetscScalar h

     47:   h = 0.0
     48:   PetscCall(VecSet(b, h, ierr))
     49: end subroutine

     51: subroutine ComputeRHS(ksp, b, dummy, ierr)
     52:   use petscksp
     53:   use petscdmda
     54:   use petscvec
     55:   implicit none

     57:   PetscErrorCode ierr
     58:   KSP ksp
     59:   Vec b
     60:   integer dummy(*)
     61:   PetscScalar h, Hx, Hy
     62:   PetscInt mx, my
     63:   DM dm

     65:   PetscCall(KSPGetDM(ksp, dm, ierr))
     66:   PetscCall(DMDAGetInfo(dm, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))

     68:   Hx = 1.0/real(mx - 1)
     69:   Hy = 1.0/real(my - 1)
     70:   h = Hx*Hy
     71:   PetscCall(VecSet(b, h, ierr))
     72: end subroutine

     74: subroutine ComputeMatrix(ksp, A, B, dummy, ierr)
     75:   use petscmat
     76:   use petscksp
     77:   use petscdmda
     78:   implicit none
     79:   PetscErrorCode ierr
     80:   KSP ksp
     81:   Mat A, B
     82:   integer dummy(*)
     83:   DM dm

     85:   PetscInt i, j, mx, my, xm
     86:   PetscInt ym, xs, ys, i1, i5
     87:   PetscScalar v(5), Hx, Hy
     88:   PetscScalar HxdHy, HydHx
     89:   MatStencil row(1), col(5)

     91:   i1 = 1
     92:   i5 = 5
     93:   PetscCall(KSPGetDM(ksp, dm, ierr))
     94:   PetscCall(DMDAGetInfo(dm, PETSC_NULL_INTEGER, mx, my, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMBOUNDARYTYPE, PETSC_NULL_DMDASTENCILTYPE, ierr))

     96:   Hx = 1.0/real(mx - 1)
     97:   Hy = 1.0/real(my - 1)
     98:   HxdHy = Hx/Hy
     99:   HydHx = Hy/Hx
    100:   PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr))
    101:   do 10, j = ys, ys + ym - 1
    102:     do 20, i = xs, xs + xm - 1
    103:       row(1)%i = i
    104:       row(1)%j = j
    105:       if (i == 0 .or. j == 0 .or. i == mx - 1 .or. j == my - 1) then
    106:         v(1) = 2.0*(HxdHy + HydHx)
    107:         PetscCall(MatSetValuesStencil(B, i1, row, i1, row, v, INSERT_VALUES, ierr))
    108:       else
    109:         v(1) = -HxdHy
    110:         col(1)%i = i
    111:         col(1)%j = j - 1
    112:         v(2) = -HydHx
    113:         col(2)%i = i - 1
    114:         col(2)%j = j
    115:         v(3) = 2.0*(HxdHy + HydHx)
    116:         col(3)%i = i
    117:         col(3)%j = j
    118:         v(4) = -HydHx
    119:         col(4)%i = i + 1
    120:         col(4)%j = j
    121:         v(5) = -HxdHy
    122:         col(5)%i = i
    123:         col(5)%j = j + 1
    124:         PetscCall(MatSetValuesStencil(B, i1, row, i5, col, v, INSERT_VALUES, ierr))
    125:       end if
    126: 20    continue
    127: 10    continue
    128:       PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
    129:       PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
    130:       if (A /= B) then
    131:         PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
    132:         PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
    133:       end if
    134:       end subroutine

    136: !/*TEST
    137: !
    138: !   test:
    139: !      nsize: 4
    140: !      args: -ksp_monitor_short -da_refine 5 -pc_type mg -pc_mg_levels 5 -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 2 -mg_levels_pc_type jacobi -ksp_pc_side right
    141: !
    142: !TEST*/