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 PetscFunctionReturn(0); 231 } 232 233 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) 234 { 235 PetscInt i,n1; 236 PetscErrorCode ierr; 237 PetscScalar *x; 238 const PetscScalar *b; 239 240 PetscFunctionBegin; 241 *its = m; 242 *reason = PCRICHARDSON_CONVERGED_ITS; 243 ierr = VecGetSize(bb,&n1);CHKERRQ(ierr); n1--; 244 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 245 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 246 while (m--) { 247 x[0] = .5*(x[1] + b[0]); 248 for (i=1; i<n1; i++) x[i] = .5*(x[i+1] + x[i-1] + b[i]); 249 x[n1] = .5*(x[n1-1] + b[n1]); 250 for (i=n1-1; i>0; i--) x[i] = .5*(x[i+1] + x[i-1] + b[i]); 251 x[0] = .5*(x[1] + b[0]); 252 } 253 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 254 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 /* --------------------------------------------------------------------- */ 258 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) 259 { 260 PetscInt i,n,n1; 261 PetscErrorCode ierr; 262 PetscScalar *r,*x; 263 const PetscScalar *b; 264 265 PetscFunctionBegin; 266 *its = m; 267 *reason = PCRICHARDSON_CONVERGED_ITS; 268 ierr = VecGetSize(bb,&n);CHKERRQ(ierr); n1 = n - 1; 269 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 270 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 271 ierr = VecGetArray(w,&r);CHKERRQ(ierr); 272 273 while (m--) { 274 r[0] = .5*(x[1] + b[0]); 275 for (i=1; i<n1; i++) r[i] = .5*(x[i+1] + x[i-1] + b[i]); 276 r[n1] = .5*(x[n1-1] + b[n1]); 277 for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0; 278 } 279 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 280 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 281 ierr = VecRestoreArray(w,&r);CHKERRQ(ierr); 282 PetscFunctionReturn(0); 283 } 284 /* 285 We know for this application that yy and zz are the same 286 */ 287 /* --------------------------------------------------------------------- */ 288 PetscErrorCode interpolate(Mat mat,Vec xx,Vec yy,Vec zz) 289 { 290 PetscInt i,n,N,i2; 291 PetscErrorCode ierr; 292 PetscScalar *y; 293 const PetscScalar *x; 294 295 PetscFunctionBegin; 296 ierr = VecGetSize(yy,&N);CHKERRQ(ierr); 297 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 298 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 299 n = N/2; 300 for (i=0; i<n; i++) { 301 i2 = 2*i; 302 y[i2] += .5*x[i]; 303 y[i2+1] += x[i]; 304 y[i2+2] += .5*x[i]; 305 } 306 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 307 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 308 PetscFunctionReturn(0); 309 } 310 /* --------------------------------------------------------------------- */ 311 PetscErrorCode restrct(Mat mat,Vec rr,Vec bb) 312 { 313 PetscInt i,n,N,i2; 314 PetscErrorCode ierr; 315 PetscScalar *b; 316 const PetscScalar *r; 317 318 PetscFunctionBegin; 319 ierr = VecGetSize(rr,&N);CHKERRQ(ierr); 320 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 321 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 322 n = N/2; 323 324 for (i=0; i<n; i++) { 325 i2 = 2*i; 326 b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]); 327 } 328 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 329 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 330 PetscFunctionReturn(0); 331 } 332 /* --------------------------------------------------------------------- */ 333 PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat) 334 { 335 PetscScalar mone = -1.0,two = 2.0; 336 PetscInt i,idx; 337 PetscErrorCode ierr; 338 339 PetscFunctionBegin; 340 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,mat);CHKERRQ(ierr); 341 342 idx = n-1; 343 ierr = MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);CHKERRQ(ierr); 344 for (i=0; i<n-1; i++) { 345 ierr = MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);CHKERRQ(ierr); 346 idx = i+1; 347 ierr = MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);CHKERRQ(ierr); 348 ierr = MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);CHKERRQ(ierr); 349 } 350 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 351 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 352 PetscFunctionReturn(0); 353 } 354 /* --------------------------------------------------------------------- */ 355 PetscErrorCode CalculateRhs(Vec u) 356 { 357 PetscErrorCode ierr; 358 PetscInt i,n; 359 PetscReal h,x = 0.0; 360 PetscScalar uu; 361 362 PetscFunctionBegin; 363 ierr = VecGetSize(u,&n);CHKERRQ(ierr); 364 h = 1.0/((PetscReal)(n+1)); 365 for (i=0; i<n; i++) { 366 x += h; uu = 2.0*h*h; 367 ierr = VecSetValues(u,1,&i,&uu,INSERT_VALUES);CHKERRQ(ierr); 368 } 369 PetscFunctionReturn(0); 370 } 371 /* --------------------------------------------------------------------- */ 372 PetscErrorCode CalculateSolution(PetscInt n,Vec *solution) 373 { 374 PetscErrorCode ierr; 375 PetscInt i; 376 PetscReal h,x = 0.0; 377 PetscScalar uu; 378 379 PetscFunctionBegin; 380 ierr = VecCreateSeq(PETSC_COMM_SELF,n,solution);CHKERRQ(ierr); 381 h = 1.0/((PetscReal)(n+1)); 382 for (i=0; i<n; i++) { 383 x += h; uu = x*(1.-x); 384 ierr = VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);CHKERRQ(ierr); 385 } 386 PetscFunctionReturn(0); 387 } 388 /* --------------------------------------------------------------------- */ 389 PetscErrorCode CalculateError(Vec solution,Vec u,Vec r,PetscReal *e) 390 { 391 PetscErrorCode ierr; 392 393 PetscFunctionBegin; 394 ierr = VecNorm(r,NORM_2,e+2);CHKERRQ(ierr); 395 ierr = VecWAXPY(r,-1.0,u,solution);CHKERRQ(ierr); 396 ierr = VecNorm(r,NORM_2,e);CHKERRQ(ierr); 397 ierr = VecNorm(r,NORM_1,e+1);CHKERRQ(ierr); 398 PetscFunctionReturn(0); 399 } 400 401 /*TEST 402 403 test: 404 405 TEST*/ 406