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