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
40: subroutine ComputeInitialGuess(ksp,b,ctx,ierr)
41: use petscksp
42: implicit none
43: PetscErrorCode ierr
44: KSP ksp
45: PetscInt ctx(*)
46: Vec b
47: PetscScalar h
49: h=0.0
50: call VecSet(b,h,ierr)
51: end subroutine
53: subroutine ComputeRHS(ksp,b,dummy,ierr)
54: use petscksp
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: call KSPGetDM(ksp,dm,ierr)
66: call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
67: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
68: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
70: Hx = 1.0 / real(mx-1)
71: Hy = 1.0 / real(my-1)
72: h = Hx*Hy
73: call VecSet(b,h,ierr)
74: end subroutine
76: subroutine ComputeMatrix(ksp,A,B,dummy,ierr)
77: use petscksp
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(4),col(4,5)
91: i1 = 1
92: i5 = 5
93: call KSPGetDM(ksp,dm,ierr)
94: call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
95: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
96: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
98: Hx = 1.0 / real(mx-1)
99: Hy = 1.0 / real(my-1)
100: HxdHy = Hx/Hy
101: HydHx = Hy/Hx
102: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)
103: do 10,j=ys,ys+ym-1
104: do 20,i=xs,xs+xm-1
105: row(MatStencil_i) = i
106: row(MatStencil_j) = j
107: if (i.eq.0 .or. j.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1) then
108: v(1) = 2.0*(HxdHy + HydHx)
109: call MatSetValuesStencil(B,i1,row,i1,row,v,INSERT_VALUES,ierr)
110: else
111: v(1) = -HxdHy
112: col(MatStencil_i,1) = i
113: col(MatStencil_j,1) = j-1
114: v(2) = -HydHx
115: col(MatStencil_i,2) = i-1
116: col(MatStencil_j,2) = j
117: v(3) = 2.0*(HxdHy + HydHx)
118: col(MatStencil_i,3) = i
119: col(MatStencil_j,3) = j
120: v(4) = -HydHx
121: col(MatStencil_i,4) = i+1
122: col(MatStencil_j,4) = j
123: v(5) = -HxdHy
124: col(MatStencil_i,5) = i
125: col(MatStencil_j,5) = j+1
126: call MatSetValuesStencil(B,i1,row,i5,col,v,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
138: !/*TEST
139: !
140: ! test:
141: ! nsize: 4
142: ! 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
143: !
144: !TEST*/