Actual source code: ex45f.F
petsc-3.4.5 2014-06-29
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