Actual source code: classical.c

  1: #include <../src/ksp/pc/impls/gamg/gamg.h>
  2: #include <petscsf.h>

  4: PetscFunctionList PCGAMGClassicalProlongatorList    = NULL;
  5: PetscBool         PCGAMGClassicalPackageInitialized = PETSC_FALSE;

  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;

 13: /*@C
 14:    PCGAMGClassicalSetType - Sets the type of classical interpolation to use

 16:    Collective on PC

 18:    Input Parameters:
 19: .  pc - the preconditioner context

 21:    Options Database Key:
 22: .  -pc_gamg_classical_type <direct,standard> - set type of Classical AMG prolongation

 24:    Level: intermediate

 26: .seealso: ()
 27: @*/
 28: PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type)
 29: {
 31:   PetscTryMethod(pc,"PCGAMGClassicalSetType_C",(PC,PCGAMGClassicalType),(pc,type));
 32:   return 0;
 33: }

 35: /*@C
 36:    PCGAMGClassicalGetType - Gets the type of classical interpolation to use

 38:    Collective on PC

 40:    Input Parameter:
 41: .  pc - the preconditioner context

 43:    Output Parameter:
 44: .  type - the type used

 46:    Level: intermediate

 48: .seealso: ()
 49: @*/
 50: PetscErrorCode PCGAMGClassicalGetType(PC pc, PCGAMGClassicalType *type)
 51: {
 53:   PetscUseMethod(pc,"PCGAMGClassicalGetType_C",(PC,PCGAMGClassicalType*),(pc,type));
 54:   return 0;
 55: }

 57: static PetscErrorCode PCGAMGClassicalSetType_GAMG(PC pc, PCGAMGClassicalType type)
 58: {
 59:   PC_MG             *mg          = (PC_MG*)pc->data;
 60:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
 61:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

 63:   PetscStrcpy(cls->prolongtype,type);
 64:   return 0;
 65: }

 67: static PetscErrorCode PCGAMGClassicalGetType_GAMG(PC pc, PCGAMGClassicalType *type)
 68: {
 69:   PC_MG             *mg          = (PC_MG*)pc->data;
 70:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
 71:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

 73:   *type = cls->prolongtype;
 74:   return 0;
 75: }

 77: PetscErrorCode PCGAMGGraph_Classical(PC pc,Mat A,Mat *G)
 78: {
 79:   PetscInt          s,f,n,idx,lidx,gidx;
 80:   PetscInt          r,c,ncols;
 81:   const PetscInt    *rcol;
 82:   const PetscScalar *rval;
 83:   PetscInt          *gcol;
 84:   PetscScalar       *gval;
 85:   PetscReal         rmax;
 86:   PetscInt          cmax = 0;
 87:   PC_MG             *mg = (PC_MG *)pc->data;
 88:   PC_GAMG           *gamg = (PC_GAMG *)mg->innerctx;
 89:   PetscInt          *gsparse,*lsparse;
 90:   PetscScalar       *Amax;
 91:   MatType           mtype;

 93:   MatGetOwnershipRange(A,&s,&f);
 94:   n=f-s;
 95:   PetscMalloc3(n,&lsparse,n,&gsparse,n,&Amax);

 97:   for (r = 0;r < n;r++) {
 98:     lsparse[r] = 0;
 99:     gsparse[r] = 0;
100:   }

102:   for (r = s;r < f;r++) {
103:     /* determine the maximum off-diagonal in each row */
104:     rmax = 0.;
105:     MatGetRow(A,r,&ncols,&rcol,&rval);
106:     for (c = 0; c < ncols; c++) {
107:       if (PetscRealPart(-rval[c]) > rmax && rcol[c] != r) {
108:         rmax = PetscRealPart(-rval[c]);
109:       }
110:     }
111:     Amax[r-s] = rmax;
112:     if (ncols > cmax) cmax = ncols;
113:     lidx = 0;
114:     gidx = 0;
115:     /* create the local and global sparsity patterns */
116:     for (c = 0; c < ncols; c++) {
117:       if (PetscRealPart(-rval[c]) > gamg->threshold[0]*PetscRealPart(Amax[r-s]) || rcol[c] == r) {
118:         if (rcol[c] < f && rcol[c] >= s) {
119:           lidx++;
120:         } else {
121:           gidx++;
122:         }
123:       }
124:     }
125:     MatRestoreRow(A,r,&ncols,&rcol,&rval);
126:     lsparse[r-s] = lidx;
127:     gsparse[r-s] = gidx;
128:   }
129:   PetscMalloc2(cmax,&gval,cmax,&gcol);

131:   MatCreate(PetscObjectComm((PetscObject)A),G);
132:   MatGetType(A,&mtype);
133:   MatSetType(*G,mtype);
134:   MatSetSizes(*G,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
135:   MatMPIAIJSetPreallocation(*G,0,lsparse,0,gsparse);
136:   MatSeqAIJSetPreallocation(*G,0,lsparse);
137:   for (r = s;r < f;r++) {
138:     MatGetRow(A,r,&ncols,&rcol,&rval);
139:     idx = 0;
140:     for (c = 0; c < ncols; c++) {
141:       /* classical strength of connection */
142:       if (PetscRealPart(-rval[c]) > gamg->threshold[0]*PetscRealPart(Amax[r-s]) || rcol[c] == r) {
143:         gcol[idx] = rcol[c];
144:         gval[idx] = rval[c];
145:         idx++;
146:       }
147:     }
148:     MatSetValues(*G,1,&r,idx,gcol,gval,INSERT_VALUES);
149:     MatRestoreRow(A,r,&ncols,&rcol,&rval);
150:   }
151:   MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY);
152:   MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY);

154:   PetscFree2(gval,gcol);
155:   PetscFree3(lsparse,gsparse,Amax);
156:   return 0;
157: }

159: PetscErrorCode PCGAMGCoarsen_Classical(PC pc,Mat *G,PetscCoarsenData **agg_lists)
160: {
161:   MatCoarsen       crs;
162:   MPI_Comm         fcomm = ((PetscObject)pc)->comm;


166:   MatCoarsenCreate(fcomm,&crs);
167:   MatCoarsenSetFromOptions(crs);
168:   MatCoarsenSetAdjacency(crs,*G);
169:   MatCoarsenSetStrictAggs(crs,PETSC_TRUE);
170:   MatCoarsenApply(crs);
171:   MatCoarsenGetData(crs,agg_lists);
172:   MatCoarsenDestroy(&crs);
173:   return 0;
174: }

176: PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
177: {
178:   PC_MG             *mg          = (PC_MG*)pc->data;
179:   PC_GAMG           *gamg        = (PC_GAMG*)mg->innerctx;
180:   PetscBool         iscoarse,isMPIAIJ,isSEQAIJ;
181:   PetscInt          fn,cn,fs,fe,cs,ce,i,j,ncols,col,row_f,row_c,cmax=0,idx,noff;
182:   PetscInt          *lcid,*gcid,*lsparse,*gsparse,*colmap,*pcols;
183:   const PetscInt    *rcol;
184:   PetscReal         *Amax_pos,*Amax_neg;
185:   PetscScalar       g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta,pij;
186:   PetscScalar       *pvals;
187:   const PetscScalar *rval;
188:   Mat               lA,gA=NULL;
189:   MatType           mtype;
190:   Vec               C,lvec;
191:   PetscLayout       clayout;
192:   PetscSF           sf;
193:   Mat_MPIAIJ        *mpiaij;

195:   MatGetOwnershipRange(A,&fs,&fe);
196:   fn = fe-fs;
197:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
198:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSEQAIJ);
200:   if (isMPIAIJ) {
201:     mpiaij = (Mat_MPIAIJ*)A->data;
202:     lA = mpiaij->A;
203:     gA = mpiaij->B;
204:     lvec = mpiaij->lvec;
205:     VecGetSize(lvec,&noff);
206:     colmap = mpiaij->garray;
207:     MatGetLayouts(A,NULL,&clayout);
208:     PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
209:     PetscSFSetGraphLayout(sf,clayout,noff,NULL,PETSC_COPY_VALUES,colmap);
210:     PetscMalloc1(noff,&gcid);
211:   } else {
212:     lA = A;
213:   }
214:   PetscMalloc5(fn,&lsparse,fn,&gsparse,fn,&lcid,fn,&Amax_pos,fn,&Amax_neg);

216:   /* count the number of coarse unknowns */
217:   cn = 0;
218:   for (i=0;i<fn;i++) {
219:     /* filter out singletons */
220:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
221:     lcid[i] = -1;
222:     if (!iscoarse) {
223:       cn++;
224:     }
225:   }

227:    /* create the coarse vector */
228:   VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&C);
229:   VecGetOwnershipRange(C,&cs,&ce);

231:   cn = 0;
232:   for (i=0;i<fn;i++) {
233:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
234:     if (!iscoarse) {
235:       lcid[i] = cs+cn;
236:       cn++;
237:     } else {
238:       lcid[i] = -1;
239:     }
240:   }

242:   if (gA) {
243:     PetscSFBcastBegin(sf,MPIU_INT,lcid,gcid,MPI_REPLACE);
244:     PetscSFBcastEnd(sf,MPIU_INT,lcid,gcid,MPI_REPLACE);
245:   }

247:   /* determine the largest off-diagonal entries in each row */
248:   for (i=fs;i<fe;i++) {
249:     Amax_pos[i-fs] = 0.;
250:     Amax_neg[i-fs] = 0.;
251:     MatGetRow(A,i,&ncols,&rcol,&rval);
252:     for (j=0;j<ncols;j++) {
253:       if ((PetscRealPart(-rval[j]) > Amax_neg[i-fs]) && i != rcol[j]) Amax_neg[i-fs] = PetscAbsScalar(rval[j]);
254:       if ((PetscRealPart(rval[j])  > Amax_pos[i-fs]) && i != rcol[j]) Amax_pos[i-fs] = PetscAbsScalar(rval[j]);
255:     }
256:     if (ncols > cmax) cmax = ncols;
257:     MatRestoreRow(A,i,&ncols,&rcol,&rval);
258:   }
259:   PetscMalloc2(cmax,&pcols,cmax,&pvals);
260:   VecDestroy(&C);

262:   /* count the on and off processor sparsity patterns for the prolongator */
263:   for (i=0;i<fn;i++) {
264:     /* on */
265:     lsparse[i] = 0;
266:     gsparse[i] = 0;
267:     if (lcid[i] >= 0) {
268:       lsparse[i] = 1;
269:       gsparse[i] = 0;
270:     } else {
271:       MatGetRow(lA,i,&ncols,&rcol,&rval);
272:       for (j = 0;j < ncols;j++) {
273:         col = rcol[j];
274:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
275:           lsparse[i] += 1;
276:         }
277:       }
278:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);
279:       /* off */
280:       if (gA) {
281:         MatGetRow(gA,i,&ncols,&rcol,&rval);
282:         for (j = 0; j < ncols; j++) {
283:           col = rcol[j];
284:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
285:             gsparse[i] += 1;
286:           }
287:         }
288:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
289:       }
290:     }
291:   }

293:   /* preallocate and create the prolongator */
294:   MatCreate(PetscObjectComm((PetscObject)A),P);
295:   MatGetType(G,&mtype);
296:   MatSetType(*P,mtype);
297:   MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);
298:   MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);
299:   MatSeqAIJSetPreallocation(*P,0,lsparse);

301:   /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
302:   for (i = 0;i < fn;i++) {
303:     /* determine on or off */
304:     row_f = i + fs;
305:     row_c = lcid[i];
306:     if (row_c >= 0) {
307:       pij = 1.;
308:       MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);
309:     } else {
310:       g_pos = 0.;
311:       g_neg = 0.;
312:       a_pos = 0.;
313:       a_neg = 0.;
314:       diag  = 0.;

316:       /* local connections */
317:       MatGetRow(lA,i,&ncols,&rcol,&rval);
318:       for (j = 0; j < ncols; j++) {
319:         col = rcol[j];
320:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
321:           if (PetscRealPart(rval[j]) > 0.) {
322:             g_pos += rval[j];
323:           } else {
324:             g_neg += rval[j];
325:           }
326:         }
327:         if (col != i) {
328:           if (PetscRealPart(rval[j]) > 0.) {
329:             a_pos += rval[j];
330:           } else {
331:             a_neg += rval[j];
332:           }
333:         } else {
334:           diag = rval[j];
335:         }
336:       }
337:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);

339:       /* ghosted connections */
340:       if (gA) {
341:         MatGetRow(gA,i,&ncols,&rcol,&rval);
342:         for (j = 0; j < ncols; j++) {
343:           col = rcol[j];
344:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
345:             if (PetscRealPart(rval[j]) > 0.) {
346:               g_pos += rval[j];
347:             } else {
348:               g_neg += rval[j];
349:             }
350:           }
351:           if (PetscRealPart(rval[j]) > 0.) {
352:             a_pos += rval[j];
353:           } else {
354:             a_neg += rval[j];
355:           }
356:         }
357:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
358:       }

360:       if (g_neg == 0.) {
361:         alpha = 0.;
362:       } else {
363:         alpha = -a_neg/g_neg;
364:       }

366:       if (g_pos == 0.) {
367:         diag += a_pos;
368:         beta = 0.;
369:       } else {
370:         beta = -a_pos/g_pos;
371:       }
372:       if (diag == 0.) {
373:         invdiag = 0.;
374:       } else invdiag = 1. / diag;
375:       /* on */
376:       MatGetRow(lA,i,&ncols,&rcol,&rval);
377:       idx = 0;
378:       for (j = 0;j < ncols;j++) {
379:         col = rcol[j];
380:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
381:           row_f = i + fs;
382:           row_c = lcid[col];
383:           /* set the values for on-processor ones */
384:           if (PetscRealPart(rval[j]) < 0.) {
385:             pij = rval[j]*alpha*invdiag;
386:           } else {
387:             pij = rval[j]*beta*invdiag;
388:           }
389:           if (PetscAbsScalar(pij) != 0.) {
390:             pvals[idx] = pij;
391:             pcols[idx] = row_c;
392:             idx++;
393:           }
394:         }
395:       }
396:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);
397:       /* off */
398:       if (gA) {
399:         MatGetRow(gA,i,&ncols,&rcol,&rval);
400:         for (j = 0; j < ncols; j++) {
401:           col = rcol[j];
402:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
403:             row_f = i + fs;
404:             row_c = gcid[col];
405:             /* set the values for on-processor ones */
406:             if (PetscRealPart(rval[j]) < 0.) {
407:               pij = rval[j]*alpha*invdiag;
408:             } else {
409:               pij = rval[j]*beta*invdiag;
410:             }
411:             if (PetscAbsScalar(pij) != 0.) {
412:               pvals[idx] = pij;
413:               pcols[idx] = row_c;
414:               idx++;
415:             }
416:           }
417:         }
418:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
419:       }
420:       MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);
421:     }
422:   }

424:   MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
425:   MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);

427:   PetscFree5(lsparse,gsparse,lcid,Amax_pos,Amax_neg);

429:   PetscFree2(pcols,pvals);
430:   if (gA) {
431:     PetscSFDestroy(&sf);
432:     PetscFree(gcid);
433:   }
434:   return 0;
435: }

437: PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc,Mat *P)
438: {
439:   PetscInt          j,i,ps,pf,pn,pcs,pcf,pcn,idx,cmax;
440:   const PetscScalar *pval;
441:   const PetscInt    *pcol;
442:   PetscScalar       *pnval;
443:   PetscInt          *pncol;
444:   PetscInt          ncols;
445:   Mat               Pnew;
446:   PetscInt          *lsparse,*gsparse;
447:   PetscReal         pmax_pos,pmax_neg,ptot_pos,ptot_neg,pthresh_pos,pthresh_neg;
448:   PC_MG             *mg          = (PC_MG*)pc->data;
449:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
450:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
451:   MatType           mtype;

453:   /* trim and rescale with reallocation */
454:   MatGetOwnershipRange(*P,&ps,&pf);
455:   MatGetOwnershipRangeColumn(*P,&pcs,&pcf);
456:   pn = pf-ps;
457:   pcn = pcf-pcs;
458:   PetscMalloc2(pn,&lsparse,pn,&gsparse);
459:   /* allocate */
460:   cmax = 0;
461:   for (i=ps;i<pf;i++) {
462:     lsparse[i-ps] = 0;
463:     gsparse[i-ps] = 0;
464:     MatGetRow(*P,i,&ncols,&pcol,&pval);
465:     if (ncols > cmax) {
466:       cmax = ncols;
467:     }
468:     pmax_pos = 0.;
469:     pmax_neg = 0.;
470:     for (j=0;j<ncols;j++) {
471:       if (PetscRealPart(pval[j]) > pmax_pos) {
472:         pmax_pos = PetscRealPart(pval[j]);
473:       } else if (PetscRealPart(pval[j]) < pmax_neg) {
474:         pmax_neg = PetscRealPart(pval[j]);
475:       }
476:     }
477:     for (j=0;j<ncols;j++) {
478:       if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
479:         if (pcol[j] >= pcs && pcol[j] < pcf) {
480:           lsparse[i-ps]++;
481:         } else {
482:           gsparse[i-ps]++;
483:         }
484:       }
485:     }
486:     MatRestoreRow(*P,i,&ncols,&pcol,&pval);
487:   }

489:   PetscMalloc2(cmax,&pnval,cmax,&pncol);

491:   MatGetType(*P,&mtype);
492:   MatCreate(PetscObjectComm((PetscObject)*P),&Pnew);
493:   MatSetType(Pnew, mtype);
494:   MatSetSizes(Pnew,pn,pcn,PETSC_DETERMINE,PETSC_DETERMINE);
495:   MatSeqAIJSetPreallocation(Pnew,0,lsparse);
496:   MatMPIAIJSetPreallocation(Pnew,0,lsparse,0,gsparse);

498:   for (i=ps;i<pf;i++) {
499:     MatGetRow(*P,i,&ncols,&pcol,&pval);
500:     pmax_pos = 0.;
501:     pmax_neg = 0.;
502:     for (j=0;j<ncols;j++) {
503:       if (PetscRealPart(pval[j]) > pmax_pos) {
504:         pmax_pos = PetscRealPart(pval[j]);
505:       } else if (PetscRealPart(pval[j]) < pmax_neg) {
506:         pmax_neg = PetscRealPart(pval[j]);
507:       }
508:     }
509:     pthresh_pos = 0.;
510:     pthresh_neg = 0.;
511:     ptot_pos = 0.;
512:     ptot_neg = 0.;
513:     for (j=0;j<ncols;j++) {
514:       if (PetscRealPart(pval[j]) >= cls->interp_threshold*pmax_pos) {
515:         pthresh_pos += PetscRealPart(pval[j]);
516:       } else if (PetscRealPart(pval[j]) <= cls->interp_threshold*pmax_neg) {
517:         pthresh_neg += PetscRealPart(pval[j]);
518:       }
519:       if (PetscRealPart(pval[j]) > 0.) {
520:         ptot_pos += PetscRealPart(pval[j]);
521:       } else {
522:         ptot_neg += PetscRealPart(pval[j]);
523:       }
524:     }
525:     if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
526:     if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
527:     idx=0;
528:     for (j=0;j<ncols;j++) {
529:       if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold) {
530:         pnval[idx] = ptot_pos*pval[j];
531:         pncol[idx] = pcol[j];
532:         idx++;
533:       } else if (PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
534:         pnval[idx] = ptot_neg*pval[j];
535:         pncol[idx] = pcol[j];
536:         idx++;
537:       }
538:     }
539:     MatRestoreRow(*P,i,&ncols,&pcol,&pval);
540:     MatSetValues(Pnew,1,&i,idx,pncol,pnval,INSERT_VALUES);
541:   }

543:   MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY);
544:   MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY);
545:   MatDestroy(P);

547:   *P = Pnew;
548:   PetscFree2(lsparse,gsparse);
549:   PetscFree2(pnval,pncol);
550:   return 0;
551: }

553: PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
554: {
555:   Mat               lA,*lAs;
556:   MatType           mtype;
557:   Vec               cv;
558:   PetscInt          *gcid,*lcid,*lsparse,*gsparse,*picol;
559:   PetscInt          fs,fe,cs,ce,nl,i,j,k,li,lni,ci,ncols,maxcols,fn,cn,cid;
560:   PetscMPIInt       size;
561:   const PetscInt    *lidx,*icol,*gidx;
562:   PetscBool         iscoarse;
563:   PetscScalar       vi,pentry,pjentry;
564:   PetscScalar       *pcontrib,*pvcol;
565:   const PetscScalar *vcol;
566:   PetscReal         diag,jdiag,jwttotal;
567:   PetscInt          pncols;
568:   PetscSF           sf;
569:   PetscLayout       clayout;
570:   IS                lis;

572:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
573:   MatGetOwnershipRange(A,&fs,&fe);
574:   fn = fe-fs;
575:   ISCreateStride(PETSC_COMM_SELF,fe-fs,fs,1,&lis);
576:   if (size > 1) {
577:     MatGetLayouts(A,NULL,&clayout);
578:     /* increase the overlap by two to get neighbors of neighbors */
579:     MatIncreaseOverlap(A,1,&lis,2);
580:     ISSort(lis);
581:     /* get the local part of A */
582:     MatCreateSubMatrices(A,1,&lis,&lis,MAT_INITIAL_MATRIX,&lAs);
583:     lA = lAs[0];
584:     /* build an SF out of it */
585:     ISGetLocalSize(lis,&nl);
586:     PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
587:     ISGetIndices(lis,&lidx);
588:     PetscSFSetGraphLayout(sf,clayout,nl,NULL,PETSC_COPY_VALUES,lidx);
589:     ISRestoreIndices(lis,&lidx);
590:   } else {
591:     lA = A;
592:     nl = fn;
593:   }
594:   /* create a communication structure for the overlapped portion and transmit coarse indices */
595:   PetscMalloc3(fn,&lsparse,fn,&gsparse,nl,&pcontrib);
596:   /* create coarse vector */
597:   cn = 0;
598:   for (i=0;i<fn;i++) {
599:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
600:     if (!iscoarse) {
601:       cn++;
602:     }
603:   }
604:   PetscMalloc1(fn,&gcid);
605:   VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&cv);
606:   VecGetOwnershipRange(cv,&cs,&ce);
607:   cn = 0;
608:   for (i=0;i<fn;i++) {
609:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
610:     if (!iscoarse) {
611:       gcid[i] = cs+cn;
612:       cn++;
613:     } else {
614:       gcid[i] = -1;
615:     }
616:   }
617:   if (size > 1) {
618:     PetscMalloc1(nl,&lcid);
619:     PetscSFBcastBegin(sf,MPIU_INT,gcid,lcid,MPI_REPLACE);
620:     PetscSFBcastEnd(sf,MPIU_INT,gcid,lcid,MPI_REPLACE);
621:   } else {
622:     lcid = gcid;
623:   }
624:   /* count to preallocate the prolongator */
625:   ISGetIndices(lis,&gidx);
626:   maxcols = 0;
627:   /* count the number of unique contributing coarse cells for each fine */
628:   for (i=0;i<nl;i++) {
629:     pcontrib[i] = 0.;
630:     MatGetRow(lA,i,&ncols,&icol,NULL);
631:     if (gidx[i] >= fs && gidx[i] < fe) {
632:       li = gidx[i] - fs;
633:       lsparse[li] = 0;
634:       gsparse[li] = 0;
635:       cid = lcid[i];
636:       if (cid >= 0) {
637:         lsparse[li] = 1;
638:       } else {
639:         for (j=0;j<ncols;j++) {
640:           if (lcid[icol[j]] >= 0) {
641:             pcontrib[icol[j]] = 1.;
642:           } else {
643:             ci = icol[j];
644:             MatRestoreRow(lA,i,&ncols,&icol,NULL);
645:             MatGetRow(lA,ci,&ncols,&icol,NULL);
646:             for (k=0;k<ncols;k++) {
647:               if (lcid[icol[k]] >= 0) {
648:                 pcontrib[icol[k]] = 1.;
649:               }
650:             }
651:             MatRestoreRow(lA,ci,&ncols,&icol,NULL);
652:             MatGetRow(lA,i,&ncols,&icol,NULL);
653:           }
654:         }
655:         for (j=0;j<ncols;j++) {
656:           if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
657:             lni = lcid[icol[j]];
658:             if (lni >= cs && lni < ce) {
659:               lsparse[li]++;
660:             } else {
661:               gsparse[li]++;
662:             }
663:             pcontrib[icol[j]] = 0.;
664:           } else {
665:             ci = icol[j];
666:             MatRestoreRow(lA,i,&ncols,&icol,NULL);
667:             MatGetRow(lA,ci,&ncols,&icol,NULL);
668:             for (k=0;k<ncols;k++) {
669:               if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
670:                 lni = lcid[icol[k]];
671:                 if (lni >= cs && lni < ce) {
672:                   lsparse[li]++;
673:                 } else {
674:                   gsparse[li]++;
675:                 }
676:                 pcontrib[icol[k]] = 0.;
677:               }
678:             }
679:             MatRestoreRow(lA,ci,&ncols,&icol,NULL);
680:             MatGetRow(lA,i,&ncols,&icol,NULL);
681:           }
682:         }
683:       }
684:       if (lsparse[li] + gsparse[li] > maxcols) maxcols = lsparse[li]+gsparse[li];
685:     }
686:     MatRestoreRow(lA,i,&ncols,&icol,&vcol);
687:   }
688:   PetscMalloc2(maxcols,&picol,maxcols,&pvcol);
689:   MatCreate(PetscObjectComm((PetscObject)A),P);
690:   MatGetType(A,&mtype);
691:   MatSetType(*P,mtype);
692:   MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);
693:   MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);
694:   MatSeqAIJSetPreallocation(*P,0,lsparse);
695:   for (i=0;i<nl;i++) {
696:     diag = 0.;
697:     if (gidx[i] >= fs && gidx[i] < fe) {
698:       pncols=0;
699:       cid = lcid[i];
700:       if (cid >= 0) {
701:         pncols = 1;
702:         picol[0] = cid;
703:         pvcol[0] = 1.;
704:       } else {
705:         MatGetRow(lA,i,&ncols,&icol,&vcol);
706:         for (j=0;j<ncols;j++) {
707:           pentry = vcol[j];
708:           if (lcid[icol[j]] >= 0) {
709:             /* coarse neighbor */
710:             pcontrib[icol[j]] += pentry;
711:           } else if (icol[j] != i) {
712:             /* the neighbor is a strongly connected fine node */
713:             ci = icol[j];
714:             vi = vcol[j];
715:             MatRestoreRow(lA,i,&ncols,&icol,&vcol);
716:             MatGetRow(lA,ci,&ncols,&icol,&vcol);
717:             jwttotal=0.;
718:             jdiag = 0.;
719:             for (k=0;k<ncols;k++) {
720:               if (ci == icol[k]) {
721:                 jdiag = PetscRealPart(vcol[k]);
722:               }
723:             }
724:             for (k=0;k<ncols;k++) {
725:               if (lcid[icol[k]] >= 0 && jdiag*PetscRealPart(vcol[k]) < 0.) {
726:                 pjentry = vcol[k];
727:                 jwttotal += PetscRealPart(pjentry);
728:               }
729:             }
730:             if (jwttotal != 0.) {
731:               jwttotal = PetscRealPart(vi)/jwttotal;
732:               for (k=0;k<ncols;k++) {
733:                 if (lcid[icol[k]] >= 0 && jdiag*PetscRealPart(vcol[k]) < 0.) {
734:                   pjentry = vcol[k]*jwttotal;
735:                   pcontrib[icol[k]] += pjentry;
736:                 }
737:               }
738:             } else {
739:               diag += PetscRealPart(vi);
740:             }
741:             MatRestoreRow(lA,ci,&ncols,&icol,&vcol);
742:             MatGetRow(lA,i,&ncols,&icol,&vcol);
743:           } else {
744:             diag += PetscRealPart(vcol[j]);
745:           }
746:         }
747:         if (diag != 0.) {
748:           diag = 1./diag;
749:           for (j=0;j<ncols;j++) {
750:             if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
751:               /* the neighbor is a coarse node */
752:               if (PetscAbsScalar(pcontrib[icol[j]]) > 0.0) {
753:                 lni = lcid[icol[j]];
754:                 pvcol[pncols] = -pcontrib[icol[j]]*diag;
755:                 picol[pncols] = lni;
756:                 pncols++;
757:               }
758:               pcontrib[icol[j]] = 0.;
759:             } else {
760:               /* the neighbor is a strongly connected fine node */
761:               ci = icol[j];
762:               MatRestoreRow(lA,i,&ncols,&icol,&vcol);
763:               MatGetRow(lA,ci,&ncols,&icol,&vcol);
764:               for (k=0;k<ncols;k++) {
765:                 if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
766:                   if (PetscAbsScalar(pcontrib[icol[k]]) > 0.0) {
767:                     lni = lcid[icol[k]];
768:                     pvcol[pncols] = -pcontrib[icol[k]]*diag;
769:                     picol[pncols] = lni;
770:                     pncols++;
771:                   }
772:                   pcontrib[icol[k]] = 0.;
773:                 }
774:               }
775:               MatRestoreRow(lA,ci,&ncols,&icol,&vcol);
776:               MatGetRow(lA,i,&ncols,&icol,&vcol);
777:             }
778:             pcontrib[icol[j]] = 0.;
779:           }
780:           MatRestoreRow(lA,i,&ncols,&icol,&vcol);
781:         }
782:       }
783:       ci = gidx[i];
784:       if (pncols > 0) {
785:         MatSetValues(*P,1,&ci,pncols,picol,pvcol,INSERT_VALUES);
786:       }
787:     }
788:   }
789:   ISRestoreIndices(lis,&gidx);
790:   PetscFree2(picol,pvcol);
791:   PetscFree3(lsparse,gsparse,pcontrib);
792:   ISDestroy(&lis);
793:   PetscFree(gcid);
794:   if (size > 1) {
795:     PetscFree(lcid);
796:     MatDestroyMatrices(1,&lAs);
797:     PetscSFDestroy(&sf);
798:   }
799:   VecDestroy(&cv);
800:   MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
801:   MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);
802:   return 0;
803: }

805: PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc,Mat A,Mat *P)
806: {

808:   PetscInt          f,s,n,cf,cs,i,idx;
809:   PetscInt          *coarserows;
810:   PetscInt          ncols;
811:   const PetscInt    *pcols;
812:   const PetscScalar *pvals;
813:   Mat               Pnew;
814:   Vec               diag;
815:   PC_MG             *mg          = (PC_MG*)pc->data;
816:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
817:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

819:   if (cls->nsmooths == 0) {
820:     PCGAMGTruncateProlongator_Private(pc,P);
821:     return 0;
822:   }
823:   MatGetOwnershipRange(*P,&s,&f);
824:   n = f-s;
825:   MatGetOwnershipRangeColumn(*P,&cs,&cf);
826:   PetscMalloc1(n,&coarserows);
827:   /* identify the rows corresponding to coarse unknowns */
828:   idx = 0;
829:   for (i=s;i<f;i++) {
830:     MatGetRow(*P,i,&ncols,&pcols,&pvals);
831:     /* assume, for now, that it's a coarse unknown if it has a single unit entry */
832:     if (ncols == 1) {
833:       if (pvals[0] == 1.) {
834:         coarserows[idx] = i;
835:         idx++;
836:       }
837:     }
838:     MatRestoreRow(*P,i,&ncols,&pcols,&pvals);
839:   }
840:   MatCreateVecs(A,&diag,NULL);
841:   MatGetDiagonal(A,diag);
842:   VecReciprocal(diag);
843:   for (i=0;i<cls->nsmooths;i++) {
844:     MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew);
845:     MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL);
846:     MatDiagonalScale(Pnew,diag,NULL);
847:     MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN);
848:     MatDestroy(P);
849:     *P  = Pnew;
850:     Pnew = NULL;
851:   }
852:   VecDestroy(&diag);
853:   PetscFree(coarserows);
854:   PCGAMGTruncateProlongator_Private(pc,P);
855:   return 0;
856: }

858: PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
859: {
860:   PetscErrorCode    (*f)(PC,Mat,Mat,PetscCoarsenData*,Mat*);
861:   PC_MG             *mg          = (PC_MG*)pc->data;
862:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
863:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

865:   PetscFunctionListFind(PCGAMGClassicalProlongatorList,cls->prolongtype,&f);
867:   (*f)(pc,A,G,agg_lists,P);
868:   return 0;
869: }

871: PetscErrorCode PCGAMGDestroy_Classical(PC pc)
872: {
873:   PC_MG          *mg          = (PC_MG*)pc->data;
874:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;

876:   PetscFree(pc_gamg->subctx);
877:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",NULL);
878:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",NULL);
879:   return 0;
880: }

882: PetscErrorCode PCGAMGSetFromOptions_Classical(PetscOptionItems *PetscOptionsObject,PC pc)
883: {
884:   PC_MG             *mg          = (PC_MG*)pc->data;
885:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
886:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
887:   char              tname[256];
888:   PetscBool         flg;

890:   PetscOptionsHead(PetscOptionsObject,"GAMG-Classical options");
891:   PetscOptionsFList("-pc_gamg_classical_type","Type of Classical AMG prolongation","PCGAMGClassicalSetType",PCGAMGClassicalProlongatorList,cls->prolongtype, tname, sizeof(tname), &flg);
892:   if (flg) {
893:     PCGAMGClassicalSetType(pc,tname);
894:   }
895:   PetscOptionsReal("-pc_gamg_classical_interp_threshold","Threshold for classical interpolator entries","",cls->interp_threshold,&cls->interp_threshold,NULL);
896:   PetscOptionsInt("-pc_gamg_classical_nsmooths","Threshold for classical interpolator entries","",cls->nsmooths,&cls->nsmooths,NULL);
897:   PetscOptionsTail();
898:   return 0;
899: }

901: PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
902: {
903:   PC_MG          *mg      = (PC_MG*)pc->data;
904:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;

906:   /* no data for classical AMG */
907:   pc_gamg->data           = NULL;
908:   pc_gamg->data_cell_cols = 0;
909:   pc_gamg->data_cell_rows = 0;
910:   pc_gamg->data_sz        = 0;
911:   return 0;
912: }

914: PetscErrorCode PCGAMGClassicalFinalizePackage(void)
915: {
916:   PCGAMGClassicalPackageInitialized = PETSC_FALSE;
917:   PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList);
918:   return 0;
919: }

921: PetscErrorCode PCGAMGClassicalInitializePackage(void)
922: {
923:   if (PCGAMGClassicalPackageInitialized) return 0;
924:   PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALDIRECT,PCGAMGProlongator_Classical_Direct);
925:   PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALSTANDARD,PCGAMGProlongator_Classical_Standard);
926:   PetscRegisterFinalize(PCGAMGClassicalFinalizePackage);
927:   return 0;
928: }

930: /* -------------------------------------------------------------------------- */
931: /*
932:    PCCreateGAMG_Classical

934: */
935: PetscErrorCode  PCCreateGAMG_Classical(PC pc)
936: {
937:   PC_MG             *mg      = (PC_MG*)pc->data;
938:   PC_GAMG           *pc_gamg = (PC_GAMG*)mg->innerctx;
939:   PC_GAMG_Classical *pc_gamg_classical;

941:   PCGAMGClassicalInitializePackage();
942:   if (pc_gamg->subctx) {
943:     /* call base class */
944:     PCDestroy_GAMG(pc);
945:   }

947:   /* create sub context for SA */
948:   PetscNewLog(pc,&pc_gamg_classical);
949:   pc_gamg->subctx = pc_gamg_classical;
950:   pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
951:   /* reset does not do anything; setup not virtual */

953:   /* set internal function pointers */
954:   pc_gamg->ops->destroy        = PCGAMGDestroy_Classical;
955:   pc_gamg->ops->graph          = PCGAMGGraph_Classical;
956:   pc_gamg->ops->coarsen        = PCGAMGCoarsen_Classical;
957:   pc_gamg->ops->prolongator    = PCGAMGProlongator_Classical;
958:   pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
959:   pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;

961:   pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
962:   pc_gamg_classical->interp_threshold = 0.2;
963:   pc_gamg_classical->nsmooths         = 0;
964:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",PCGAMGClassicalSetType_GAMG);
965:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",PCGAMGClassicalGetType_GAMG);
966:   PCGAMGClassicalSetType(pc,PCGAMGCLASSICALSTANDARD);
967:   return 0;
968: }