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