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 /* --------------------------------------------------------------------- */ 175 PetscErrorCode residual(Mat mat, Vec bb, Vec xx, Vec rr) { 176 PetscInt i, n1; 177 PetscScalar *x, *r; 178 const PetscScalar *b; 179 180 PetscFunctionBegin; 181 PetscCall(VecGetSize(bb, &n1)); 182 PetscCall(VecGetArrayRead(bb, &b)); 183 PetscCall(VecGetArray(xx, &x)); 184 PetscCall(VecGetArray(rr, &r)); 185 n1--; 186 r[0] = b[0] + x[1] - 2.0 * x[0]; 187 r[n1] = b[n1] + x[n1 - 1] - 2.0 * x[n1]; 188 for (i = 1; i < n1; i++) r[i] = b[i] + x[i + 1] + x[i - 1] - 2.0 * x[i]; 189 PetscCall(VecRestoreArrayRead(bb, &b)); 190 PetscCall(VecRestoreArray(xx, &x)); 191 PetscCall(VecRestoreArray(rr, &r)); 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode amult(Mat mat, Vec xx, Vec yy) { 196 PetscInt i, n1; 197 PetscScalar *y; 198 const PetscScalar *x; 199 200 PetscFunctionBegin; 201 PetscCall(VecGetSize(xx, &n1)); 202 PetscCall(VecGetArrayRead(xx, &x)); 203 PetscCall(VecGetArray(yy, &y)); 204 n1--; 205 y[0] = -x[1] + 2.0 * x[0]; 206 y[n1] = -x[n1 - 1] + 2.0 * x[n1]; 207 for (i = 1; i < n1; i++) y[i] = -x[i + 1] - x[i - 1] + 2.0 * x[i]; 208 PetscCall(VecRestoreArrayRead(xx, &x)); 209 PetscCall(VecRestoreArray(yy, &y)); 210 PetscFunctionReturn(0); 211 } 212 213 /* --------------------------------------------------------------------- */ 214 PetscErrorCode apply_pc(PC pc, Vec bb, Vec xx) { 215 PetscFunctionBegin; 216 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented"); 217 } 218 219 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) { 220 PetscInt i, n1; 221 PetscScalar *x; 222 const PetscScalar *b; 223 224 PetscFunctionBegin; 225 *its = m; 226 *reason = PCRICHARDSON_CONVERGED_ITS; 227 PetscCall(VecGetSize(bb, &n1)); 228 n1--; 229 PetscCall(VecGetArrayRead(bb, &b)); 230 PetscCall(VecGetArray(xx, &x)); 231 while (m--) { 232 x[0] = .5 * (x[1] + b[0]); 233 for (i = 1; i < n1; i++) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 234 x[n1] = .5 * (x[n1 - 1] + b[n1]); 235 for (i = n1 - 1; i > 0; i--) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 236 x[0] = .5 * (x[1] + b[0]); 237 } 238 PetscCall(VecRestoreArrayRead(bb, &b)); 239 PetscCall(VecRestoreArray(xx, &x)); 240 PetscFunctionReturn(0); 241 } 242 /* --------------------------------------------------------------------- */ 243 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) { 244 PetscInt i, n, n1; 245 PetscScalar *r, *x; 246 const PetscScalar *b; 247 248 PetscFunctionBegin; 249 *its = m; 250 *reason = PCRICHARDSON_CONVERGED_ITS; 251 PetscCall(VecGetSize(bb, &n)); 252 n1 = n - 1; 253 PetscCall(VecGetArrayRead(bb, &b)); 254 PetscCall(VecGetArray(xx, &x)); 255 PetscCall(VecGetArray(w, &r)); 256 257 while (m--) { 258 r[0] = .5 * (x[1] + b[0]); 259 for (i = 1; i < n1; i++) r[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 260 r[n1] = .5 * (x[n1 - 1] + b[n1]); 261 for (i = 0; i < n; i++) x[i] = (2.0 * r[i] + x[i]) / 3.0; 262 } 263 PetscCall(VecRestoreArrayRead(bb, &b)); 264 PetscCall(VecRestoreArray(xx, &x)); 265 PetscCall(VecRestoreArray(w, &r)); 266 PetscFunctionReturn(0); 267 } 268 /* 269 We know for this application that yy and zz are the same 270 */ 271 /* --------------------------------------------------------------------- */ 272 PetscErrorCode interpolate(Mat mat, Vec xx, Vec yy, Vec zz) { 273 PetscInt i, n, N, i2; 274 PetscScalar *y; 275 const PetscScalar *x; 276 277 PetscFunctionBegin; 278 PetscCall(VecGetSize(yy, &N)); 279 PetscCall(VecGetArrayRead(xx, &x)); 280 PetscCall(VecGetArray(yy, &y)); 281 n = N / 2; 282 for (i = 0; i < n; i++) { 283 i2 = 2 * i; 284 y[i2] += .5 * x[i]; 285 y[i2 + 1] += x[i]; 286 y[i2 + 2] += .5 * x[i]; 287 } 288 PetscCall(VecRestoreArrayRead(xx, &x)); 289 PetscCall(VecRestoreArray(yy, &y)); 290 PetscFunctionReturn(0); 291 } 292 /* --------------------------------------------------------------------- */ 293 PetscErrorCode restrct(Mat mat, Vec rr, Vec bb) { 294 PetscInt i, n, N, i2; 295 PetscScalar *b; 296 const PetscScalar *r; 297 298 PetscFunctionBegin; 299 PetscCall(VecGetSize(rr, &N)); 300 PetscCall(VecGetArrayRead(rr, &r)); 301 PetscCall(VecGetArray(bb, &b)); 302 n = N / 2; 303 304 for (i = 0; i < n; i++) { 305 i2 = 2 * i; 306 b[i] = (r[i2] + 2.0 * r[i2 + 1] + r[i2 + 2]); 307 } 308 PetscCall(VecRestoreArrayRead(rr, &r)); 309 PetscCall(VecRestoreArray(bb, &b)); 310 PetscFunctionReturn(0); 311 } 312 /* --------------------------------------------------------------------- */ 313 PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat) { 314 PetscScalar mone = -1.0, two = 2.0; 315 PetscInt i, idx; 316 317 PetscFunctionBegin; 318 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat)); 319 320 idx = n - 1; 321 PetscCall(MatSetValues(*mat, 1, &idx, 1, &idx, &two, INSERT_VALUES)); 322 for (i = 0; i < n - 1; i++) { 323 PetscCall(MatSetValues(*mat, 1, &i, 1, &i, &two, INSERT_VALUES)); 324 idx = i + 1; 325 PetscCall(MatSetValues(*mat, 1, &idx, 1, &i, &mone, INSERT_VALUES)); 326 PetscCall(MatSetValues(*mat, 1, &i, 1, &idx, &mone, INSERT_VALUES)); 327 } 328 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 329 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 330 PetscFunctionReturn(0); 331 } 332 /* --------------------------------------------------------------------- */ 333 PetscErrorCode CalculateRhs(Vec u) { 334 PetscInt i, n; 335 PetscReal h; 336 PetscScalar uu; 337 338 PetscFunctionBegin; 339 PetscCall(VecGetSize(u, &n)); 340 h = 1.0 / ((PetscReal)(n + 1)); 341 for (i = 0; i < n; i++) { 342 uu = 2.0 * h * h; 343 PetscCall(VecSetValues(u, 1, &i, &uu, INSERT_VALUES)); 344 } 345 PetscFunctionReturn(0); 346 } 347 /* --------------------------------------------------------------------- */ 348 PetscErrorCode CalculateSolution(PetscInt n, Vec *solution) { 349 PetscInt i; 350 PetscReal h, x = 0.0; 351 PetscScalar uu; 352 353 PetscFunctionBegin; 354 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, solution)); 355 h = 1.0 / ((PetscReal)(n + 1)); 356 for (i = 0; i < n; i++) { 357 x += h; 358 uu = x * (1. - x); 359 PetscCall(VecSetValues(*solution, 1, &i, &uu, INSERT_VALUES)); 360 } 361 PetscFunctionReturn(0); 362 } 363 /* --------------------------------------------------------------------- */ 364 PetscErrorCode CalculateError(Vec solution, Vec u, Vec r, PetscReal *e) { 365 PetscFunctionBegin; 366 PetscCall(VecNorm(r, NORM_2, e + 2)); 367 PetscCall(VecWAXPY(r, -1.0, u, solution)); 368 PetscCall(VecNorm(r, NORM_2, e)); 369 PetscCall(VecNorm(r, NORM_1, e + 1)); 370 PetscFunctionReturn(0); 371 } 372 373 /*TEST 374 375 test: 376 377 TEST*/ 378