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