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