1 static char help[] = "Tests the multigrid code. The input parameters are:\n\ 2 -x N Use a mesh in the x direction of N. \n\ 3 -c N Use N V-cycles. \n\ 4 -l N Use N Levels. \n\ 5 -smooths N Use N pre smooths and N post smooths. \n\ 6 -j Use Jacobi smoother. \n\ 7 -a use additive multigrid \n\ 8 -f use full multigrid (preconditioner variant) \n\ 9 This example also demonstrates matrix-free methods\n\n"; 10 11 /* 12 This is not a good example to understand the use of multigrid with PETSc. 13 */ 14 15 #include <petscksp.h> 16 17 PetscErrorCode residual(Mat, Vec, Vec, Vec); 18 PetscErrorCode gauss_seidel(PC, Vec, Vec, Vec, PetscReal, PetscReal, PetscReal, PetscInt, PetscBool, PetscInt *, PCRichardsonConvergedReason *); 19 PetscErrorCode jacobi_smoother(PC, Vec, Vec, Vec, PetscReal, PetscReal, PetscReal, PetscInt, PetscBool, PetscInt *, PCRichardsonConvergedReason *); 20 PetscErrorCode interpolate(Mat, Vec, Vec, Vec); 21 PetscErrorCode restrct(Mat, Vec, Vec); 22 PetscErrorCode Create1dLaplacian(PetscInt, Mat *); 23 PetscErrorCode CalculateRhs(Vec); 24 PetscErrorCode CalculateError(Vec, Vec, Vec, PetscReal *); 25 PetscErrorCode CalculateSolution(PetscInt, Vec *); 26 PetscErrorCode amult(Mat, Vec, Vec); 27 PetscErrorCode apply_pc(PC, Vec, Vec); 28 29 int main(int Argc, char **Args) 30 { 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 { 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(PETSC_SUCCESS); 193 } 194 195 PetscErrorCode amult(Mat mat, Vec xx, Vec yy) 196 { 197 PetscInt i, n1; 198 PetscScalar *y; 199 const PetscScalar *x; 200 201 PetscFunctionBegin; 202 PetscCall(VecGetSize(xx, &n1)); 203 PetscCall(VecGetArrayRead(xx, &x)); 204 PetscCall(VecGetArray(yy, &y)); 205 n1--; 206 y[0] = -x[1] + 2.0 * x[0]; 207 y[n1] = -x[n1 - 1] + 2.0 * x[n1]; 208 for (i = 1; i < n1; i++) y[i] = -x[i + 1] - x[i - 1] + 2.0 * x[i]; 209 PetscCall(VecRestoreArrayRead(xx, &x)); 210 PetscCall(VecRestoreArray(yy, &y)); 211 PetscFunctionReturn(PETSC_SUCCESS); 212 } 213 214 PetscErrorCode apply_pc(PC pc, Vec bb, Vec xx) 215 { 216 PetscFunctionBegin; 217 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented"); 218 } 219 220 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) 221 { 222 PetscInt i, n1; 223 PetscScalar *x; 224 const PetscScalar *b; 225 226 PetscFunctionBegin; 227 *its = m; 228 *reason = PCRICHARDSON_CONVERGED_ITS; 229 PetscCall(VecGetSize(bb, &n1)); 230 n1--; 231 PetscCall(VecGetArrayRead(bb, &b)); 232 PetscCall(VecGetArray(xx, &x)); 233 while (m--) { 234 x[0] = .5 * (x[1] + b[0]); 235 for (i = 1; i < n1; i++) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 236 x[n1] = .5 * (x[n1 - 1] + b[n1]); 237 for (i = n1 - 1; i > 0; i--) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 238 x[0] = .5 * (x[1] + b[0]); 239 } 240 PetscCall(VecRestoreArrayRead(bb, &b)); 241 PetscCall(VecRestoreArray(xx, &x)); 242 PetscFunctionReturn(PETSC_SUCCESS); 243 } 244 245 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) 246 { 247 PetscInt i, n, n1; 248 PetscScalar *r, *x; 249 const PetscScalar *b; 250 251 PetscFunctionBegin; 252 *its = m; 253 *reason = PCRICHARDSON_CONVERGED_ITS; 254 PetscCall(VecGetSize(bb, &n)); 255 n1 = n - 1; 256 PetscCall(VecGetArrayRead(bb, &b)); 257 PetscCall(VecGetArray(xx, &x)); 258 PetscCall(VecGetArray(w, &r)); 259 260 while (m--) { 261 r[0] = .5 * (x[1] + b[0]); 262 for (i = 1; i < n1; i++) r[i] = .5 * (x[i + 1] + x[i - 1] + b[i]); 263 r[n1] = .5 * (x[n1 - 1] + b[n1]); 264 for (i = 0; i < n; i++) x[i] = (2.0 * r[i] + x[i]) / 3.0; 265 } 266 PetscCall(VecRestoreArrayRead(bb, &b)); 267 PetscCall(VecRestoreArray(xx, &x)); 268 PetscCall(VecRestoreArray(w, &r)); 269 PetscFunctionReturn(PETSC_SUCCESS); 270 } 271 /* 272 We know for this application that yy and zz are the same 273 */ 274 275 PetscErrorCode interpolate(Mat mat, Vec xx, Vec yy, Vec zz) 276 { 277 PetscInt i, n, N, i2; 278 PetscScalar *y; 279 const PetscScalar *x; 280 281 PetscFunctionBegin; 282 PetscCall(VecGetSize(yy, &N)); 283 PetscCall(VecGetArrayRead(xx, &x)); 284 PetscCall(VecGetArray(yy, &y)); 285 n = N / 2; 286 for (i = 0; i < n; i++) { 287 i2 = 2 * i; 288 y[i2] += .5 * x[i]; 289 y[i2 + 1] += x[i]; 290 y[i2 + 2] += .5 * x[i]; 291 } 292 PetscCall(VecRestoreArrayRead(xx, &x)); 293 PetscCall(VecRestoreArray(yy, &y)); 294 PetscFunctionReturn(PETSC_SUCCESS); 295 } 296 297 PetscErrorCode restrct(Mat mat, Vec rr, Vec bb) 298 { 299 PetscInt i, n, N, i2; 300 PetscScalar *b; 301 const PetscScalar *r; 302 303 PetscFunctionBegin; 304 PetscCall(VecGetSize(rr, &N)); 305 PetscCall(VecGetArrayRead(rr, &r)); 306 PetscCall(VecGetArray(bb, &b)); 307 n = N / 2; 308 309 for (i = 0; i < n; i++) { 310 i2 = 2 * i; 311 b[i] = (r[i2] + 2.0 * r[i2 + 1] + r[i2 + 2]); 312 } 313 PetscCall(VecRestoreArrayRead(rr, &r)); 314 PetscCall(VecRestoreArray(bb, &b)); 315 PetscFunctionReturn(PETSC_SUCCESS); 316 } 317 318 PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat) 319 { 320 PetscScalar mone = -1.0, two = 2.0; 321 PetscInt i, idx; 322 323 PetscFunctionBegin; 324 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat)); 325 326 idx = n - 1; 327 PetscCall(MatSetValues(*mat, 1, &idx, 1, &idx, &two, INSERT_VALUES)); 328 for (i = 0; i < n - 1; i++) { 329 PetscCall(MatSetValues(*mat, 1, &i, 1, &i, &two, INSERT_VALUES)); 330 idx = i + 1; 331 PetscCall(MatSetValues(*mat, 1, &idx, 1, &i, &mone, INSERT_VALUES)); 332 PetscCall(MatSetValues(*mat, 1, &i, 1, &idx, &mone, INSERT_VALUES)); 333 } 334 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 335 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 336 PetscFunctionReturn(PETSC_SUCCESS); 337 } 338 339 PetscErrorCode CalculateRhs(Vec u) 340 { 341 PetscInt i, n; 342 PetscReal h; 343 PetscScalar uu; 344 345 PetscFunctionBegin; 346 PetscCall(VecGetSize(u, &n)); 347 h = 1.0 / ((PetscReal)(n + 1)); 348 for (i = 0; i < n; i++) { 349 uu = 2.0 * h * h; 350 PetscCall(VecSetValues(u, 1, &i, &uu, INSERT_VALUES)); 351 } 352 PetscFunctionReturn(PETSC_SUCCESS); 353 } 354 355 PetscErrorCode CalculateSolution(PetscInt n, Vec *solution) 356 { 357 PetscInt i; 358 PetscReal h, x = 0.0; 359 PetscScalar uu; 360 361 PetscFunctionBegin; 362 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, solution)); 363 h = 1.0 / ((PetscReal)(n + 1)); 364 for (i = 0; i < n; i++) { 365 x += h; 366 uu = x * (1. - x); 367 PetscCall(VecSetValues(*solution, 1, &i, &uu, INSERT_VALUES)); 368 } 369 PetscFunctionReturn(PETSC_SUCCESS); 370 } 371 372 PetscErrorCode CalculateError(Vec solution, Vec u, Vec r, PetscReal *e) 373 { 374 PetscFunctionBegin; 375 PetscCall(VecNorm(r, NORM_2, e + 2)); 376 PetscCall(VecWAXPY(r, -1.0, u, solution)); 377 PetscCall(VecNorm(r, NORM_2, e)); 378 PetscCall(VecNorm(r, NORM_1, e + 1)); 379 PetscFunctionReturn(PETSC_SUCCESS); 380 } 381 382 /*TEST 383 384 test: 385 386 TEST*/ 387