Actual source code: ex45f.F
petsc-3.5.4 2015-05-23
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, DM_BOUNDARY_NONE, &
25: & DM_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,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: integer dummy(*)
76: DM dm
78: PetscInt i,j,mx,my,xm
79: PetscInt ym,xs,ys,i1,i5
80: PetscScalar v(5),Hx,Hy
81: PetscScalar HxdHy,HydHx
82: MatStencil row(4),col(4,5)
84: i1 = 1
85: i5 = 5
86: call KSPGetDM(ksp,dm,ierr)
87: call DMDAGetInfo(dm,PETSC_NULL_INTEGER,mx,my,PETSC_NULL_INTEGER, &
88: & PETSC_NULL_INTEGER,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,ierr)
94: Hx = 1.d0 / (mx-1)
95: Hy = 1.d0 / (my-1)
96: HxdHy = Hx/Hy
97: HydHx = Hy/Hx
98: call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
99: & PETSC_NULL_INTEGER,ierr)
100: do 10,j=ys,ys+ym-1
101: do 20,i=xs,xs+xm-1
102: row(MatStencil_i) = i
103: row(MatStencil_j) = j
104: if (i.eq.0 .or. j.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1 ) then
105: v(1) = 2.d0*(HxdHy + HydHx)
106: call MatSetValuesStencil(B,i1,row,i1,row,v, &
107: & INSERT_VALUES,ierr)
108: else
109: v(1) = -HxdHy
110: col(MatStencil_i,1) = i
111: col(MatStencil_j,1) = j-1
112: v(2) = -HydHx
113: col(MatStencil_i,2) = i-1
114: col(MatStencil_j,2) = j
115: v(3) = 2.d0*(HxdHy + HydHx)
116: col(MatStencil_i,3) = i
117: col(MatStencil_j,3) = j
118: v(4) = -HydHx
119: col(MatStencil_i,4) = i+1
120: col(MatStencil_j,4) = j
121: v(5) = -HxdHy
122: col(MatStencil_i,5) = i
123: col(MatStencil_j,5) = j+1
124: call MatSetValuesStencil(B,i1,row,i5,col,v, &
125: & INSERT_VALUES,ierr)
126: endif
127: 20 continue
128: 10 continue
129: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr)
130: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr)
131: if ( A .ne. B) then
132: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
133: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
134: endif
135: end subroutine