Actual source code: ex8f.F90
1: !
2: ! Tests PCMGSetResidual() and PCMGSetResidualTranspose(), including on a MATSHELL with a Fortran MatMult()
3: !
4: ! -----------------------------------------------------------------------
5: #include <petsc/finclude/petscksp.h>
6: module ex8fmodule
7: use petscksp
8: implicit none
9: PetscInt :: nmult = 0, nresidual = 0, nresidualt = 0
11: contains
12: subroutine MyMult(S, x, y, ierr)
13: Mat S
14: Vec x, y
15: PetscErrorCode, intent(out) :: ierr
16: nmult = nmult + 1
17: PetscCall(VecCopy(x, y, ierr))
18: end
20: subroutine MyResidual(A, b, x, r, ierr)
21: Mat A
22: Vec b, x, r
23: PetscErrorCode, intent(out) :: ierr
24: PetscScalar, parameter :: mone = -1.0
25: nresidual = nresidual + 1
26: PetscCall(MatMult(A, x, r, ierr))
27: PetscCall(VecAYPX(r, mone, b, ierr))
28: end
30: subroutine MyResidualTranspose(A, b, x, r, ierr)
31: Mat A
32: Vec b, x, r
33: PetscErrorCode, intent(out) :: ierr
34: PetscScalar, parameter :: mone = -1.0
35: nresidualt = nresidualt + 1
36: PetscCall(MatMultTranspose(A, x, r, ierr))
37: PetscCall(VecAYPX(r, mone, b, ierr))
38: end
40: end module ex8fmodule
42: program main
43: use petscksp
44: use ex8fmodule
45: implicit none
47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: ! Variable declarations
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: !
51: ! Variables:
52: ! ksp - linear solver context
53: ! x, b, u - approx solution, right-hand side, exact solution vectors
54: ! A - matrix that defines linear system
55: ! its - iterations for convergence
56: ! norm - norm of error in solution
57: ! rctx - random number context
58: !
60: Mat A, S, P
61: Vec x, b, u
62: PC pc
63: PetscInt, parameter :: n = 6, dim = n**2, dimc = dim/2
64: PetscInt i, j, jj, ii, istart, iend
65: PetscErrorCode ierr
66: PetscScalar v
67: PetscScalar, parameter :: pfive = .5
68: KSP ksp, sksp
69: PC spc
71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: ! Beginning of program
73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: PetscCallA(PetscInitialize(ierr))
77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: ! Compute the matrix and right-hand-side vector that define
79: ! the linear system, Ax = b.
80: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: ! Create parallel matrix, specifying only its global dimensions.
83: ! When using MatCreate(), the matrix format can be specified at
84: ! runtime. Also, the parallel partitioning of the matrix is
85: ! determined by PETSc at runtime.
87: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
88: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim, ierr))
89: PetscCallA(MatSetFromOptions(A, ierr))
90: PetscCallA(MatSetUp(A, ierr))
92: ! Currently, all PETSc parallel matrix formats are partitioned by
93: ! contiguous chunks of rows across the processors. Determine which
94: ! rows of the matrix are locally owned.
96: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
98: ! Set matrix elements in parallel.
99: ! - Each processor needs to insert only elements that it owns
100: ! locally (but any non-local elements will be sent to the
101: ! appropriate processor during matrix assembly).
102: ! - Always specify global rows and columns of matrix entries.
104: do II = Istart, Iend - 1
105: v = -1.0
106: i = II/n
107: j = II - i*n
108: if (i > 0) then
109: JJ = II - n
110: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
111: end if
112: if (i < n - 1) then
113: JJ = II + n
114: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
115: end if
116: if (j > 0) then
117: JJ = II - 1
118: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
119: end if
120: if (j < n - 1) then
121: JJ = II + 1
122: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], ADD_VALUES, ierr))
123: end if
124: v = 4.0
125: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [II], [v], ADD_VALUES, ierr))
126: end do
128: ! Assemble matrix, using the 2-step process:
129: ! MatAssemblyBegin(), MatAssemblyEnd()
130: ! Computations can be done while messages are in transition
131: ! by placing code between these two statements.
133: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
134: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
136: ! Create parallel vectors.
137: ! - Here, the parallel partitioning of the vector is determined by
138: ! PETSc at runtime. We could also specify the local dimensions
139: ! if desired.
140: ! - Note: We form 1 vector from scratch and then duplicate as needed.
142: PetscCallA(VecCreate(PETSC_COMM_WORLD, u, ierr))
143: PetscCallA(VecSetSizes(u, PETSC_DECIDE, dim, ierr))
144: PetscCallA(VecSetFromOptions(u, ierr))
145: PetscCallA(VecDuplicate(u, b, ierr))
146: PetscCallA(VecDuplicate(b, x, ierr))
148: ! Set exact solution; then compute right-hand-side vector.
150: PetscCallA(VecSet(u, pfive, ierr))
151: PetscCallA(MatMult(A, u, b, ierr))
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: ! Create the linear solver and set various options
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! Create linear solver context
159: PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
160: PetscCallA(KSPGetPC(ksp, pc, ierr))
161: PetscCallA(PCSetType(pc, PCMG, ierr))
162: PetscCallA(PCMGSetLevels(pc, 2_PETSC_INT_KIND, PETSC_NULL_MPI_COMM, ierr))
163: PetscCallA(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH, ierr))
165: ! Piecewise constant interpolation from pairs of fine points
167: PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, dim, dimc, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, P, ierr))
168: PetscCallA(MatGetOwnershipRange(P, Istart, Iend, ierr))
169: v = 1.0
170: do II = Istart, Iend - 1
171: JJ = II/2
172: PetscCallA(MatSetValues(P, 1_PETSC_INT_KIND, [II], 1_PETSC_INT_KIND, [JJ], [v], INSERT_VALUES, ierr))
173: end do
174: PetscCallA(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY, ierr))
175: PetscCallA(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY, ierr))
176: PetscCallA(PCMGSetInterpolation(pc, 1_PETSC_INT_KIND, P, ierr))
178: ! The residual routines set on a MATSHELL must not replace its Fortran MatMult
180: PetscCallA(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, dim, dim, PETSC_NULL_INTEGER, S, ierr))
181: PetscCallA(MatShellSetOperation(S, MATOP_MULT, MyMult, ierr))
182: PetscCallA(PCMGSetResidual(pc, 1_PETSC_INT_KIND, MyResidual, S, ierr))
183: PetscCallA(PCMGSetResidualTranspose(pc, 1_PETSC_INT_KIND, MyResidualTranspose, S, ierr))
184: PetscCallA(MatMult(S, u, x, ierr))
185: PetscCheckA(nmult == 1, PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'The MatMult of the MATSHELL was not called')
187: PetscCallA(PCMGSetResidual(pc, 1_PETSC_INT_KIND, MyResidual, A, ierr))
188: PetscCallA(PCMGSetResidualTranspose(pc, 1_PETSC_INT_KIND, MyResidualTranspose, A, ierr))
190: ! Smoothers whose transpose can be applied
192: PetscCallA(PCMGGetSmoother(pc, 1_PETSC_INT_KIND, sksp, ierr))
193: PetscCallA(KSPSetType(sksp, KSPRICHARDSON, ierr))
194: PetscCallA(KSPGetPC(sksp, spc, ierr))
195: PetscCallA(PCSetType(spc, PCJACOBI, ierr))
197: ! Set operators. Here the matrix that defines the linear system
198: ! also serves as the matrix used to construct the preconditioner.
200: PetscCallA(KSPSetOperators(ksp, A, A, ierr))
201: PetscCallA(KSPSetUp(ksp, ierr))
203: ! Both residual routines must be called by the cycles
205: PetscCallA(PCApply(pc, b, x, ierr))
206: PetscCheckA(nresidual > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'The residual routine was not called')
207: PetscCheckA(nresidualt == 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'The transposed residual routine was called by PCApply()')
208: PetscCallA(PCApplyTranspose(pc, b, x, ierr))
209: PetscCheckA(nresidualt > 0, PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'The transposed residual routine was not called')
211: PetscCallA(MatDestroy(P, ierr))
212: PetscCallA(KSPDestroy(ksp, ierr))
213: PetscCallA(VecDestroy(u, ierr))
214: PetscCallA(VecDestroy(x, ierr))
215: PetscCallA(VecDestroy(b, ierr))
216: PetscCallA(MatDestroy(A, ierr))
217: PetscCallA(MatDestroy(S, ierr))
219: PetscCallA(PetscFinalize(ierr))
220: end
222: !/*TEST
223: !
224: ! test:
225: ! nsize: 1
226: ! output_file: output/empty.out
227: !TEST*/