1 2 static char help[] = "Tests the multigrid code. The input parameters are:\n\ 3 -x N Use a mesh in the x direction of N. \n\ 4 -c N Use N V-cycles. \n\ 5 -l N Use N Levels. \n\ 6 -smooths N Use N pre smooths and N post smooths. \n\ 7 -j Use Jacobi smoother. \n\ 8 -a use additive multigrid \n\ 9 -f use full multigrid (preconditioner variant) \n\ 10 This example also demonstrates matrix-free methods\n\n"; 11 12 /* 13 This is not a good example to understand the use of multigrid with PETSc. 14 */ 15 16 #include <petscksp.h> 17 18 PetscErrorCode residual(Mat, Vec, Vec, Vec); 19 PetscErrorCode gauss_seidel(PC, Vec, Vec, Vec, PetscReal, PetscReal, PetscReal, PetscInt, PetscBool, PetscInt *, PCRichardsonConvergedReason *); 20 PetscErrorCode jacobi_smoother(PC, Vec, Vec, Vec, PetscReal, PetscReal, PetscReal, PetscInt, PetscBool, PetscInt *, PCRichardsonConvergedReason *); 21 PetscErrorCode interpolate(Mat, Vec, Vec, Vec); 22 PetscErrorCode restrct(Mat, Vec, Vec); 23 PetscErrorCode Create1dLaplacian(PetscInt, Mat *); 24 PetscErrorCode CalculateRhs(Vec); 25 PetscErrorCode CalculateError(Vec, Vec, Vec, PetscReal *); 26 PetscErrorCode CalculateSolution(PetscInt, Vec *); 27 PetscErrorCode amult(Mat, Vec, Vec); 28 PetscErrorCode apply_pc(PC, Vec, Vec); 29 30 int main(int Argc, char **Args) { 31 PetscInt x_mesh = 15, levels = 3, cycles = 1, use_jacobi = 0; 32 PetscInt i, smooths = 1, *N, its; 33 PCMGType am = PC_MG_MULTIPLICATIVE; 34 Mat cmat, mat[20], fmat; 35 KSP cksp, ksp[20], kspmg; 36 PetscReal e[3]; /* l_2 error,max error, residual */ 37 const char *shellname; 38 Vec x, solution, X[20], R[20], B[20]; 39 PC pcmg, pc; 40 PetscBool flg; 41 42 PetscCall(PetscInitialize(&Argc, &Args, (char *)0, help)); 43 PetscCall(PetscOptionsGetInt(NULL, NULL, "-x", &x_mesh, NULL)); 44 PetscCall(PetscOptionsGetInt(NULL, NULL, "-l", &levels, NULL)); 45 PetscCall(PetscOptionsGetInt(NULL, NULL, "-c", &cycles, NULL)); 46 PetscCall(PetscOptionsGetInt(NULL, NULL, "-smooths", &smooths, NULL)); 47 PetscCall(PetscOptionsHasName(NULL, NULL, "-a", &flg)); 48 49 if (flg) am = PC_MG_ADDITIVE; 50 PetscCall(PetscOptionsHasName(NULL, NULL, "-f", &flg)); 51 if (flg) am = PC_MG_FULL; 52 PetscCall(PetscOptionsHasName(NULL, NULL, "-j", &flg)); 53 if (flg) use_jacobi = 1; 54 55 PetscCall(PetscMalloc1(levels, &N)); 56 N[0] = x_mesh; 57 for (i = 1; i < levels; i++) { 58 N[i] = N[i - 1] / 2; 59 PetscCheck(N[i] >= 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "Too many levels or N is not large enough"); 60 } 61 62 PetscCall(Create1dLaplacian(N[levels - 1], &cmat)); 63 64 PetscCall(KSPCreate(PETSC_COMM_WORLD, &kspmg)); 65 PetscCall(KSPGetPC(kspmg, &pcmg)); 66 PetscCall(KSPSetFromOptions(kspmg)); 67 PetscCall(PCSetType(pcmg, PCMG)); 68 PetscCall(PCMGSetLevels(pcmg, levels, NULL)); 69 PetscCall(PCMGSetType(pcmg, am)); 70 71 PetscCall(PCMGGetCoarseSolve(pcmg, &cksp)); 72 PetscCall(KSPSetOperators(cksp, cmat, cmat)); 73 PetscCall(KSPGetPC(cksp, &pc)); 74 PetscCall(PCSetType(pc, PCLU)); 75 PetscCall(KSPSetType(cksp, KSPPREONLY)); 76 77 /* zero is finest level */ 78 for (i = 0; i < levels - 1; i++) { 79 Mat dummy; 80 81 PetscCall(PCMGSetResidual(pcmg, levels - 1 - i, residual, NULL)); 82 PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[i + 1], N[i], N[i + 1], N[i], NULL, &mat[i])); 83 PetscCall(MatShellSetOperation(mat[i], MATOP_MULT, (void (*)(void))restrct)); 84 PetscCall(MatShellSetOperation(mat[i], MATOP_MULT_TRANSPOSE_ADD, (void (*)(void))interpolate)); 85 PetscCall(PCMGSetInterpolation(pcmg, levels - 1 - i, mat[i])); 86 PetscCall(PCMGSetRestriction(pcmg, levels - 1 - i, mat[i])); 87 PetscCall(PCMGSetCycleTypeOnLevel(pcmg, levels - 1 - i, (PCMGCycleType)cycles)); 88 89 /* set smoother */ 90 PetscCall(PCMGGetSmoother(pcmg, levels - 1 - i, &ksp[i])); 91 PetscCall(KSPGetPC(ksp[i], &pc)); 92 PetscCall(PCSetType(pc, PCSHELL)); 93 PetscCall(PCShellSetName(pc, "user_precond")); 94 PetscCall(PCShellGetName(pc, &shellname)); 95 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "level=%" PetscInt_FMT ", PCShell name is %s\n", i, shellname)); 96 97 /* this is not used unless different options are passed to the solver */ 98 PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[i], N[i], N[i], N[i], NULL, &dummy)); 99 PetscCall(MatShellSetOperation(dummy, MATOP_MULT, (void (*)(void))amult)); 100 PetscCall(KSPSetOperators(ksp[i], dummy, dummy)); 101 PetscCall(MatDestroy(&dummy)); 102 103 /* 104 We override the matrix passed in by forcing it to use Richardson with 105 a user provided application. This is non-standard and this practice 106 should be avoided. 107 */ 108 PetscCall(PCShellSetApply(pc, apply_pc)); 109 PetscCall(PCShellSetApplyRichardson(pc, gauss_seidel)); 110 if (use_jacobi) PetscCall(PCShellSetApplyRichardson(pc, jacobi_smoother)); 111 PetscCall(KSPSetType(ksp[i], KSPRICHARDSON)); 112 PetscCall(KSPSetInitialGuessNonzero(ksp[i], PETSC_TRUE)); 113 PetscCall(KSPSetTolerances(ksp[i], PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, smooths)); 114 115 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x)); 116 117 X[levels - 1 - i] = x; 118 if (i > 0) PetscCall(PCMGSetX(pcmg, levels - 1 - i, x)); 119 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x)); 120 121 B[levels - 1 - i] = x; 122 if (i > 0) PetscCall(PCMGSetRhs(pcmg, levels - 1 - i, x)); 123 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x)); 124 125 R[levels - 1 - i] = x; 126 127 PetscCall(PCMGSetR(pcmg, levels - 1 - i, x)); 128 } 129 /* create coarse level vectors */ 130 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[levels - 1], &x)); 131 PetscCall(PCMGSetX(pcmg, 0, x)); 132 X[0] = x; 133 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[levels - 1], &x)); 134 PetscCall(PCMGSetRhs(pcmg, 0, x)); 135 B[0] = x; 136 137 /* create matrix multiply for finest level */ 138 PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[0], N[0], N[0], N[0], NULL, &fmat)); 139 PetscCall(MatShellSetOperation(fmat, MATOP_MULT, (void (*)(void))amult)); 140 PetscCall(KSPSetOperators(kspmg, fmat, fmat)); 141 142 PetscCall(CalculateSolution(N[0], &solution)); 143 PetscCall(CalculateRhs(B[levels - 1])); 144 PetscCall(VecSet(X[levels - 1], 0.0)); 145 146 PetscCall(residual((Mat)0, B[levels - 1], X[levels - 1], R[levels - 1])); 147 PetscCall(CalculateError(solution, X[levels - 1], R[levels - 1], e)); 148 PetscCall(PetscPrintf(PETSC_COMM_SELF, "l_2 error %g max error %g resi %g\n", (double)e[0], (double)e[1], (double)e[2])); 149 150 PetscCall(KSPSolve(kspmg, B[levels - 1], X[levels - 1])); 151 PetscCall(KSPGetIterationNumber(kspmg, &its)); 152 PetscCall(residual((Mat)0, B[levels - 1], X[levels - 1], R[levels - 1])); 153 PetscCall(CalculateError(solution, X[levels - 1], R[levels - 1], e)); 154 PetscCall(PetscPrintf(PETSC_COMM_SELF, "its %" PetscInt_FMT " l_2 error %g max error %g resi %g\n", its, (double)e[0], (double)e[1], (double)e[2])); 155 156 PetscCall(PetscFree(N)); 157 PetscCall(VecDestroy(&solution)); 158 159 /* note we have to keep a list of all vectors allocated, this is 160 not ideal, but putting it in MGDestroy is not so good either*/ 161 for (i = 0; i < levels; i++) { 162 PetscCall(VecDestroy(&X[i])); 163 PetscCall(VecDestroy(&B[i])); 164 if (i) PetscCall(VecDestroy(&R[i])); 165 } 166 for (i = 0; i < levels - 1; i++) PetscCall(MatDestroy(&mat[i])); 167 PetscCall(MatDestroy(&cmat)); 168 PetscCall(MatDestroy(&fmat)); 169 PetscCall(KSPDestroy(&kspmg)); 170 PetscCall(PetscFinalize()); 171 return 0; 172 } 173 174 PetscErrorCode residual(Mat mat, Vec bb, Vec xx, Vec rr) { 175 PetscInt i, n1; 176 PetscScalar *x, *r; 177 const PetscScalar *b; 178 179 PetscFunctionBegin; 180 PetscCall(VecGetSize(bb, &n1)); 181 PetscCall(VecGetArrayRead(bb, &b)); 182 PetscCall(VecGetArray(xx, &x)); 183 PetscCall(VecGetArray(rr, &r)); 184 n1--; 185 r[0] = b[0] + x[1] - 2.0 * x[0]; 186 r[n1] = b[n1] + x[n1 - 1] - 2.0 * x[n1]; 187 for (i = 1; i < n1; i++) r[i] = b[i] + x[i + 1] + x[i - 1] - 2.0 * x[i]; 188 PetscCall(VecRestoreArrayRead(bb, &b)); 189 PetscCall(VecRestoreArray(xx, &x)); 190 PetscCall(VecRestoreArray(rr, &r)); 191 PetscFunctionReturn(0); 192 } 193 194 PetscErrorCode amult(Mat mat, Vec xx, Vec yy) { 195 PetscInt i, n1; 196 PetscScalar *y; 197 const PetscScalar *x; 198 199 PetscFunctionBegin; 200 PetscCall(VecGetSize(xx, &n1)); 201 PetscCall(VecGetArrayRead(xx, &x)); 202 PetscCall(VecGetArray(yy, &y)); 203 n1--; 204 y[0] = -x[1] + 2.0 * x[0]; 205 y[n1] = -x[n1 - 1] + 2.0 * x[n1]; 206 for (i = 1; i < n1; i++) y[i] = -x[i + 1] - x[i - 1] + 2.0 * x[i]; 207 PetscCall(VecRestoreArrayRead(xx, &x)); 208 PetscCall(VecRestoreArray(yy, &y)); 209 PetscFunctionReturn(0); 210 } 211 212 PetscErrorCode apply_pc(PC pc, Vec bb, Vec xx) { 213 PetscFunctionBegin; 214 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented"); 215 } 216 217 PetscErrorCode gauss_seidel(PC pc, Vec bb, Vec xx, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt m, PetscBool guesszero, PetscInt *its, PCRichardsonConvergedReason *reason) { 218 PetscInt i, n1; 219 PetscScalar *x; 220 const PetscScalar *b; 221 222 PetscFunctionBegin; 223 *its = m; 224 *reason = PCRICHARDSON_CONVERGED_ITS; 225 PetscCall(VecGetSize(bb, &n1)); 226 n1--; 227 PetscCall(VecGetArrayRead(bb, &b)); 228 PetscCall(VecGetArray(xx, &x)); 229 while (m--) { 230 x[0] = .5 * (x[1] + b[0]); 231 for (i = 1; i < n1; i++) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 232 x[n1] = .5 * (x[n1 - 1] + b[n1]); 233 for (i = n1 - 1; i > 0; i--) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 234 x[0] = .5 * (x[1] + b[0]); 235 } 236 PetscCall(VecRestoreArrayRead(bb, &b)); 237 PetscCall(VecRestoreArray(xx, &x)); 238 PetscFunctionReturn(0); 239 } 240 241 PetscErrorCode jacobi_smoother(PC pc, Vec bb, Vec xx, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt m, PetscBool guesszero, PetscInt *its, PCRichardsonConvergedReason *reason) { 242 PetscInt i, n, n1; 243 PetscScalar *r, *x; 244 const PetscScalar *b; 245 246 PetscFunctionBegin; 247 *its = m; 248 *reason = PCRICHARDSON_CONVERGED_ITS; 249 PetscCall(VecGetSize(bb, &n)); 250 n1 = n - 1; 251 PetscCall(VecGetArrayRead(bb, &b)); 252 PetscCall(VecGetArray(xx, &x)); 253 PetscCall(VecGetArray(w, &r)); 254 255 while (m--) { 256 r[0] = .5 * (x[1] + b[0]); 257 for (i = 1; i < n1; i++) r[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 258 r[n1] = .5 * (x[n1 - 1] + b[n1]); 259 for (i = 0; i < n; i++) x[i] = (2.0 * r[i] + x[i]) / 3.0; 260 } 261 PetscCall(VecRestoreArrayRead(bb, &b)); 262 PetscCall(VecRestoreArray(xx, &x)); 263 PetscCall(VecRestoreArray(w, &r)); 264 PetscFunctionReturn(0); 265 } 266 /* 267 We know for this application that yy and zz are the same 268 */ 269 270 PetscErrorCode interpolate(Mat mat, Vec xx, Vec yy, Vec zz) { 271 PetscInt i, n, N, i2; 272 PetscScalar *y; 273 const PetscScalar *x; 274 275 PetscFunctionBegin; 276 PetscCall(VecGetSize(yy, &N)); 277 PetscCall(VecGetArrayRead(xx, &x)); 278 PetscCall(VecGetArray(yy, &y)); 279 n = N / 2; 280 for (i = 0; i < n; i++) { 281 i2 = 2 * i; 282 y[i2] += .5 * x[i]; 283 y[i2 + 1] += x[i]; 284 y[i2 + 2] += .5 * x[i]; 285 } 286 PetscCall(VecRestoreArrayRead(xx, &x)); 287 PetscCall(VecRestoreArray(yy, &y)); 288 PetscFunctionReturn(0); 289 } 290 291 PetscErrorCode restrct(Mat mat, Vec rr, Vec bb) { 292 PetscInt i, n, N, i2; 293 PetscScalar *b; 294 const PetscScalar *r; 295 296 PetscFunctionBegin; 297 PetscCall(VecGetSize(rr, &N)); 298 PetscCall(VecGetArrayRead(rr, &r)); 299 PetscCall(VecGetArray(bb, &b)); 300 n = N / 2; 301 302 for (i = 0; i < n; i++) { 303 i2 = 2 * i; 304 b[i] = (r[i2] + 2.0 * r[i2 + 1] + r[i2 + 2]); 305 } 306 PetscCall(VecRestoreArrayRead(rr, &r)); 307 PetscCall(VecRestoreArray(bb, &b)); 308 PetscFunctionReturn(0); 309 } 310 311 PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat) { 312 PetscScalar mone = -1.0, two = 2.0; 313 PetscInt i, idx; 314 315 PetscFunctionBegin; 316 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat)); 317 318 idx = n - 1; 319 PetscCall(MatSetValues(*mat, 1, &idx, 1, &idx, &two, INSERT_VALUES)); 320 for (i = 0; i < n - 1; i++) { 321 PetscCall(MatSetValues(*mat, 1, &i, 1, &i, &two, INSERT_VALUES)); 322 idx = i + 1; 323 PetscCall(MatSetValues(*mat, 1, &idx, 1, &i, &mone, INSERT_VALUES)); 324 PetscCall(MatSetValues(*mat, 1, &i, 1, &idx, &mone, INSERT_VALUES)); 325 } 326 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 327 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 328 PetscFunctionReturn(0); 329 } 330 331 PetscErrorCode CalculateRhs(Vec u) { 332 PetscInt i, n; 333 PetscReal h; 334 PetscScalar uu; 335 336 PetscFunctionBegin; 337 PetscCall(VecGetSize(u, &n)); 338 h = 1.0 / ((PetscReal)(n + 1)); 339 for (i = 0; i < n; i++) { 340 uu = 2.0 * h * h; 341 PetscCall(VecSetValues(u, 1, &i, &uu, INSERT_VALUES)); 342 } 343 PetscFunctionReturn(0); 344 } 345 346 PetscErrorCode CalculateSolution(PetscInt n, Vec *solution) { 347 PetscInt i; 348 PetscReal h, x = 0.0; 349 PetscScalar uu; 350 351 PetscFunctionBegin; 352 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, solution)); 353 h = 1.0 / ((PetscReal)(n + 1)); 354 for (i = 0; i < n; i++) { 355 x += h; 356 uu = x * (1. - x); 357 PetscCall(VecSetValues(*solution, 1, &i, &uu, INSERT_VALUES)); 358 } 359 PetscFunctionReturn(0); 360 } 361 362 PetscErrorCode CalculateError(Vec solution, Vec u, Vec r, PetscReal *e) { 363 PetscFunctionBegin; 364 PetscCall(VecNorm(r, NORM_2, e + 2)); 365 PetscCall(VecWAXPY(r, -1.0, u, solution)); 366 PetscCall(VecNorm(r, NORM_2, e)); 367 PetscCall(VecNorm(r, NORM_1, e + 1)); 368 PetscFunctionReturn(0); 369 } 370 371 /*TEST 372 373 test: 374 375 TEST*/ 376