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

 24:    Level: intermediate

 26: .seealso: ()
 27: @*/
 28: PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type)
 29: {

 34:   PetscTryMethod(pc,"PCGAMGClassicalSetType_C",(PC,PCGAMGClassicalType),(pc,type));
 35:   return(0);
 36: }

 38: /*@C
 39:    PCGAMGClassicalGetType - Gets the type of classical interpolation to use

 41:    Collective on PC

 43:    Input Parameter:
 44: .  pc - the preconditioner context

 46:    Output Parameter:
 47: .  type - the type used

 49:    Level: intermediate

 51: .seealso: ()
 52: @*/
 53: PetscErrorCode PCGAMGClassicalGetType(PC pc, PCGAMGClassicalType *type)
 54: {

 59:   PetscUseMethod(pc,"PCGAMGClassicalGetType_C",(PC,PCGAMGClassicalType*),(pc,type));
 60:   return(0);
 61: }

 63: static PetscErrorCode PCGAMGClassicalSetType_GAMG(PC pc, PCGAMGClassicalType type)
 64: {
 65:   PetscErrorCode    ierr;
 66:   PC_MG             *mg          = (PC_MG*)pc->data;
 67:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
 68:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

 71:   PetscStrcpy(cls->prolongtype,type);
 72:   return(0);
 73: }

 75: static PetscErrorCode PCGAMGClassicalGetType_GAMG(PC pc, PCGAMGClassicalType *type)
 76: {
 77:   PC_MG             *mg          = (PC_MG*)pc->data;
 78:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
 79:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

 82:   *type = cls->prolongtype;
 83:   return(0);
 84: }

 86: PetscErrorCode PCGAMGGraph_Classical(PC pc,Mat A,Mat *G)
 87: {
 88:   PetscInt          s,f,n,idx,lidx,gidx;
 89:   PetscInt          r,c,ncols;
 90:   const PetscInt    *rcol;
 91:   const PetscScalar *rval;
 92:   PetscInt          *gcol;
 93:   PetscScalar       *gval;
 94:   PetscReal         rmax;
 95:   PetscInt          cmax = 0;
 96:   PC_MG             *mg = (PC_MG *)pc->data;
 97:   PC_GAMG           *gamg = (PC_GAMG *)mg->innerctx;
 98:   PetscErrorCode    ierr;
 99:   PetscInt          *gsparse,*lsparse;
100:   PetscScalar       *Amax;
101:   MatType           mtype;

104:   MatGetOwnershipRange(A,&s,&f);
105:   n=f-s;
106:   PetscMalloc3(n,&lsparse,n,&gsparse,n,&Amax);

108:   for (r = 0;r < n;r++) {
109:     lsparse[r] = 0;
110:     gsparse[r] = 0;
111:   }

113:   for (r = s;r < f;r++) {
114:     /* determine the maximum off-diagonal in each row */
115:     rmax = 0.;
116:     MatGetRow(A,r,&ncols,&rcol,&rval);
117:     for (c = 0; c < ncols; c++) {
118:       if (PetscRealPart(-rval[c]) > rmax && rcol[c] != r) {
119:         rmax = PetscRealPart(-rval[c]);
120:       }
121:     }
122:     Amax[r-s] = rmax;
123:     if (ncols > cmax) cmax = ncols;
124:     lidx = 0;
125:     gidx = 0;
126:     /* create the local and global sparsity patterns */
127:     for (c = 0; c < ncols; c++) {
128:       if (PetscRealPart(-rval[c]) > gamg->threshold[0]*PetscRealPart(Amax[r-s]) || rcol[c] == r) {
129:         if (rcol[c] < f && rcol[c] >= s) {
130:           lidx++;
131:         } else {
132:           gidx++;
133:         }
134:       }
135:     }
136:     MatRestoreRow(A,r,&ncols,&rcol,&rval);
137:     lsparse[r-s] = lidx;
138:     gsparse[r-s] = gidx;
139:   }
140:   PetscMalloc2(cmax,&gval,cmax,&gcol);

142:   MatCreate(PetscObjectComm((PetscObject)A),G);
143:   MatGetType(A,&mtype);
144:   MatSetType(*G,mtype);
145:   MatSetSizes(*G,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
146:   MatMPIAIJSetPreallocation(*G,0,lsparse,0,gsparse);
147:   MatSeqAIJSetPreallocation(*G,0,lsparse);
148:   for (r = s;r < f;r++) {
149:     MatGetRow(A,r,&ncols,&rcol,&rval);
150:     idx = 0;
151:     for (c = 0; c < ncols; c++) {
152:       /* classical strength of connection */
153:       if (PetscRealPart(-rval[c]) > gamg->threshold[0]*PetscRealPart(Amax[r-s]) || rcol[c] == r) {
154:         gcol[idx] = rcol[c];
155:         gval[idx] = rval[c];
156:         idx++;
157:       }
158:     }
159:     MatSetValues(*G,1,&r,idx,gcol,gval,INSERT_VALUES);
160:     MatRestoreRow(A,r,&ncols,&rcol,&rval);
161:   }
162:   MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY);
163:   MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY);

165:   PetscFree2(gval,gcol);
166:   PetscFree3(lsparse,gsparse,Amax);
167:   return(0);
168: }

170: PetscErrorCode PCGAMGCoarsen_Classical(PC pc,Mat *G,PetscCoarsenData **agg_lists)
171: {
172:   PetscErrorCode   ierr;
173:   MatCoarsen       crs;
174:   MPI_Comm         fcomm = ((PetscObject)pc)->comm;

177:   if (!G) SETERRQ(fcomm,PETSC_ERR_ARG_WRONGSTATE,"Must set Graph in PC in PCGAMG before coarsening");

179:   MatCoarsenCreate(fcomm,&crs);
180:   MatCoarsenSetFromOptions(crs);
181:   MatCoarsenSetAdjacency(crs,*G);
182:   MatCoarsenSetStrictAggs(crs,PETSC_TRUE);
183:   MatCoarsenApply(crs);
184:   MatCoarsenGetData(crs,agg_lists);
185:   MatCoarsenDestroy(&crs);
186:   return(0);
187: }

189: PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
190: {
191:   PetscErrorCode    ierr;
192:   PC_MG             *mg          = (PC_MG*)pc->data;
193:   PC_GAMG           *gamg        = (PC_GAMG*)mg->innerctx;
194:   PetscBool         iscoarse,isMPIAIJ,isSEQAIJ;
195:   PetscInt          fn,cn,fs,fe,cs,ce,i,j,ncols,col,row_f,row_c,cmax=0,idx,noff;
196:   PetscInt          *lcid,*gcid,*lsparse,*gsparse,*colmap,*pcols;
197:   const PetscInt    *rcol;
198:   PetscReal         *Amax_pos,*Amax_neg;
199:   PetscScalar       g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta,pij;
200:   PetscScalar       *pvals;
201:   const PetscScalar *rval;
202:   Mat               lA,gA=NULL;
203:   MatType           mtype;
204:   Vec               C,lvec;
205:   PetscLayout       clayout;
206:   PetscSF           sf;
207:   Mat_MPIAIJ        *mpiaij;

210:   MatGetOwnershipRange(A,&fs,&fe);
211:   fn = fe-fs;
212:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
213:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSEQAIJ);
214:   if (!isMPIAIJ && !isSEQAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Classical AMG requires MPIAIJ matrix");
215:   if (isMPIAIJ) {
216:     mpiaij = (Mat_MPIAIJ*)A->data;
217:     lA = mpiaij->A;
218:     gA = mpiaij->B;
219:     lvec = mpiaij->lvec;
220:     VecGetSize(lvec,&noff);
221:     colmap = mpiaij->garray;
222:     MatGetLayouts(A,NULL,&clayout);
223:     PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
224:     PetscSFSetGraphLayout(sf,clayout,noff,NULL,PETSC_COPY_VALUES,colmap);
225:     PetscMalloc1(noff,&gcid);
226:   } else {
227:     lA = A;
228:   }
229:   PetscMalloc5(fn,&lsparse,fn,&gsparse,fn,&lcid,fn,&Amax_pos,fn,&Amax_neg);

231:   /* count the number of coarse unknowns */
232:   cn = 0;
233:   for (i=0;i<fn;i++) {
234:     /* filter out singletons */
235:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
236:     lcid[i] = -1;
237:     if (!iscoarse) {
238:       cn++;
239:     }
240:   }

242:    /* create the coarse vector */
243:   VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&C);
244:   VecGetOwnershipRange(C,&cs,&ce);

246:   cn = 0;
247:   for (i=0;i<fn;i++) {
248:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
249:     if (!iscoarse) {
250:       lcid[i] = cs+cn;
251:       cn++;
252:     } else {
253:       lcid[i] = -1;
254:     }
255:   }

257:   if (gA) {
258:     PetscSFBcastBegin(sf,MPIU_INT,lcid,gcid,MPI_REPLACE);
259:     PetscSFBcastEnd(sf,MPIU_INT,lcid,gcid,MPI_REPLACE);
260:   }

262:   /* determine the largest off-diagonal entries in each row */
263:   for (i=fs;i<fe;i++) {
264:     Amax_pos[i-fs] = 0.;
265:     Amax_neg[i-fs] = 0.;
266:     MatGetRow(A,i,&ncols,&rcol,&rval);
267:     for (j=0;j<ncols;j++) {
268:       if ((PetscRealPart(-rval[j]) > Amax_neg[i-fs]) && i != rcol[j]) Amax_neg[i-fs] = PetscAbsScalar(rval[j]);
269:       if ((PetscRealPart(rval[j])  > Amax_pos[i-fs]) && i != rcol[j]) Amax_pos[i-fs] = PetscAbsScalar(rval[j]);
270:     }
271:     if (ncols > cmax) cmax = ncols;
272:     MatRestoreRow(A,i,&ncols,&rcol,&rval);
273:   }
274:   PetscMalloc2(cmax,&pcols,cmax,&pvals);
275:   VecDestroy(&C);

277:   /* count the on and off processor sparsity patterns for the prolongator */
278:   for (i=0;i<fn;i++) {
279:     /* on */
280:     lsparse[i] = 0;
281:     gsparse[i] = 0;
282:     if (lcid[i] >= 0) {
283:       lsparse[i] = 1;
284:       gsparse[i] = 0;
285:     } else {
286:       MatGetRow(lA,i,&ncols,&rcol,&rval);
287:       for (j = 0;j < ncols;j++) {
288:         col = rcol[j];
289:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
290:           lsparse[i] += 1;
291:         }
292:       }
293:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);
294:       /* off */
295:       if (gA) {
296:         MatGetRow(gA,i,&ncols,&rcol,&rval);
297:         for (j = 0; j < ncols; j++) {
298:           col = rcol[j];
299:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
300:             gsparse[i] += 1;
301:           }
302:         }
303:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
304:       }
305:     }
306:   }

308:   /* preallocate and create the prolongator */
309:   MatCreate(PetscObjectComm((PetscObject)A),P);
310:   MatGetType(G,&mtype);
311:   MatSetType(*P,mtype);
312:   MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);
313:   MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);
314:   MatSeqAIJSetPreallocation(*P,0,lsparse);

316:   /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
317:   for (i = 0;i < fn;i++) {
318:     /* determine on or off */
319:     row_f = i + fs;
320:     row_c = lcid[i];
321:     if (row_c >= 0) {
322:       pij = 1.;
323:       MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);
324:     } else {
325:       g_pos = 0.;
326:       g_neg = 0.;
327:       a_pos = 0.;
328:       a_neg = 0.;
329:       diag  = 0.;

331:       /* local connections */
332:       MatGetRow(lA,i,&ncols,&rcol,&rval);
333:       for (j = 0; j < ncols; j++) {
334:         col = rcol[j];
335:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
336:           if (PetscRealPart(rval[j]) > 0.) {
337:             g_pos += rval[j];
338:           } else {
339:             g_neg += rval[j];
340:           }
341:         }
342:         if (col != i) {
343:           if (PetscRealPart(rval[j]) > 0.) {
344:             a_pos += rval[j];
345:           } else {
346:             a_neg += rval[j];
347:           }
348:         } else {
349:           diag = rval[j];
350:         }
351:       }
352:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);

354:       /* ghosted connections */
355:       if (gA) {
356:         MatGetRow(gA,i,&ncols,&rcol,&rval);
357:         for (j = 0; j < ncols; j++) {
358:           col = rcol[j];
359:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
360:             if (PetscRealPart(rval[j]) > 0.) {
361:               g_pos += rval[j];
362:             } else {
363:               g_neg += rval[j];
364:             }
365:           }
366:           if (PetscRealPart(rval[j]) > 0.) {
367:             a_pos += rval[j];
368:           } else {
369:             a_neg += rval[j];
370:           }
371:         }
372:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
373:       }

375:       if (g_neg == 0.) {
376:         alpha = 0.;
377:       } else {
378:         alpha = -a_neg/g_neg;
379:       }

381:       if (g_pos == 0.) {
382:         diag += a_pos;
383:         beta = 0.;
384:       } else {
385:         beta = -a_pos/g_pos;
386:       }
387:       if (diag == 0.) {
388:         invdiag = 0.;
389:       } else invdiag = 1. / diag;
390:       /* on */
391:       MatGetRow(lA,i,&ncols,&rcol,&rval);
392:       idx = 0;
393:       for (j = 0;j < ncols;j++) {
394:         col = rcol[j];
395:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
396:           row_f = i + fs;
397:           row_c = lcid[col];
398:           /* set the values for on-processor ones */
399:           if (PetscRealPart(rval[j]) < 0.) {
400:             pij = rval[j]*alpha*invdiag;
401:           } else {
402:             pij = rval[j]*beta*invdiag;
403:           }
404:           if (PetscAbsScalar(pij) != 0.) {
405:             pvals[idx] = pij;
406:             pcols[idx] = row_c;
407:             idx++;
408:           }
409:         }
410:       }
411:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);
412:       /* off */
413:       if (gA) {
414:         MatGetRow(gA,i,&ncols,&rcol,&rval);
415:         for (j = 0; j < ncols; j++) {
416:           col = rcol[j];
417:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
418:             row_f = i + fs;
419:             row_c = gcid[col];
420:             /* set the values for on-processor ones */
421:             if (PetscRealPart(rval[j]) < 0.) {
422:               pij = rval[j]*alpha*invdiag;
423:             } else {
424:               pij = rval[j]*beta*invdiag;
425:             }
426:             if (PetscAbsScalar(pij) != 0.) {
427:               pvals[idx] = pij;
428:               pcols[idx] = row_c;
429:               idx++;
430:             }
431:           }
432:         }
433:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
434:       }
435:       MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);
436:     }
437:   }

439:   MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
440:   MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);

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

444:   PetscFree2(pcols,pvals);
445:   if (gA) {
446:     PetscSFDestroy(&sf);
447:     PetscFree(gcid);
448:   }
449:   return(0);
450: }

452: PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc,Mat *P)
453: {
454:   PetscInt          j,i,ps,pf,pn,pcs,pcf,pcn,idx,cmax;
455:   PetscErrorCode    ierr;
456:   const PetscScalar *pval;
457:   const PetscInt    *pcol;
458:   PetscScalar       *pnval;
459:   PetscInt          *pncol;
460:   PetscInt          ncols;
461:   Mat               Pnew;
462:   PetscInt          *lsparse,*gsparse;
463:   PetscReal         pmax_pos,pmax_neg,ptot_pos,ptot_neg,pthresh_pos,pthresh_neg;
464:   PC_MG             *mg          = (PC_MG*)pc->data;
465:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
466:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
467:   MatType           mtype;

470:   /* trim and rescale with reallocation */
471:   MatGetOwnershipRange(*P,&ps,&pf);
472:   MatGetOwnershipRangeColumn(*P,&pcs,&pcf);
473:   pn = pf-ps;
474:   pcn = pcf-pcs;
475:   PetscMalloc2(pn,&lsparse,pn,&gsparse);
476:   /* allocate */
477:   cmax = 0;
478:   for (i=ps;i<pf;i++) {
479:     lsparse[i-ps] = 0;
480:     gsparse[i-ps] = 0;
481:     MatGetRow(*P,i,&ncols,&pcol,&pval);
482:     if (ncols > cmax) {
483:       cmax = ncols;
484:     }
485:     pmax_pos = 0.;
486:     pmax_neg = 0.;
487:     for (j=0;j<ncols;j++) {
488:       if (PetscRealPart(pval[j]) > pmax_pos) {
489:         pmax_pos = PetscRealPart(pval[j]);
490:       } else if (PetscRealPart(pval[j]) < pmax_neg) {
491:         pmax_neg = PetscRealPart(pval[j]);
492:       }
493:     }
494:     for (j=0;j<ncols;j++) {
495:       if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
496:         if (pcol[j] >= pcs && pcol[j] < pcf) {
497:           lsparse[i-ps]++;
498:         } else {
499:           gsparse[i-ps]++;
500:         }
501:       }
502:     }
503:     MatRestoreRow(*P,i,&ncols,&pcol,&pval);
504:   }

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

508:   MatGetType(*P,&mtype);
509:   MatCreate(PetscObjectComm((PetscObject)*P),&Pnew);
510:   MatSetType(Pnew, mtype);
511:   MatSetSizes(Pnew,pn,pcn,PETSC_DETERMINE,PETSC_DETERMINE);
512:   MatSeqAIJSetPreallocation(Pnew,0,lsparse);
513:   MatMPIAIJSetPreallocation(Pnew,0,lsparse,0,gsparse);

515:   for (i=ps;i<pf;i++) {
516:     MatGetRow(*P,i,&ncols,&pcol,&pval);
517:     pmax_pos = 0.;
518:     pmax_neg = 0.;
519:     for (j=0;j<ncols;j++) {
520:       if (PetscRealPart(pval[j]) > pmax_pos) {
521:         pmax_pos = PetscRealPart(pval[j]);
522:       } else if (PetscRealPart(pval[j]) < pmax_neg) {
523:         pmax_neg = PetscRealPart(pval[j]);
524:       }
525:     }
526:     pthresh_pos = 0.;
527:     pthresh_neg = 0.;
528:     ptot_pos = 0.;
529:     ptot_neg = 0.;
530:     for (j=0;j<ncols;j++) {
531:       if (PetscRealPart(pval[j]) >= cls->interp_threshold*pmax_pos) {
532:         pthresh_pos += PetscRealPart(pval[j]);
533:       } else if (PetscRealPart(pval[j]) <= cls->interp_threshold*pmax_neg) {
534:         pthresh_neg += PetscRealPart(pval[j]);
535:       }
536:       if (PetscRealPart(pval[j]) > 0.) {
537:         ptot_pos += PetscRealPart(pval[j]);
538:       } else {
539:         ptot_neg += PetscRealPart(pval[j]);
540:       }
541:     }
542:     if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
543:     if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
544:     idx=0;
545:     for (j=0;j<ncols;j++) {
546:       if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold) {
547:         pnval[idx] = ptot_pos*pval[j];
548:         pncol[idx] = pcol[j];
549:         idx++;
550:       } else if (PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
551:         pnval[idx] = ptot_neg*pval[j];
552:         pncol[idx] = pcol[j];
553:         idx++;
554:       }
555:     }
556:     MatRestoreRow(*P,i,&ncols,&pcol,&pval);
557:     MatSetValues(Pnew,1,&i,idx,pncol,pnval,INSERT_VALUES);
558:   }

560:   MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY);
561:   MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY);
562:   MatDestroy(P);

564:   *P = Pnew;
565:   PetscFree2(lsparse,gsparse);
566:   PetscFree2(pnval,pncol);
567:   return(0);
568: }

570: PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
571: {
572:   PetscErrorCode    ierr;
573:   Mat               lA,*lAs;
574:   MatType           mtype;
575:   Vec               cv;
576:   PetscInt          *gcid,*lcid,*lsparse,*gsparse,*picol;
577:   PetscInt          fs,fe,cs,ce,nl,i,j,k,li,lni,ci,ncols,maxcols,fn,cn,cid;
578:   PetscMPIInt       size;
579:   const PetscInt    *lidx,*icol,*gidx;
580:   PetscBool         iscoarse;
581:   PetscScalar       vi,pentry,pjentry;
582:   PetscScalar       *pcontrib,*pvcol;
583:   const PetscScalar *vcol;
584:   PetscReal         diag,jdiag,jwttotal;
585:   PetscInt          pncols;
586:   PetscSF           sf;
587:   PetscLayout       clayout;
588:   IS                lis;

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

824: PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc,Mat A,Mat *P)
825: {

827:   PetscErrorCode    ierr;
828:   PetscInt          f,s,n,cf,cs,i,idx;
829:   PetscInt          *coarserows;
830:   PetscInt          ncols;
831:   const PetscInt    *pcols;
832:   const PetscScalar *pvals;
833:   Mat               Pnew;
834:   Vec               diag;
835:   PC_MG             *mg          = (PC_MG*)pc->data;
836:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
837:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

840:   if (cls->nsmooths == 0) {
841:     PCGAMGTruncateProlongator_Private(pc,P);
842:     return(0);
843:   }
844:   MatGetOwnershipRange(*P,&s,&f);
845:   n = f-s;
846:   MatGetOwnershipRangeColumn(*P,&cs,&cf);
847:   PetscMalloc1(n,&coarserows);
848:   /* identify the rows corresponding to coarse unknowns */
849:   idx = 0;
850:   for (i=s;i<f;i++) {
851:     MatGetRow(*P,i,&ncols,&pcols,&pvals);
852:     /* assume, for now, that it's a coarse unknown if it has a single unit entry */
853:     if (ncols == 1) {
854:       if (pvals[0] == 1.) {
855:         coarserows[idx] = i;
856:         idx++;
857:       }
858:     }
859:     MatRestoreRow(*P,i,&ncols,&pcols,&pvals);
860:   }
861:   MatCreateVecs(A,&diag,NULL);
862:   MatGetDiagonal(A,diag);
863:   VecReciprocal(diag);
864:   for (i=0;i<cls->nsmooths;i++) {
865:     MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew);
866:     MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL);
867:     MatDiagonalScale(Pnew,diag,NULL);
868:     MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN);
869:     MatDestroy(P);
870:     *P  = Pnew;
871:     Pnew = NULL;
872:   }
873:   VecDestroy(&diag);
874:   PetscFree(coarserows);
875:   PCGAMGTruncateProlongator_Private(pc,P);
876:   return(0);
877: }

879: PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
880: {
881:   PetscErrorCode    ierr;
882:   PetscErrorCode    (*f)(PC,Mat,Mat,PetscCoarsenData*,Mat*);
883:   PC_MG             *mg          = (PC_MG*)pc->data;
884:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
885:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

888:   PetscFunctionListFind(PCGAMGClassicalProlongatorList,cls->prolongtype,&f);
889:   if (!f)SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot find PCGAMG Classical prolongator type");
890:   (*f)(pc,A,G,agg_lists,P);
891:   return(0);
892: }

894: PetscErrorCode PCGAMGDestroy_Classical(PC pc)
895: {
897:   PC_MG          *mg          = (PC_MG*)pc->data;
898:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;

901:   PetscFree(pc_gamg->subctx);
902:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",NULL);
903:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",NULL);
904:   return(0);
905: }

907: PetscErrorCode PCGAMGSetFromOptions_Classical(PetscOptionItems *PetscOptionsObject,PC pc)
908: {
909:   PC_MG             *mg          = (PC_MG*)pc->data;
910:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
911:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
912:   char              tname[256];
913:   PetscErrorCode    ierr;
914:   PetscBool         flg;

917:   PetscOptionsHead(PetscOptionsObject,"GAMG-Classical options");
918:   PetscOptionsFList("-pc_gamg_classical_type","Type of Classical AMG prolongation","PCGAMGClassicalSetType",PCGAMGClassicalProlongatorList,cls->prolongtype, tname, sizeof(tname), &flg);
919:   if (flg) {
920:     PCGAMGClassicalSetType(pc,tname);
921:   }
922:   PetscOptionsReal("-pc_gamg_classical_interp_threshold","Threshold for classical interpolator entries","",cls->interp_threshold,&cls->interp_threshold,NULL);
923:   PetscOptionsInt("-pc_gamg_classical_nsmooths","Threshold for classical interpolator entries","",cls->nsmooths,&cls->nsmooths,NULL);
924:   PetscOptionsTail();
925:   return(0);
926: }

928: PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
929: {
930:   PC_MG          *mg      = (PC_MG*)pc->data;
931:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;

934:   /* no data for classical AMG */
935:   pc_gamg->data           = NULL;
936:   pc_gamg->data_cell_cols = 0;
937:   pc_gamg->data_cell_rows = 0;
938:   pc_gamg->data_sz        = 0;
939:   return(0);
940: }

942: PetscErrorCode PCGAMGClassicalFinalizePackage(void)
943: {

947:   PCGAMGClassicalPackageInitialized = PETSC_FALSE;
948:   PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList);
949:   return(0);
950: }

952: PetscErrorCode PCGAMGClassicalInitializePackage(void)
953: {

957:   if (PCGAMGClassicalPackageInitialized) return(0);
958:   PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALDIRECT,PCGAMGProlongator_Classical_Direct);
959:   PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALSTANDARD,PCGAMGProlongator_Classical_Standard);
960:   PetscRegisterFinalize(PCGAMGClassicalFinalizePackage);
961:   return(0);
962: }

964: /* -------------------------------------------------------------------------- */
965: /*
966:    PCCreateGAMG_Classical

968: */
969: PetscErrorCode  PCCreateGAMG_Classical(PC pc)
970: {
972:   PC_MG             *mg      = (PC_MG*)pc->data;
973:   PC_GAMG           *pc_gamg = (PC_GAMG*)mg->innerctx;
974:   PC_GAMG_Classical *pc_gamg_classical;

977:   PCGAMGClassicalInitializePackage();
978:   if (pc_gamg->subctx) {
979:     /* call base class */
980:     PCDestroy_GAMG(pc);
981:   }

983:   /* create sub context for SA */
984:   PetscNewLog(pc,&pc_gamg_classical);
985:   pc_gamg->subctx = pc_gamg_classical;
986:   pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
987:   /* reset does not do anything; setup not virtual */

989:   /* set internal function pointers */
990:   pc_gamg->ops->destroy        = PCGAMGDestroy_Classical;
991:   pc_gamg->ops->graph          = PCGAMGGraph_Classical;
992:   pc_gamg->ops->coarsen        = PCGAMGCoarsen_Classical;
993:   pc_gamg->ops->prolongator    = PCGAMGProlongator_Classical;
994:   pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
995:   pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;

997:   pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
998:   pc_gamg_classical->interp_threshold = 0.2;
999:   pc_gamg_classical->nsmooths         = 0;
1000:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",PCGAMGClassicalSetType_GAMG);
1001:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",PCGAMGClassicalGetType_GAMG);
1002:   PCGAMGClassicalSetType(pc,PCGAMGCLASSICALSTANDARD);
1003:   return(0);
1004: }