Actual source code: ex45f.F

petsc-3.8.3 2017-12-09
Report Typos and Errors
  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(MPI_COMM_WORLD,ksp,ierr)
 23:        call DMDACreate2D(MPI_COMM_WORLD, DM_BOUNDARY_NONE,              &
 24:      &    DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,three,three,              &
 25:      &    PETSC_DECIDE,PETSC_DECIDE,one,one, PETSC_NULL_INTEGER,        &
 26:      &    PETSC_NULL_INTEGER, dm, ierr)
 27:        call DMSetFromOptions(dm,ierr)
 28:        call DMSetUp(dm,ierr)
 29:        call KSPSetDM(ksp,dm,ierr)
 30:        call KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,             &
 31:      &                                0,ierr)
 32:        call KSPSetComputeRHS(ksp,ComputeRHS,0,ierr)
 33:        call KSPSetComputeOperators(ksp,ComputeMatrix,                   &
 34:      &      0,ierr)
 35:        call DMDAGetCorners(dm,is,js,PETSC_NULL_INTEGER,iw,jw,             &
 36:      &                     PETSC_NULL_INTEGER,ierr)
 37:        call KSPSetFromOptions(ksp,ierr)
 38:        call KSPSetUp(ksp,ierr)
 39:        call KSPSolve(ksp,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
 40:        call KSPDestroy(ksp,ierr)
 41:        call DMDestroy(dm,ierr)
 42:        call PetscFinalize(ierr)
 43:        end


 46:        subroutine ComputeInitialGuess(ksp,b,ctx,ierr)
 47:        use petscksp
 48:        implicit none
 49:        PetscErrorCode  ierr
 50:        KSP ksp
 51:        PetscInt ctx(*)
 52:        Vec b
 53:        PetscScalar  h

 55:        h=0.0
 56:        call VecSet(b,h,ierr)
 57:        end subroutine

 59:        subroutine ComputeRHS(ksp,b,dummy,ierr)
 60:        use petscksp
 61:        implicit none

 63:        PetscErrorCode  ierr
 64:        KSP ksp
 65:        Vec b
 66:        integer dummy(*)
 67:        PetscScalar  h,Hx,Hy
 68:        PetscInt  mx,my
 69:        DM dm

 71:        call KSPGetDM(ksp,dm,ierr)
 72:        call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,  &
 73:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 74:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 75:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 76:      &                     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
 77:      &                     PETSC_NULL_INTEGER,ierr)

 79:        Hx = 1.0 / real(mx-1)
 80:        Hy = 1.0 / real(my-1)
 81:        h = Hx*Hy
 82:        call VecSet(b,h,ierr)
 83:        end subroutine

 85:       subroutine ComputeMatrix(ksp,A,B,dummy,ierr)
 86:       use petscksp
 87:        implicit none
 88:        PetscErrorCode  ierr
 89:        KSP ksp
 90:        Mat A,B
 91:        integer dummy(*)
 92:        DM dm

 94:       PetscInt    i,j,mx,my,xm
 95:       PetscInt    ym,xs,ys,i1,i5
 96:       PetscScalar  v(5),Hx,Hy
 97:       PetscScalar  HxdHy,HydHx
 98:       MatStencil   row(4),col(4,5)

100:       i1 = 1
101:       i5 = 5
102:       call KSPGetDM(ksp,dm,ierr)
103:       call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,  &
104:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
105:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
106:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
107:      &               PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,             &
108:      &               PETSC_NULL_INTEGER,ierr)

110:       Hx = 1.0 / real(mx-1)
111:       Hy = 1.0 / real(my-1)
112:       HxdHy = Hx/Hy
113:       HydHx = Hy/Hx
114:       call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,            &
115:      &     PETSC_NULL_INTEGER,ierr)
116:       do 10,j=ys,ys+ym-1
117:         do 20,i=xs,xs+xm-1
118:           row(MatStencil_i) = i
119:           row(MatStencil_j) = j
120:           if (i.eq.0 .or. j.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1 ) then
121:             v(1) = 2.0*(HxdHy + HydHx)
122:             call MatSetValuesStencil(B,i1,row,i1,row,v,                 &
123:      &           INSERT_VALUES,ierr)
124:           else
125:             v(1) = -HxdHy
126:             col(MatStencil_i,1) = i
127:             col(MatStencil_j,1) = j-1
128:             v(2) = -HydHx
129:             col(MatStencil_i,2) = i-1
130:             col(MatStencil_j,2) = j
131:             v(3) = 2.0*(HxdHy + HydHx)
132:             col(MatStencil_i,3) = i
133:             col(MatStencil_j,3) = j
134:             v(4) = -HydHx
135:             col(MatStencil_i,4) = i+1
136:             col(MatStencil_j,4) = j
137:             v(5) = -HxdHy
138:             col(MatStencil_i,5) = i
139:             col(MatStencil_j,5) = j+1
140:             call MatSetValuesStencil(B,i1,row,i5,col,v,                 &
141:      &           INSERT_VALUES,ierr)
142:             endif
143:  20      continue
144:  10   continue
145:        call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
146:        call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
147:        if ( A .ne. B) then
148:          call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
149:          call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
150:        endif
151:        end subroutine