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) { 120 PetscCall(PCMGSetX(pcmg,levels - 1 - i,x)); 121 } 122 PetscCall(VecCreateSeq(PETSC_COMM_SELF,N[i],&x)); 123 124 B[levels -1 - i] = x; 125 if (i > 0) { 126 PetscCall(PCMGSetRhs(pcmg,levels - 1 - i,x)); 127 } 128 PetscCall(VecCreateSeq(PETSC_COMM_SELF,N[i],&x)); 129 130 R[levels - 1 - i] = x; 131 132 PetscCall(PCMGSetR(pcmg,levels - 1 - i,x)); 133 } 134 /* create coarse level vectors */ 135 PetscCall(VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x)); 136 PetscCall(PCMGSetX(pcmg,0,x)); X[0] = x; 137 PetscCall(VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x)); 138 PetscCall(PCMGSetRhs(pcmg,0,x)); B[0] = x; 139 140 /* create matrix multiply for finest level */ 141 PetscCall(MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],NULL,&fmat)); 142 PetscCall(MatShellSetOperation(fmat,MATOP_MULT,(void (*)(void))amult)); 143 PetscCall(KSPSetOperators(kspmg,fmat,fmat)); 144 145 PetscCall(CalculateSolution(N[0],&solution)); 146 PetscCall(CalculateRhs(B[levels-1])); 147 PetscCall(VecSet(X[levels-1],0.0)); 148 149 PetscCall(residual((Mat)0,B[levels-1],X[levels-1],R[levels-1])); 150 PetscCall(CalculateError(solution,X[levels-1],R[levels-1],e)); 151 PetscCall(PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",(double)e[0],(double)e[1],(double)e[2])); 152 153 PetscCall(KSPSolve(kspmg,B[levels-1],X[levels-1])); 154 PetscCall(KSPGetIterationNumber(kspmg,&its)); 155 PetscCall(residual((Mat)0,B[levels-1],X[levels-1],R[levels-1])); 156 PetscCall(CalculateError(solution,X[levels-1],R[levels-1],e)); 157 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])); 158 159 PetscCall(PetscFree(N)); 160 PetscCall(VecDestroy(&solution)); 161 162 /* note we have to keep a list of all vectors allocated, this is 163 not ideal, but putting it in MGDestroy is not so good either*/ 164 for (i=0; i<levels; i++) { 165 PetscCall(VecDestroy(&X[i])); 166 PetscCall(VecDestroy(&B[i])); 167 if (i) PetscCall(VecDestroy(&R[i])); 168 } 169 for (i=0; i<levels-1; i++) { 170 PetscCall(MatDestroy(&mat[i])); 171 } 172 PetscCall(MatDestroy(&cmat)); 173 PetscCall(MatDestroy(&fmat)); 174 PetscCall(KSPDestroy(&kspmg)); 175 PetscCall(PetscFinalize()); 176 return 0; 177 } 178 179 /* --------------------------------------------------------------------- */ 180 PetscErrorCode residual(Mat mat,Vec bb,Vec xx,Vec rr) 181 { 182 PetscInt i,n1; 183 PetscScalar *x,*r; 184 const PetscScalar *b; 185 186 PetscFunctionBegin; 187 PetscCall(VecGetSize(bb,&n1)); 188 PetscCall(VecGetArrayRead(bb,&b)); 189 PetscCall(VecGetArray(xx,&x)); 190 PetscCall(VecGetArray(rr,&r)); 191 n1--; 192 r[0] = b[0] + x[1] - 2.0*x[0]; 193 r[n1] = b[n1] + x[n1-1] - 2.0*x[n1]; 194 for (i=1; i<n1; i++) r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i]; 195 PetscCall(VecRestoreArrayRead(bb,&b)); 196 PetscCall(VecRestoreArray(xx,&x)); 197 PetscCall(VecRestoreArray(rr,&r)); 198 PetscFunctionReturn(0); 199 } 200 201 PetscErrorCode amult(Mat mat,Vec xx,Vec yy) 202 { 203 PetscInt i,n1; 204 PetscScalar *y; 205 const PetscScalar *x; 206 207 PetscFunctionBegin; 208 PetscCall(VecGetSize(xx,&n1)); 209 PetscCall(VecGetArrayRead(xx,&x)); 210 PetscCall(VecGetArray(yy,&y)); 211 n1--; 212 y[0] = -x[1] + 2.0*x[0]; 213 y[n1] = -x[n1-1] + 2.0*x[n1]; 214 for (i=1; i<n1; i++) y[i] = -x[i+1] - x[i-1] + 2.0*x[i]; 215 PetscCall(VecRestoreArrayRead(xx,&x)); 216 PetscCall(VecRestoreArray(yy,&y)); 217 PetscFunctionReturn(0); 218 } 219 220 /* --------------------------------------------------------------------- */ 221 PetscErrorCode apply_pc(PC pc,Vec bb,Vec xx) 222 { 223 PetscFunctionBegin; 224 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not implemented"); 225 } 226 227 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) 228 { 229 PetscInt i,n1; 230 PetscScalar *x; 231 const PetscScalar *b; 232 233 PetscFunctionBegin; 234 *its = m; 235 *reason = PCRICHARDSON_CONVERGED_ITS; 236 PetscCall(VecGetSize(bb,&n1)); n1--; 237 PetscCall(VecGetArrayRead(bb,&b)); 238 PetscCall(VecGetArray(xx,&x)); 239 while (m--) { 240 x[0] = .5*(x[1] + b[0]); 241 for (i=1; i<n1; i++) x[i] = .5*(x[i+1] + x[i-1] + b[i]); 242 x[n1] = .5*(x[n1-1] + b[n1]); 243 for (i=n1-1; i>0; i--) x[i] = .5*(x[i+1] + x[i-1] + b[i]); 244 x[0] = .5*(x[1] + b[0]); 245 } 246 PetscCall(VecRestoreArrayRead(bb,&b)); 247 PetscCall(VecRestoreArray(xx,&x)); 248 PetscFunctionReturn(0); 249 } 250 /* --------------------------------------------------------------------- */ 251 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) 252 { 253 PetscInt i,n,n1; 254 PetscScalar *r,*x; 255 const PetscScalar *b; 256 257 PetscFunctionBegin; 258 *its = m; 259 *reason = PCRICHARDSON_CONVERGED_ITS; 260 PetscCall(VecGetSize(bb,&n)); n1 = n - 1; 261 PetscCall(VecGetArrayRead(bb,&b)); 262 PetscCall(VecGetArray(xx,&x)); 263 PetscCall(VecGetArray(w,&r)); 264 265 while (m--) { 266 r[0] = .5*(x[1] + b[0]); 267 for (i=1; i<n1; i++) r[i] = .5*(x[i+1] + x[i-1] + b[i]); 268 r[n1] = .5*(x[n1-1] + b[n1]); 269 for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0; 270 } 271 PetscCall(VecRestoreArrayRead(bb,&b)); 272 PetscCall(VecRestoreArray(xx,&x)); 273 PetscCall(VecRestoreArray(w,&r)); 274 PetscFunctionReturn(0); 275 } 276 /* 277 We know for this application that yy and zz are the same 278 */ 279 /* --------------------------------------------------------------------- */ 280 PetscErrorCode interpolate(Mat mat,Vec xx,Vec yy,Vec zz) 281 { 282 PetscInt i,n,N,i2; 283 PetscScalar *y; 284 const PetscScalar *x; 285 286 PetscFunctionBegin; 287 PetscCall(VecGetSize(yy,&N)); 288 PetscCall(VecGetArrayRead(xx,&x)); 289 PetscCall(VecGetArray(yy,&y)); 290 n = N/2; 291 for (i=0; i<n; i++) { 292 i2 = 2*i; 293 y[i2] += .5*x[i]; 294 y[i2+1] += x[i]; 295 y[i2+2] += .5*x[i]; 296 } 297 PetscCall(VecRestoreArrayRead(xx,&x)); 298 PetscCall(VecRestoreArray(yy,&y)); 299 PetscFunctionReturn(0); 300 } 301 /* --------------------------------------------------------------------- */ 302 PetscErrorCode restrct(Mat mat,Vec rr,Vec bb) 303 { 304 PetscInt i,n,N,i2; 305 PetscScalar *b; 306 const PetscScalar *r; 307 308 PetscFunctionBegin; 309 PetscCall(VecGetSize(rr,&N)); 310 PetscCall(VecGetArrayRead(rr,&r)); 311 PetscCall(VecGetArray(bb,&b)); 312 n = N/2; 313 314 for (i=0; i<n; i++) { 315 i2 = 2*i; 316 b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]); 317 } 318 PetscCall(VecRestoreArrayRead(rr,&r)); 319 PetscCall(VecRestoreArray(bb,&b)); 320 PetscFunctionReturn(0); 321 } 322 /* --------------------------------------------------------------------- */ 323 PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat) 324 { 325 PetscScalar mone = -1.0,two = 2.0; 326 PetscInt i,idx; 327 328 PetscFunctionBegin; 329 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,mat)); 330 331 idx = n-1; 332 PetscCall(MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES)); 333 for (i=0; i<n-1; i++) { 334 PetscCall(MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES)); 335 idx = i+1; 336 PetscCall(MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES)); 337 PetscCall(MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES)); 338 } 339 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 340 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 341 PetscFunctionReturn(0); 342 } 343 /* --------------------------------------------------------------------- */ 344 PetscErrorCode CalculateRhs(Vec u) 345 { 346 PetscInt i,n; 347 PetscReal h; 348 PetscScalar uu; 349 350 PetscFunctionBegin; 351 PetscCall(VecGetSize(u,&n)); 352 h = 1.0/((PetscReal)(n+1)); 353 for (i=0; i<n; i++) { 354 uu = 2.0*h*h; 355 PetscCall(VecSetValues(u,1,&i,&uu,INSERT_VALUES)); 356 } 357 PetscFunctionReturn(0); 358 } 359 /* --------------------------------------------------------------------- */ 360 PetscErrorCode CalculateSolution(PetscInt n,Vec *solution) 361 { 362 PetscInt i; 363 PetscReal h,x = 0.0; 364 PetscScalar uu; 365 366 PetscFunctionBegin; 367 PetscCall(VecCreateSeq(PETSC_COMM_SELF,n,solution)); 368 h = 1.0/((PetscReal)(n+1)); 369 for (i=0; i<n; i++) { 370 x += h; uu = x*(1.-x); 371 PetscCall(VecSetValues(*solution,1,&i,&uu,INSERT_VALUES)); 372 } 373 PetscFunctionReturn(0); 374 } 375 /* --------------------------------------------------------------------- */ 376 PetscErrorCode CalculateError(Vec solution,Vec u,Vec r,PetscReal *e) 377 { 378 PetscFunctionBegin; 379 PetscCall(VecNorm(r,NORM_2,e+2)); 380 PetscCall(VecWAXPY(r,-1.0,u,solution)); 381 PetscCall(VecNorm(r,NORM_2,e)); 382 PetscCall(VecNorm(r,NORM_1,e+1)); 383 PetscFunctionReturn(0); 384 } 385 386 /*TEST 387 388 test: 389 390 TEST*/ 391