Actual source code: ex45f.F
petsc-3.4.2 2013-07-02

      1:       program main
      2:       implicit none

      4: #include <finclude/petscsys.h>
      5: #include <finclude/petscvec.h>
      6: #include <finclude/petscmat.h>
      7: #include <finclude/petscpc.h>
      8: #include <finclude/petscksp.h>
      9: #include <finclude/petscdm.h>
     10: #include <finclude/petscdmda.h>

     12:        PetscInt is,js,iw,jw
     13:        PetscInt one,three
     14:        PetscErrorCode ierr
     15:        KSP ksp
     16:        DM dm
     17:        external ComputeRHS,ComputeMatrix,ComputeInitialGuess

     19:        one = 1
     20:        three = 3

     22:        call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     23:        call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
     24:        call DMDACreate2D(MPI_COMM_WORLD, DMDA_BOUNDARY_NONE,            &
     25:      &    DMDA_BOUNDARY_NONE, DMDA_STENCIL_STAR,-three,-three,          &
     26:      &    PETSC_DECIDE,PETSC_DECIDE,one,one, PETSC_NULL_INTEGER,        &
     27:      &    PETSC_NULL_INTEGER, dm, ierr)
     28:        call KSPSetDM(ksp,dm,ierr)
     29:        call KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,             &
     30:      &                                PETSC_NULL_OBJECT,ierr)
     31:        call KSPSetComputeRHS(ksp,ComputeRHS,PETSC_NULL_OBJECT,ierr)
     32:        call KSPSetComputeOperators(ksp,ComputeMatrix,                   &
     33:      &      PETSC_NULL_OBJECT,ierr)
     34:        call DMDAGetCorners(dm,is,js,PETSC_NULL_INTEGER,iw,jw,             &
     35:      &                     PETSC_NULL_INTEGER,ierr)
     36:        call KSPSetFromOptions(ksp,ierr)
     37:        call KSPSetUp(ksp,ierr)
     38:        call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
     39:        call KSPDestroy(ksp,ierr)
     40:        call DMDestroy(dm,ierr)
     41:        call PetscFinalize(ierr)
     42:        end

     45:        subroutine ComputeInitialGuess(ksp,b,ctx,ierr)
     46:        implicit none
     47:        PetscErrorCode  ierr
     48:        KSP ksp
     49:        PetscInt ctx(*)
     50:        Vec b
     51:        PetscScalar  h
     52:        h=0.0
     53:        call VecSet(b,h,ierr)
     54:        end subroutine

     56:        subroutine ComputeRHS(ksp,b,dummy,ierr)
     57:        implicit none
     58:        PetscErrorCode  ierr
     59:        KSP ksp
     60:        Vec b
     61:        integer dummy(*)
     62:        PetscScalar  h
     63:        h=1.0
     64:        call VecSet(b,h,ierr)
     65:        end subroutine

     67:        subroutine ComputeMatrix(ksp,A,B,str,dummy,ierr)
     68:        implicit none
     69: #include <finclude/petscsys.h>
     70: #include <finclude/petscvec.h>
     71: #include <finclude/petscmat.h>
     72:        PetscErrorCode  ierr
     73:        KSP ksp
     74:        Mat A,B
     75:        MatStructure str
     76:        integer dummy(*)
     77:        DM dm

     79:       PetscInt    i,j,mx,my,xm
     80:       PetscInt    ym,xs,ys,i1,i5
     81:       PetscScalar  v(5),Hx,Hy
     82:       PetscScalar  HxdHy,HydHx
     83:       MatStencil   row(4),col(4,5)

     85:       i1 = 1
     86:       i5 = 5
     87:       call KSPGetDM(ksp,dm,ierr)
     88:       call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,  &
     89:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
     90:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
     91:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
     92:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
     93:      &               PETSC_NULL_INTEGER,ierr)

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