Actual source code: ex45f.F90

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

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

     14:        one = 1
     15:        three = 3

     17:        call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     18:        if (ierr .ne. 0) then
     19:          print*,'Unable to initialize PETSc'
     20:          stop
     21:        endif
     22:        call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
     23:        call DMDACreate2D(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,three,three,              &
     24:      &                   PETSC_DECIDE,PETSC_DECIDE,one,one, PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, dm, ierr)
     25:        call DMSetFromOptions(dm,ierr)
     26:        call DMSetUp(dm,ierr)
     27:        call KSPSetDM(ksp,dm,ierr)
     28:        call KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,0,ierr)
     29:        call KSPSetComputeRHS(ksp,ComputeRHS,0,ierr)
     30:        call KSPSetComputeOperators(ksp,ComputeMatrix,0,ierr)
     31:        call DMDAGetCorners(dm,is,js,PETSC_NULL_INTEGER,iw,jw,PETSC_NULL_INTEGER,ierr)
     32:        call KSPSetFromOptions(ksp,ierr)
     33:        call KSPSetUp(ksp,ierr)
     34:        call KSPSolve(ksp,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
     35:        call KSPDestroy(ksp,ierr)
     36:        call DMDestroy(dm,ierr)
     37:        call PetscFinalize(ierr)
     38:        end

     41:        subroutine ComputeInitialGuess(ksp,b,ctx,ierr)
     42:        use petscksp
     43:        implicit none
     44:        PetscErrorCode  ierr
     45:        KSP ksp
     46:        PetscInt ctx(*)
     47:        Vec b
     48:        PetscScalar  h

     50:        h=0.0
     51:        call VecSet(b,h,ierr)
     52:        end subroutine

     54:        subroutine ComputeRHS(ksp,b,dummy,ierr)
     55:        use petscksp
     56:        implicit none

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

     66:        call KSPGetDM(ksp,dm,ierr)
     67:        call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
     68:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
     69:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)

     71:        Hx = 1.0 / real(mx-1)
     72:        Hy = 1.0 / real(my-1)
     73:        h = Hx*Hy
     74:        call VecSet(b,h,ierr)
     75:        end subroutine

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

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

     92:       i1 = 1
     93:       i5 = 5
     94:       call KSPGetDM(ksp,dm,ierr)
     95:       call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
     96:      &                 PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
     97:      &                 PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)

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

    139: !/*TEST
    140: !
    141: !   test:
    142: !      nsize: 4
    143: !      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
    144: !
    145: !TEST*/