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