1 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 2 #include <petscsf.h> 3 4 PetscFunctionList PCGAMGClassicalProlongatorList = NULL; 5 PetscBool PCGAMGClassicalPackageInitialized = PETSC_FALSE; 6 7 typedef struct { 8 PetscReal interp_threshold; /* interpolation threshold */ 9 char prolongtype[256]; 10 PetscInt nsmooths; /* number of jacobi smoothings on the prolongator */ 11 } PC_GAMG_Classical; 12 13 /*@C 14 PCGAMGClassicalSetType - Sets the type of classical interpolation to use 15 16 Collective on PC 17 18 Input Parameters: 19 . pc - the preconditioner context 20 21 Options Database Key: 22 . -pc_gamg_classical_type <direct,standard> - set type of Classical AMG prolongation 23 24 Level: intermediate 25 26 .seealso: () 27 @*/ 28 PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type) 29 { 30 PetscFunctionBegin; 31 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 32 CHKERRQ(PetscTryMethod(pc,"PCGAMGClassicalSetType_C",(PC,PCGAMGClassicalType),(pc,type))); 33 PetscFunctionReturn(0); 34 } 35 36 /*@C 37 PCGAMGClassicalGetType - Gets the type of classical interpolation to use 38 39 Collective on PC 40 41 Input Parameter: 42 . pc - the preconditioner context 43 44 Output Parameter: 45 . type - the type used 46 47 Level: intermediate 48 49 .seealso: () 50 @*/ 51 PetscErrorCode PCGAMGClassicalGetType(PC pc, PCGAMGClassicalType *type) 52 { 53 PetscFunctionBegin; 54 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 55 CHKERRQ(PetscUseMethod(pc,"PCGAMGClassicalGetType_C",(PC,PCGAMGClassicalType*),(pc,type))); 56 PetscFunctionReturn(0); 57 } 58 59 static PetscErrorCode PCGAMGClassicalSetType_GAMG(PC pc, PCGAMGClassicalType type) 60 { 61 PC_MG *mg = (PC_MG*)pc->data; 62 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 63 PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx; 64 65 PetscFunctionBegin; 66 CHKERRQ(PetscStrcpy(cls->prolongtype,type)); 67 PetscFunctionReturn(0); 68 } 69 70 static PetscErrorCode PCGAMGClassicalGetType_GAMG(PC pc, PCGAMGClassicalType *type) 71 { 72 PC_MG *mg = (PC_MG*)pc->data; 73 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 74 PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx; 75 76 PetscFunctionBegin; 77 *type = cls->prolongtype; 78 PetscFunctionReturn(0); 79 } 80 81 PetscErrorCode PCGAMGGraph_Classical(PC pc,Mat A,Mat *G) 82 { 83 PetscInt s,f,n,idx,lidx,gidx; 84 PetscInt r,c,ncols; 85 const PetscInt *rcol; 86 const PetscScalar *rval; 87 PetscInt *gcol; 88 PetscScalar *gval; 89 PetscReal rmax; 90 PetscInt cmax = 0; 91 PC_MG *mg = (PC_MG *)pc->data; 92 PC_GAMG *gamg = (PC_GAMG *)mg->innerctx; 93 PetscInt *gsparse,*lsparse; 94 PetscScalar *Amax; 95 MatType mtype; 96 97 PetscFunctionBegin; 98 CHKERRQ(MatGetOwnershipRange(A,&s,&f)); 99 n=f-s; 100 CHKERRQ(PetscMalloc3(n,&lsparse,n,&gsparse,n,&Amax)); 101 102 for (r = 0;r < n;r++) { 103 lsparse[r] = 0; 104 gsparse[r] = 0; 105 } 106 107 for (r = s;r < f;r++) { 108 /* determine the maximum off-diagonal in each row */ 109 rmax = 0.; 110 CHKERRQ(MatGetRow(A,r,&ncols,&rcol,&rval)); 111 for (c = 0; c < ncols; c++) { 112 if (PetscRealPart(-rval[c]) > rmax && rcol[c] != r) { 113 rmax = PetscRealPart(-rval[c]); 114 } 115 } 116 Amax[r-s] = rmax; 117 if (ncols > cmax) cmax = ncols; 118 lidx = 0; 119 gidx = 0; 120 /* create the local and global sparsity patterns */ 121 for (c = 0; c < ncols; c++) { 122 if (PetscRealPart(-rval[c]) > gamg->threshold[0]*PetscRealPart(Amax[r-s]) || rcol[c] == r) { 123 if (rcol[c] < f && rcol[c] >= s) { 124 lidx++; 125 } else { 126 gidx++; 127 } 128 } 129 } 130 CHKERRQ(MatRestoreRow(A,r,&ncols,&rcol,&rval)); 131 lsparse[r-s] = lidx; 132 gsparse[r-s] = gidx; 133 } 134 CHKERRQ(PetscMalloc2(cmax,&gval,cmax,&gcol)); 135 136 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),G)); 137 CHKERRQ(MatGetType(A,&mtype)); 138 CHKERRQ(MatSetType(*G,mtype)); 139 CHKERRQ(MatSetSizes(*G,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 140 CHKERRQ(MatMPIAIJSetPreallocation(*G,0,lsparse,0,gsparse)); 141 CHKERRQ(MatSeqAIJSetPreallocation(*G,0,lsparse)); 142 for (r = s;r < f;r++) { 143 CHKERRQ(MatGetRow(A,r,&ncols,&rcol,&rval)); 144 idx = 0; 145 for (c = 0; c < ncols; c++) { 146 /* classical strength of connection */ 147 if (PetscRealPart(-rval[c]) > gamg->threshold[0]*PetscRealPart(Amax[r-s]) || rcol[c] == r) { 148 gcol[idx] = rcol[c]; 149 gval[idx] = rval[c]; 150 idx++; 151 } 152 } 153 CHKERRQ(MatSetValues(*G,1,&r,idx,gcol,gval,INSERT_VALUES)); 154 CHKERRQ(MatRestoreRow(A,r,&ncols,&rcol,&rval)); 155 } 156 CHKERRQ(MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY)); 157 CHKERRQ(MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY)); 158 159 CHKERRQ(PetscFree2(gval,gcol)); 160 CHKERRQ(PetscFree3(lsparse,gsparse,Amax)); 161 PetscFunctionReturn(0); 162 } 163 164 PetscErrorCode PCGAMGCoarsen_Classical(PC pc,Mat *G,PetscCoarsenData **agg_lists) 165 { 166 MatCoarsen crs; 167 MPI_Comm fcomm = ((PetscObject)pc)->comm; 168 169 PetscFunctionBegin; 170 PetscCheck(G,fcomm,PETSC_ERR_ARG_WRONGSTATE,"Must set Graph in PC in PCGAMG before coarsening"); 171 172 CHKERRQ(MatCoarsenCreate(fcomm,&crs)); 173 CHKERRQ(MatCoarsenSetFromOptions(crs)); 174 CHKERRQ(MatCoarsenSetAdjacency(crs,*G)); 175 CHKERRQ(MatCoarsenSetStrictAggs(crs,PETSC_TRUE)); 176 CHKERRQ(MatCoarsenApply(crs)); 177 CHKERRQ(MatCoarsenGetData(crs,agg_lists)); 178 CHKERRQ(MatCoarsenDestroy(&crs)); 179 PetscFunctionReturn(0); 180 } 181 182 PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P) 183 { 184 PC_MG *mg = (PC_MG*)pc->data; 185 PC_GAMG *gamg = (PC_GAMG*)mg->innerctx; 186 PetscBool iscoarse,isMPIAIJ,isSEQAIJ; 187 PetscInt fn,cn,fs,fe,cs,ce,i,j,ncols,col,row_f,row_c,cmax=0,idx,noff; 188 PetscInt *lcid,*gcid,*lsparse,*gsparse,*colmap,*pcols; 189 const PetscInt *rcol; 190 PetscReal *Amax_pos,*Amax_neg; 191 PetscScalar g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta,pij; 192 PetscScalar *pvals; 193 const PetscScalar *rval; 194 Mat lA,gA=NULL; 195 MatType mtype; 196 Vec C,lvec; 197 PetscLayout clayout; 198 PetscSF sf; 199 Mat_MPIAIJ *mpiaij; 200 201 PetscFunctionBegin; 202 CHKERRQ(MatGetOwnershipRange(A,&fs,&fe)); 203 fn = fe-fs; 204 CHKERRQ(PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ)); 205 CHKERRQ(PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSEQAIJ)); 206 PetscCheckFalse(!isMPIAIJ && !isSEQAIJ,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Classical AMG requires MPIAIJ matrix"); 207 if (isMPIAIJ) { 208 mpiaij = (Mat_MPIAIJ*)A->data; 209 lA = mpiaij->A; 210 gA = mpiaij->B; 211 lvec = mpiaij->lvec; 212 CHKERRQ(VecGetSize(lvec,&noff)); 213 colmap = mpiaij->garray; 214 CHKERRQ(MatGetLayouts(A,NULL,&clayout)); 215 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject)A),&sf)); 216 CHKERRQ(PetscSFSetGraphLayout(sf,clayout,noff,NULL,PETSC_COPY_VALUES,colmap)); 217 CHKERRQ(PetscMalloc1(noff,&gcid)); 218 } else { 219 lA = A; 220 } 221 CHKERRQ(PetscMalloc5(fn,&lsparse,fn,&gsparse,fn,&lcid,fn,&Amax_pos,fn,&Amax_neg)); 222 223 /* count the number of coarse unknowns */ 224 cn = 0; 225 for (i=0;i<fn;i++) { 226 /* filter out singletons */ 227 CHKERRQ(PetscCDEmptyAt(agg_lists,i,&iscoarse)); 228 lcid[i] = -1; 229 if (!iscoarse) { 230 cn++; 231 } 232 } 233 234 /* create the coarse vector */ 235 CHKERRQ(VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&C)); 236 CHKERRQ(VecGetOwnershipRange(C,&cs,&ce)); 237 238 cn = 0; 239 for (i=0;i<fn;i++) { 240 CHKERRQ(PetscCDEmptyAt(agg_lists,i,&iscoarse)); 241 if (!iscoarse) { 242 lcid[i] = cs+cn; 243 cn++; 244 } else { 245 lcid[i] = -1; 246 } 247 } 248 249 if (gA) { 250 CHKERRQ(PetscSFBcastBegin(sf,MPIU_INT,lcid,gcid,MPI_REPLACE)); 251 CHKERRQ(PetscSFBcastEnd(sf,MPIU_INT,lcid,gcid,MPI_REPLACE)); 252 } 253 254 /* determine the largest off-diagonal entries in each row */ 255 for (i=fs;i<fe;i++) { 256 Amax_pos[i-fs] = 0.; 257 Amax_neg[i-fs] = 0.; 258 CHKERRQ(MatGetRow(A,i,&ncols,&rcol,&rval)); 259 for (j=0;j<ncols;j++) { 260 if ((PetscRealPart(-rval[j]) > Amax_neg[i-fs]) && i != rcol[j]) Amax_neg[i-fs] = PetscAbsScalar(rval[j]); 261 if ((PetscRealPart(rval[j]) > Amax_pos[i-fs]) && i != rcol[j]) Amax_pos[i-fs] = PetscAbsScalar(rval[j]); 262 } 263 if (ncols > cmax) cmax = ncols; 264 CHKERRQ(MatRestoreRow(A,i,&ncols,&rcol,&rval)); 265 } 266 CHKERRQ(PetscMalloc2(cmax,&pcols,cmax,&pvals)); 267 CHKERRQ(VecDestroy(&C)); 268 269 /* count the on and off processor sparsity patterns for the prolongator */ 270 for (i=0;i<fn;i++) { 271 /* on */ 272 lsparse[i] = 0; 273 gsparse[i] = 0; 274 if (lcid[i] >= 0) { 275 lsparse[i] = 1; 276 gsparse[i] = 0; 277 } else { 278 CHKERRQ(MatGetRow(lA,i,&ncols,&rcol,&rval)); 279 for (j = 0;j < ncols;j++) { 280 col = rcol[j]; 281 if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) { 282 lsparse[i] += 1; 283 } 284 } 285 CHKERRQ(MatRestoreRow(lA,i,&ncols,&rcol,&rval)); 286 /* off */ 287 if (gA) { 288 CHKERRQ(MatGetRow(gA,i,&ncols,&rcol,&rval)); 289 for (j = 0; j < ncols; j++) { 290 col = rcol[j]; 291 if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) { 292 gsparse[i] += 1; 293 } 294 } 295 CHKERRQ(MatRestoreRow(gA,i,&ncols,&rcol,&rval)); 296 } 297 } 298 } 299 300 /* preallocate and create the prolongator */ 301 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),P)); 302 CHKERRQ(MatGetType(G,&mtype)); 303 CHKERRQ(MatSetType(*P,mtype)); 304 CHKERRQ(MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE)); 305 CHKERRQ(MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse)); 306 CHKERRQ(MatSeqAIJSetPreallocation(*P,0,lsparse)); 307 308 /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */ 309 for (i = 0;i < fn;i++) { 310 /* determine on or off */ 311 row_f = i + fs; 312 row_c = lcid[i]; 313 if (row_c >= 0) { 314 pij = 1.; 315 CHKERRQ(MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES)); 316 } else { 317 g_pos = 0.; 318 g_neg = 0.; 319 a_pos = 0.; 320 a_neg = 0.; 321 diag = 0.; 322 323 /* local connections */ 324 CHKERRQ(MatGetRow(lA,i,&ncols,&rcol,&rval)); 325 for (j = 0; j < ncols; j++) { 326 col = rcol[j]; 327 if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) { 328 if (PetscRealPart(rval[j]) > 0.) { 329 g_pos += rval[j]; 330 } else { 331 g_neg += rval[j]; 332 } 333 } 334 if (col != i) { 335 if (PetscRealPart(rval[j]) > 0.) { 336 a_pos += rval[j]; 337 } else { 338 a_neg += rval[j]; 339 } 340 } else { 341 diag = rval[j]; 342 } 343 } 344 CHKERRQ(MatRestoreRow(lA,i,&ncols,&rcol,&rval)); 345 346 /* ghosted connections */ 347 if (gA) { 348 CHKERRQ(MatGetRow(gA,i,&ncols,&rcol,&rval)); 349 for (j = 0; j < ncols; j++) { 350 col = rcol[j]; 351 if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) { 352 if (PetscRealPart(rval[j]) > 0.) { 353 g_pos += rval[j]; 354 } else { 355 g_neg += rval[j]; 356 } 357 } 358 if (PetscRealPart(rval[j]) > 0.) { 359 a_pos += rval[j]; 360 } else { 361 a_neg += rval[j]; 362 } 363 } 364 CHKERRQ(MatRestoreRow(gA,i,&ncols,&rcol,&rval)); 365 } 366 367 if (g_neg == 0.) { 368 alpha = 0.; 369 } else { 370 alpha = -a_neg/g_neg; 371 } 372 373 if (g_pos == 0.) { 374 diag += a_pos; 375 beta = 0.; 376 } else { 377 beta = -a_pos/g_pos; 378 } 379 if (diag == 0.) { 380 invdiag = 0.; 381 } else invdiag = 1. / diag; 382 /* on */ 383 CHKERRQ(MatGetRow(lA,i,&ncols,&rcol,&rval)); 384 idx = 0; 385 for (j = 0;j < ncols;j++) { 386 col = rcol[j]; 387 if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) { 388 row_f = i + fs; 389 row_c = lcid[col]; 390 /* set the values for on-processor ones */ 391 if (PetscRealPart(rval[j]) < 0.) { 392 pij = rval[j]*alpha*invdiag; 393 } else { 394 pij = rval[j]*beta*invdiag; 395 } 396 if (PetscAbsScalar(pij) != 0.) { 397 pvals[idx] = pij; 398 pcols[idx] = row_c; 399 idx++; 400 } 401 } 402 } 403 CHKERRQ(MatRestoreRow(lA,i,&ncols,&rcol,&rval)); 404 /* off */ 405 if (gA) { 406 CHKERRQ(MatGetRow(gA,i,&ncols,&rcol,&rval)); 407 for (j = 0; j < ncols; j++) { 408 col = rcol[j]; 409 if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) { 410 row_f = i + fs; 411 row_c = gcid[col]; 412 /* set the values for on-processor ones */ 413 if (PetscRealPart(rval[j]) < 0.) { 414 pij = rval[j]*alpha*invdiag; 415 } else { 416 pij = rval[j]*beta*invdiag; 417 } 418 if (PetscAbsScalar(pij) != 0.) { 419 pvals[idx] = pij; 420 pcols[idx] = row_c; 421 idx++; 422 } 423 } 424 } 425 CHKERRQ(MatRestoreRow(gA,i,&ncols,&rcol,&rval)); 426 } 427 CHKERRQ(MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES)); 428 } 429 } 430 431 CHKERRQ(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY)); 432 CHKERRQ(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY)); 433 434 CHKERRQ(PetscFree5(lsparse,gsparse,lcid,Amax_pos,Amax_neg)); 435 436 CHKERRQ(PetscFree2(pcols,pvals)); 437 if (gA) { 438 CHKERRQ(PetscSFDestroy(&sf)); 439 CHKERRQ(PetscFree(gcid)); 440 } 441 PetscFunctionReturn(0); 442 } 443 444 PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc,Mat *P) 445 { 446 PetscInt j,i,ps,pf,pn,pcs,pcf,pcn,idx,cmax; 447 const PetscScalar *pval; 448 const PetscInt *pcol; 449 PetscScalar *pnval; 450 PetscInt *pncol; 451 PetscInt ncols; 452 Mat Pnew; 453 PetscInt *lsparse,*gsparse; 454 PetscReal pmax_pos,pmax_neg,ptot_pos,ptot_neg,pthresh_pos,pthresh_neg; 455 PC_MG *mg = (PC_MG*)pc->data; 456 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 457 PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx; 458 MatType mtype; 459 460 PetscFunctionBegin; 461 /* trim and rescale with reallocation */ 462 CHKERRQ(MatGetOwnershipRange(*P,&ps,&pf)); 463 CHKERRQ(MatGetOwnershipRangeColumn(*P,&pcs,&pcf)); 464 pn = pf-ps; 465 pcn = pcf-pcs; 466 CHKERRQ(PetscMalloc2(pn,&lsparse,pn,&gsparse)); 467 /* allocate */ 468 cmax = 0; 469 for (i=ps;i<pf;i++) { 470 lsparse[i-ps] = 0; 471 gsparse[i-ps] = 0; 472 CHKERRQ(MatGetRow(*P,i,&ncols,&pcol,&pval)); 473 if (ncols > cmax) { 474 cmax = ncols; 475 } 476 pmax_pos = 0.; 477 pmax_neg = 0.; 478 for (j=0;j<ncols;j++) { 479 if (PetscRealPart(pval[j]) > pmax_pos) { 480 pmax_pos = PetscRealPart(pval[j]); 481 } else if (PetscRealPart(pval[j]) < pmax_neg) { 482 pmax_neg = PetscRealPart(pval[j]); 483 } 484 } 485 for (j=0;j<ncols;j++) { 486 if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) { 487 if (pcol[j] >= pcs && pcol[j] < pcf) { 488 lsparse[i-ps]++; 489 } else { 490 gsparse[i-ps]++; 491 } 492 } 493 } 494 CHKERRQ(MatRestoreRow(*P,i,&ncols,&pcol,&pval)); 495 } 496 497 CHKERRQ(PetscMalloc2(cmax,&pnval,cmax,&pncol)); 498 499 CHKERRQ(MatGetType(*P,&mtype)); 500 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)*P),&Pnew)); 501 CHKERRQ(MatSetType(Pnew, mtype)); 502 CHKERRQ(MatSetSizes(Pnew,pn,pcn,PETSC_DETERMINE,PETSC_DETERMINE)); 503 CHKERRQ(MatSeqAIJSetPreallocation(Pnew,0,lsparse)); 504 CHKERRQ(MatMPIAIJSetPreallocation(Pnew,0,lsparse,0,gsparse)); 505 506 for (i=ps;i<pf;i++) { 507 CHKERRQ(MatGetRow(*P,i,&ncols,&pcol,&pval)); 508 pmax_pos = 0.; 509 pmax_neg = 0.; 510 for (j=0;j<ncols;j++) { 511 if (PetscRealPart(pval[j]) > pmax_pos) { 512 pmax_pos = PetscRealPart(pval[j]); 513 } else if (PetscRealPart(pval[j]) < pmax_neg) { 514 pmax_neg = PetscRealPart(pval[j]); 515 } 516 } 517 pthresh_pos = 0.; 518 pthresh_neg = 0.; 519 ptot_pos = 0.; 520 ptot_neg = 0.; 521 for (j=0;j<ncols;j++) { 522 if (PetscRealPart(pval[j]) >= cls->interp_threshold*pmax_pos) { 523 pthresh_pos += PetscRealPart(pval[j]); 524 } else if (PetscRealPart(pval[j]) <= cls->interp_threshold*pmax_neg) { 525 pthresh_neg += PetscRealPart(pval[j]); 526 } 527 if (PetscRealPart(pval[j]) > 0.) { 528 ptot_pos += PetscRealPart(pval[j]); 529 } else { 530 ptot_neg += PetscRealPart(pval[j]); 531 } 532 } 533 if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos; 534 if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg; 535 idx=0; 536 for (j=0;j<ncols;j++) { 537 if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold) { 538 pnval[idx] = ptot_pos*pval[j]; 539 pncol[idx] = pcol[j]; 540 idx++; 541 } else if (PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) { 542 pnval[idx] = ptot_neg*pval[j]; 543 pncol[idx] = pcol[j]; 544 idx++; 545 } 546 } 547 CHKERRQ(MatRestoreRow(*P,i,&ncols,&pcol,&pval)); 548 CHKERRQ(MatSetValues(Pnew,1,&i,idx,pncol,pnval,INSERT_VALUES)); 549 } 550 551 CHKERRQ(MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY)); 552 CHKERRQ(MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY)); 553 CHKERRQ(MatDestroy(P)); 554 555 *P = Pnew; 556 CHKERRQ(PetscFree2(lsparse,gsparse)); 557 CHKERRQ(PetscFree2(pnval,pncol)); 558 PetscFunctionReturn(0); 559 } 560 561 PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P) 562 { 563 Mat lA,*lAs; 564 MatType mtype; 565 Vec cv; 566 PetscInt *gcid,*lcid,*lsparse,*gsparse,*picol; 567 PetscInt fs,fe,cs,ce,nl,i,j,k,li,lni,ci,ncols,maxcols,fn,cn,cid; 568 PetscMPIInt size; 569 const PetscInt *lidx,*icol,*gidx; 570 PetscBool iscoarse; 571 PetscScalar vi,pentry,pjentry; 572 PetscScalar *pcontrib,*pvcol; 573 const PetscScalar *vcol; 574 PetscReal diag,jdiag,jwttotal; 575 PetscInt pncols; 576 PetscSF sf; 577 PetscLayout clayout; 578 IS lis; 579 580 PetscFunctionBegin; 581 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A),&size)); 582 CHKERRQ(MatGetOwnershipRange(A,&fs,&fe)); 583 fn = fe-fs; 584 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,fe-fs,fs,1,&lis)); 585 if (size > 1) { 586 CHKERRQ(MatGetLayouts(A,NULL,&clayout)); 587 /* increase the overlap by two to get neighbors of neighbors */ 588 CHKERRQ(MatIncreaseOverlap(A,1,&lis,2)); 589 CHKERRQ(ISSort(lis)); 590 /* get the local part of A */ 591 CHKERRQ(MatCreateSubMatrices(A,1,&lis,&lis,MAT_INITIAL_MATRIX,&lAs)); 592 lA = lAs[0]; 593 /* build an SF out of it */ 594 CHKERRQ(ISGetLocalSize(lis,&nl)); 595 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject)A),&sf)); 596 CHKERRQ(ISGetIndices(lis,&lidx)); 597 CHKERRQ(PetscSFSetGraphLayout(sf,clayout,nl,NULL,PETSC_COPY_VALUES,lidx)); 598 CHKERRQ(ISRestoreIndices(lis,&lidx)); 599 } else { 600 lA = A; 601 nl = fn; 602 } 603 /* create a communication structure for the overlapped portion and transmit coarse indices */ 604 CHKERRQ(PetscMalloc3(fn,&lsparse,fn,&gsparse,nl,&pcontrib)); 605 /* create coarse vector */ 606 cn = 0; 607 for (i=0;i<fn;i++) { 608 CHKERRQ(PetscCDEmptyAt(agg_lists,i,&iscoarse)); 609 if (!iscoarse) { 610 cn++; 611 } 612 } 613 CHKERRQ(PetscMalloc1(fn,&gcid)); 614 CHKERRQ(VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&cv)); 615 CHKERRQ(VecGetOwnershipRange(cv,&cs,&ce)); 616 cn = 0; 617 for (i=0;i<fn;i++) { 618 CHKERRQ(PetscCDEmptyAt(agg_lists,i,&iscoarse)); 619 if (!iscoarse) { 620 gcid[i] = cs+cn; 621 cn++; 622 } else { 623 gcid[i] = -1; 624 } 625 } 626 if (size > 1) { 627 CHKERRQ(PetscMalloc1(nl,&lcid)); 628 CHKERRQ(PetscSFBcastBegin(sf,MPIU_INT,gcid,lcid,MPI_REPLACE)); 629 CHKERRQ(PetscSFBcastEnd(sf,MPIU_INT,gcid,lcid,MPI_REPLACE)); 630 } else { 631 lcid = gcid; 632 } 633 /* count to preallocate the prolongator */ 634 CHKERRQ(ISGetIndices(lis,&gidx)); 635 maxcols = 0; 636 /* count the number of unique contributing coarse cells for each fine */ 637 for (i=0;i<nl;i++) { 638 pcontrib[i] = 0.; 639 CHKERRQ(MatGetRow(lA,i,&ncols,&icol,NULL)); 640 if (gidx[i] >= fs && gidx[i] < fe) { 641 li = gidx[i] - fs; 642 lsparse[li] = 0; 643 gsparse[li] = 0; 644 cid = lcid[i]; 645 if (cid >= 0) { 646 lsparse[li] = 1; 647 } else { 648 for (j=0;j<ncols;j++) { 649 if (lcid[icol[j]] >= 0) { 650 pcontrib[icol[j]] = 1.; 651 } else { 652 ci = icol[j]; 653 CHKERRQ(MatRestoreRow(lA,i,&ncols,&icol,NULL)); 654 CHKERRQ(MatGetRow(lA,ci,&ncols,&icol,NULL)); 655 for (k=0;k<ncols;k++) { 656 if (lcid[icol[k]] >= 0) { 657 pcontrib[icol[k]] = 1.; 658 } 659 } 660 CHKERRQ(MatRestoreRow(lA,ci,&ncols,&icol,NULL)); 661 CHKERRQ(MatGetRow(lA,i,&ncols,&icol,NULL)); 662 } 663 } 664 for (j=0;j<ncols;j++) { 665 if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) { 666 lni = lcid[icol[j]]; 667 if (lni >= cs && lni < ce) { 668 lsparse[li]++; 669 } else { 670 gsparse[li]++; 671 } 672 pcontrib[icol[j]] = 0.; 673 } else { 674 ci = icol[j]; 675 CHKERRQ(MatRestoreRow(lA,i,&ncols,&icol,NULL)); 676 CHKERRQ(MatGetRow(lA,ci,&ncols,&icol,NULL)); 677 for (k=0;k<ncols;k++) { 678 if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) { 679 lni = lcid[icol[k]]; 680 if (lni >= cs && lni < ce) { 681 lsparse[li]++; 682 } else { 683 gsparse[li]++; 684 } 685 pcontrib[icol[k]] = 0.; 686 } 687 } 688 CHKERRQ(MatRestoreRow(lA,ci,&ncols,&icol,NULL)); 689 CHKERRQ(MatGetRow(lA,i,&ncols,&icol,NULL)); 690 } 691 } 692 } 693 if (lsparse[li] + gsparse[li] > maxcols) maxcols = lsparse[li]+gsparse[li]; 694 } 695 CHKERRQ(MatRestoreRow(lA,i,&ncols,&icol,&vcol)); 696 } 697 CHKERRQ(PetscMalloc2(maxcols,&picol,maxcols,&pvcol)); 698 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),P)); 699 CHKERRQ(MatGetType(A,&mtype)); 700 CHKERRQ(MatSetType(*P,mtype)); 701 CHKERRQ(MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE)); 702 CHKERRQ(MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse)); 703 CHKERRQ(MatSeqAIJSetPreallocation(*P,0,lsparse)); 704 for (i=0;i<nl;i++) { 705 diag = 0.; 706 if (gidx[i] >= fs && gidx[i] < fe) { 707 pncols=0; 708 cid = lcid[i]; 709 if (cid >= 0) { 710 pncols = 1; 711 picol[0] = cid; 712 pvcol[0] = 1.; 713 } else { 714 CHKERRQ(MatGetRow(lA,i,&ncols,&icol,&vcol)); 715 for (j=0;j<ncols;j++) { 716 pentry = vcol[j]; 717 if (lcid[icol[j]] >= 0) { 718 /* coarse neighbor */ 719 pcontrib[icol[j]] += pentry; 720 } else if (icol[j] != i) { 721 /* the neighbor is a strongly connected fine node */ 722 ci = icol[j]; 723 vi = vcol[j]; 724 CHKERRQ(MatRestoreRow(lA,i,&ncols,&icol,&vcol)); 725 CHKERRQ(MatGetRow(lA,ci,&ncols,&icol,&vcol)); 726 jwttotal=0.; 727 jdiag = 0.; 728 for (k=0;k<ncols;k++) { 729 if (ci == icol[k]) { 730 jdiag = PetscRealPart(vcol[k]); 731 } 732 } 733 for (k=0;k<ncols;k++) { 734 if (lcid[icol[k]] >= 0 && jdiag*PetscRealPart(vcol[k]) < 0.) { 735 pjentry = vcol[k]; 736 jwttotal += PetscRealPart(pjentry); 737 } 738 } 739 if (jwttotal != 0.) { 740 jwttotal = PetscRealPart(vi)/jwttotal; 741 for (k=0;k<ncols;k++) { 742 if (lcid[icol[k]] >= 0 && jdiag*PetscRealPart(vcol[k]) < 0.) { 743 pjentry = vcol[k]*jwttotal; 744 pcontrib[icol[k]] += pjentry; 745 } 746 } 747 } else { 748 diag += PetscRealPart(vi); 749 } 750 CHKERRQ(MatRestoreRow(lA,ci,&ncols,&icol,&vcol)); 751 CHKERRQ(MatGetRow(lA,i,&ncols,&icol,&vcol)); 752 } else { 753 diag += PetscRealPart(vcol[j]); 754 } 755 } 756 if (diag != 0.) { 757 diag = 1./diag; 758 for (j=0;j<ncols;j++) { 759 if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) { 760 /* the neighbor is a coarse node */ 761 if (PetscAbsScalar(pcontrib[icol[j]]) > 0.0) { 762 lni = lcid[icol[j]]; 763 pvcol[pncols] = -pcontrib[icol[j]]*diag; 764 picol[pncols] = lni; 765 pncols++; 766 } 767 pcontrib[icol[j]] = 0.; 768 } else { 769 /* the neighbor is a strongly connected fine node */ 770 ci = icol[j]; 771 CHKERRQ(MatRestoreRow(lA,i,&ncols,&icol,&vcol)); 772 CHKERRQ(MatGetRow(lA,ci,&ncols,&icol,&vcol)); 773 for (k=0;k<ncols;k++) { 774 if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) { 775 if (PetscAbsScalar(pcontrib[icol[k]]) > 0.0) { 776 lni = lcid[icol[k]]; 777 pvcol[pncols] = -pcontrib[icol[k]]*diag; 778 picol[pncols] = lni; 779 pncols++; 780 } 781 pcontrib[icol[k]] = 0.; 782 } 783 } 784 CHKERRQ(MatRestoreRow(lA,ci,&ncols,&icol,&vcol)); 785 CHKERRQ(MatGetRow(lA,i,&ncols,&icol,&vcol)); 786 } 787 pcontrib[icol[j]] = 0.; 788 } 789 CHKERRQ(MatRestoreRow(lA,i,&ncols,&icol,&vcol)); 790 } 791 } 792 ci = gidx[i]; 793 if (pncols > 0) { 794 CHKERRQ(MatSetValues(*P,1,&ci,pncols,picol,pvcol,INSERT_VALUES)); 795 } 796 } 797 } 798 CHKERRQ(ISRestoreIndices(lis,&gidx)); 799 CHKERRQ(PetscFree2(picol,pvcol)); 800 CHKERRQ(PetscFree3(lsparse,gsparse,pcontrib)); 801 CHKERRQ(ISDestroy(&lis)); 802 CHKERRQ(PetscFree(gcid)); 803 if (size > 1) { 804 CHKERRQ(PetscFree(lcid)); 805 CHKERRQ(MatDestroyMatrices(1,&lAs)); 806 CHKERRQ(PetscSFDestroy(&sf)); 807 } 808 CHKERRQ(VecDestroy(&cv)); 809 CHKERRQ(MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY)); 810 CHKERRQ(MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY)); 811 PetscFunctionReturn(0); 812 } 813 814 PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc,Mat A,Mat *P) 815 { 816 817 PetscInt f,s,n,cf,cs,i,idx; 818 PetscInt *coarserows; 819 PetscInt ncols; 820 const PetscInt *pcols; 821 const PetscScalar *pvals; 822 Mat Pnew; 823 Vec diag; 824 PC_MG *mg = (PC_MG*)pc->data; 825 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 826 PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx; 827 828 PetscFunctionBegin; 829 if (cls->nsmooths == 0) { 830 CHKERRQ(PCGAMGTruncateProlongator_Private(pc,P)); 831 PetscFunctionReturn(0); 832 } 833 CHKERRQ(MatGetOwnershipRange(*P,&s,&f)); 834 n = f-s; 835 CHKERRQ(MatGetOwnershipRangeColumn(*P,&cs,&cf)); 836 CHKERRQ(PetscMalloc1(n,&coarserows)); 837 /* identify the rows corresponding to coarse unknowns */ 838 idx = 0; 839 for (i=s;i<f;i++) { 840 CHKERRQ(MatGetRow(*P,i,&ncols,&pcols,&pvals)); 841 /* assume, for now, that it's a coarse unknown if it has a single unit entry */ 842 if (ncols == 1) { 843 if (pvals[0] == 1.) { 844 coarserows[idx] = i; 845 idx++; 846 } 847 } 848 CHKERRQ(MatRestoreRow(*P,i,&ncols,&pcols,&pvals)); 849 } 850 CHKERRQ(MatCreateVecs(A,&diag,NULL)); 851 CHKERRQ(MatGetDiagonal(A,diag)); 852 CHKERRQ(VecReciprocal(diag)); 853 for (i=0;i<cls->nsmooths;i++) { 854 CHKERRQ(MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew)); 855 CHKERRQ(MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL)); 856 CHKERRQ(MatDiagonalScale(Pnew,diag,NULL)); 857 CHKERRQ(MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN)); 858 CHKERRQ(MatDestroy(P)); 859 *P = Pnew; 860 Pnew = NULL; 861 } 862 CHKERRQ(VecDestroy(&diag)); 863 CHKERRQ(PetscFree(coarserows)); 864 CHKERRQ(PCGAMGTruncateProlongator_Private(pc,P)); 865 PetscFunctionReturn(0); 866 } 867 868 PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P) 869 { 870 PetscErrorCode (*f)(PC,Mat,Mat,PetscCoarsenData*,Mat*); 871 PC_MG *mg = (PC_MG*)pc->data; 872 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 873 PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx; 874 875 PetscFunctionBegin; 876 CHKERRQ(PetscFunctionListFind(PCGAMGClassicalProlongatorList,cls->prolongtype,&f)); 877 PetscCheck(f,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot find PCGAMG Classical prolongator type"); 878 CHKERRQ((*f)(pc,A,G,agg_lists,P)); 879 PetscFunctionReturn(0); 880 } 881 882 PetscErrorCode PCGAMGDestroy_Classical(PC pc) 883 { 884 PC_MG *mg = (PC_MG*)pc->data; 885 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 886 887 PetscFunctionBegin; 888 CHKERRQ(PetscFree(pc_gamg->subctx)); 889 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",NULL)); 890 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",NULL)); 891 PetscFunctionReturn(0); 892 } 893 894 PetscErrorCode PCGAMGSetFromOptions_Classical(PetscOptionItems *PetscOptionsObject,PC pc) 895 { 896 PC_MG *mg = (PC_MG*)pc->data; 897 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 898 PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx; 899 char tname[256]; 900 PetscBool flg; 901 902 PetscFunctionBegin; 903 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"GAMG-Classical options")); 904 CHKERRQ(PetscOptionsFList("-pc_gamg_classical_type","Type of Classical AMG prolongation","PCGAMGClassicalSetType",PCGAMGClassicalProlongatorList,cls->prolongtype, tname, sizeof(tname), &flg)); 905 if (flg) { 906 CHKERRQ(PCGAMGClassicalSetType(pc,tname)); 907 } 908 CHKERRQ(PetscOptionsReal("-pc_gamg_classical_interp_threshold","Threshold for classical interpolator entries","",cls->interp_threshold,&cls->interp_threshold,NULL)); 909 CHKERRQ(PetscOptionsInt("-pc_gamg_classical_nsmooths","Threshold for classical interpolator entries","",cls->nsmooths,&cls->nsmooths,NULL)); 910 CHKERRQ(PetscOptionsTail()); 911 PetscFunctionReturn(0); 912 } 913 914 PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A) 915 { 916 PC_MG *mg = (PC_MG*)pc->data; 917 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 918 919 PetscFunctionBegin; 920 /* no data for classical AMG */ 921 pc_gamg->data = NULL; 922 pc_gamg->data_cell_cols = 0; 923 pc_gamg->data_cell_rows = 0; 924 pc_gamg->data_sz = 0; 925 PetscFunctionReturn(0); 926 } 927 928 PetscErrorCode PCGAMGClassicalFinalizePackage(void) 929 { 930 PetscFunctionBegin; 931 PCGAMGClassicalPackageInitialized = PETSC_FALSE; 932 CHKERRQ(PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList)); 933 PetscFunctionReturn(0); 934 } 935 936 PetscErrorCode PCGAMGClassicalInitializePackage(void) 937 { 938 PetscFunctionBegin; 939 if (PCGAMGClassicalPackageInitialized) PetscFunctionReturn(0); 940 CHKERRQ(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALDIRECT,PCGAMGProlongator_Classical_Direct)); 941 CHKERRQ(PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALSTANDARD,PCGAMGProlongator_Classical_Standard)); 942 CHKERRQ(PetscRegisterFinalize(PCGAMGClassicalFinalizePackage)); 943 PetscFunctionReturn(0); 944 } 945 946 /* -------------------------------------------------------------------------- */ 947 /* 948 PCCreateGAMG_Classical 949 950 */ 951 PetscErrorCode PCCreateGAMG_Classical(PC pc) 952 { 953 PC_MG *mg = (PC_MG*)pc->data; 954 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 955 PC_GAMG_Classical *pc_gamg_classical; 956 957 PetscFunctionBegin; 958 CHKERRQ(PCGAMGClassicalInitializePackage()); 959 if (pc_gamg->subctx) { 960 /* call base class */ 961 CHKERRQ(PCDestroy_GAMG(pc)); 962 } 963 964 /* create sub context for SA */ 965 CHKERRQ(PetscNewLog(pc,&pc_gamg_classical)); 966 pc_gamg->subctx = pc_gamg_classical; 967 pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical; 968 /* reset does not do anything; setup not virtual */ 969 970 /* set internal function pointers */ 971 pc_gamg->ops->destroy = PCGAMGDestroy_Classical; 972 pc_gamg->ops->graph = PCGAMGGraph_Classical; 973 pc_gamg->ops->coarsen = PCGAMGCoarsen_Classical; 974 pc_gamg->ops->prolongator = PCGAMGProlongator_Classical; 975 pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi; 976 pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical; 977 978 pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical; 979 pc_gamg_classical->interp_threshold = 0.2; 980 pc_gamg_classical->nsmooths = 0; 981 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",PCGAMGClassicalSetType_GAMG)); 982 CHKERRQ(PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",PCGAMGClassicalGetType_GAMG)); 983 CHKERRQ(PCGAMGClassicalSetType(pc,PCGAMGCLASSICALSTANDARD)); 984 PetscFunctionReturn(0); 985 } 986