Actual source code: classical.c

petsc-master 2019-08-22
Report Typos and Errors
  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: }


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

880: PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
881: {
882:   PetscErrorCode    ierr;
883:   PetscErrorCode    (*f)(PC,Mat,Mat,PetscCoarsenData*,Mat*);
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;

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

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

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

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

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

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

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


944: PetscErrorCode PCGAMGClassicalFinalizePackage(void)
945: {

949:   PCGAMGClassicalPackageInitialized = PETSC_FALSE;
950:   PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList);
951:   return(0);
952: }

954: PetscErrorCode PCGAMGClassicalInitializePackage(void)
955: {

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

966: /* -------------------------------------------------------------------------- */
967: /*
968:    PCCreateGAMG_Classical

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

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

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

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

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