Actual source code: classical.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <../src/ksp/pc/impls/gamg/gamg.h>        /*I "petscpc.h" I*/
  2: #include <petsc/private/kspimpl.h>
  3: #include <petscsf.h>

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

  8: typedef struct {
  9:   PetscReal interp_threshold; /* interpolation threshold */
 10:   char      prolongtype[256];
 11:   PetscInt  nsmooths;         /* number of jacobi smoothings on the prolongator */
 12: } PC_GAMG_Classical;

 16: /*@C
 17:    PCGAMGClassicalSetType - Sets the type of classical interpolation to use

 19:    Collective on PC

 21:    Input Parameters:
 22: .  pc - the preconditioner context

 24:    Options Database Key:
 25: .  -pc_gamg_classical_type

 27:    Level: intermediate

 29: .seealso: ()
 30: @*/
 31: PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type)
 32: {

 37:   PetscTryMethod(pc,"PCGAMGClassicalSetType_C",(PC,PCGAMGType),(pc,type));
 38:   return(0);
 39: }

 43: /*@C
 44:    PCGAMGClassicalGetType - Gets the type of classical interpolation to use

 46:    Collective on PC

 48:    Input Parameter:
 49: .  pc - the preconditioner context

 51:    Output Parameter:
 52: .  type - the type used

 54:    Level: intermediate

 56: .seealso: ()
 57: @*/
 58: PetscErrorCode PCGAMGClassicalGetType(PC pc, PCGAMGClassicalType *type)
 59: {

 64:   PetscUseMethod(pc,"PCGAMGClassicalGetType_C",(PC,PCGAMGType*),(pc,type));
 65:   return(0);
 66: }

 70: static PetscErrorCode PCGAMGClassicalSetType_GAMG(PC pc, PCGAMGClassicalType type)
 71: {
 72:   PetscErrorCode    ierr;
 73:   PC_MG             *mg          = (PC_MG*)pc->data;
 74:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
 75:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

 78:   PetscStrcpy(cls->prolongtype,type);
 79:   return(0);
 80: }

 84: static PetscErrorCode PCGAMGClassicalGetType_GAMG(PC pc, PCGAMGClassicalType *type)
 85: {
 86:   PC_MG             *mg          = (PC_MG*)pc->data;
 87:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
 88:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

 91:   *type = cls->prolongtype;
 92:   return(0);
 93: }

 97: PetscErrorCode PCGAMGGraph_Classical(PC pc,Mat A,Mat *G)
 98: {
 99:   PetscInt          s,f,n,idx,lidx,gidx;
100:   PetscInt          r,c,ncols;
101:   const PetscInt    *rcol;
102:   const PetscScalar *rval;
103:   PetscInt          *gcol;
104:   PetscScalar       *gval;
105:   PetscReal         rmax;
106:   PetscInt          cmax = 0;
107:   PC_MG             *mg;
108:   PC_GAMG           *gamg;
109:   PetscErrorCode    ierr;
110:   PetscInt          *gsparse,*lsparse;
111:   PetscScalar       *Amax;
112:   MatType           mtype;

115:   mg   = (PC_MG *)pc->data;
116:   gamg = (PC_GAMG *)mg->innerctx;

118:   MatGetOwnershipRange(A,&s,&f);
119:   n=f-s;
120:   PetscMalloc3(n,&lsparse,n,&gsparse,n,&Amax);

122:   for (r = 0;r < n;r++) {
123:     lsparse[r] = 0;
124:     gsparse[r] = 0;
125:   }

127:   for (r = s;r < f;r++) {
128:     /* determine the maximum off-diagonal in each row */
129:     rmax = 0.;
130:     MatGetRow(A,r,&ncols,&rcol,&rval);
131:     for (c = 0; c < ncols; c++) {
132:       if (PetscRealPart(-rval[c]) > rmax && rcol[c] != r) {
133:         rmax = PetscRealPart(-rval[c]);
134:       }
135:     }
136:     Amax[r-s] = rmax;
137:     if (ncols > cmax) cmax = ncols;
138:     lidx = 0;
139:     gidx = 0;
140:     /* create the local and global sparsity patterns */
141:     for (c = 0; c < ncols; c++) {
142:       if (PetscRealPart(-rval[c]) > gamg->threshold*PetscRealPart(Amax[r-s]) || rcol[c] == r) {
143:         if (rcol[c] < f && rcol[c] >= s) {
144:           lidx++;
145:         } else {
146:           gidx++;
147:         }
148:       }
149:     }
150:     MatRestoreRow(A,r,&ncols,&rcol,&rval);
151:     lsparse[r-s] = lidx;
152:     gsparse[r-s] = gidx;
153:   }
154:   PetscMalloc2(cmax,&gval,cmax,&gcol);

156:   MatCreate(PetscObjectComm((PetscObject)A),G);
157:   MatGetType(A,&mtype);
158:   MatSetType(*G,mtype);
159:   MatSetSizes(*G,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
160:   MatMPIAIJSetPreallocation(*G,0,lsparse,0,gsparse);
161:   MatSeqAIJSetPreallocation(*G,0,lsparse);
162:   for (r = s;r < f;r++) {
163:     MatGetRow(A,r,&ncols,&rcol,&rval);
164:     idx = 0;
165:     for (c = 0; c < ncols; c++) {
166:       /* classical strength of connection */
167:       if (PetscRealPart(-rval[c]) > gamg->threshold*PetscRealPart(Amax[r-s]) || rcol[c] == r) {
168:         gcol[idx] = rcol[c];
169:         gval[idx] = rval[c];
170:         idx++;
171:       }
172:     }
173:     MatSetValues(*G,1,&r,idx,gcol,gval,INSERT_VALUES);
174:     MatRestoreRow(A,r,&ncols,&rcol,&rval);
175:   }
176:   MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY);
177:   MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY);

179:   PetscFree2(gval,gcol);
180:   PetscFree3(lsparse,gsparse,Amax);
181:   return(0);
182: }


187: PetscErrorCode PCGAMGCoarsen_Classical(PC pc,Mat *G,PetscCoarsenData **agg_lists)
188: {
189:   PetscErrorCode   ierr;
190:   MatCoarsen       crs;
191:   MPI_Comm         fcomm = ((PetscObject)pc)->comm;

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

196:   MatCoarsenCreate(fcomm,&crs);
197:   MatCoarsenSetFromOptions(crs);
198:   MatCoarsenSetAdjacency(crs,*G);
199:   MatCoarsenSetStrictAggs(crs,PETSC_TRUE);
200:   MatCoarsenApply(crs);
201:   MatCoarsenGetData(crs,agg_lists);
202:   MatCoarsenDestroy(&crs);

204:   return(0);
205: }

209: PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
210: {
211:   PetscErrorCode    ierr;
212:   PC_MG             *mg          = (PC_MG*)pc->data;
213:   PC_GAMG           *gamg        = (PC_GAMG*)mg->innerctx;
214:   PetscBool         iscoarse,isMPIAIJ,isSEQAIJ;
215:   PetscInt          fn,cn,fs,fe,cs,ce,i,j,ncols,col,row_f,row_c,cmax=0,idx,noff;
216:   PetscInt          *lcid,*gcid,*lsparse,*gsparse,*colmap,*pcols;
217:   const PetscInt    *rcol;
218:   PetscReal         *Amax_pos,*Amax_neg;
219:   PetscScalar       g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta,pij;
220:   PetscScalar       *pvals;
221:   const PetscScalar *rval;
222:   Mat               lA,gA=NULL;
223:   MatType           mtype;
224:   Vec               C,lvec;
225:   PetscLayout       clayout;
226:   PetscSF           sf;
227:   Mat_MPIAIJ        *mpiaij;

230:   MatGetOwnershipRange(A,&fs,&fe);
231:   fn = fe-fs;
232:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
233:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSEQAIJ);
234:   if (!isMPIAIJ && !isSEQAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Classical AMG requires MPIAIJ matrix");
235:   if (isMPIAIJ) {
236:     mpiaij = (Mat_MPIAIJ*)A->data;
237:     lA = mpiaij->A;
238:     gA = mpiaij->B;
239:     lvec = mpiaij->lvec;
240:     VecGetSize(lvec,&noff);
241:     colmap = mpiaij->garray;
242:     MatGetLayouts(A,NULL,&clayout);
243:     PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
244:     PetscSFSetGraphLayout(sf,clayout,noff,NULL,PETSC_COPY_VALUES,colmap);
245:     PetscMalloc1(noff,&gcid);
246:   } else {
247:     lA = A;
248:   }
249:   PetscMalloc5(fn,&lsparse,fn,&gsparse,fn,&lcid,fn,&Amax_pos,fn,&Amax_neg);

251:   /* count the number of coarse unknowns */
252:   cn = 0;
253:   for (i=0;i<fn;i++) {
254:     /* filter out singletons */
255:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
256:     lcid[i] = -1;
257:     if (!iscoarse) {
258:       cn++;
259:     }
260:   }

262:    /* create the coarse vector */
263:   VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&C);
264:   VecGetOwnershipRange(C,&cs,&ce);

266:   cn = 0;
267:   for (i=0;i<fn;i++) {
268:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
269:     if (!iscoarse) {
270:       lcid[i] = cs+cn;
271:       cn++;
272:     } else {
273:       lcid[i] = -1;
274:     }
275:   }

277:   if (gA) {
278:     PetscSFBcastBegin(sf,MPIU_INT,lcid,gcid);
279:     PetscSFBcastEnd(sf,MPIU_INT,lcid,gcid);
280:   }

282:   /* determine the biggest off-diagonal entries in each row */
283:   for (i=fs;i<fe;i++) {
284:     Amax_pos[i-fs] = 0.;
285:     Amax_neg[i-fs] = 0.;
286:     MatGetRow(A,i,&ncols,&rcol,&rval);
287:     for(j=0;j<ncols;j++){
288:       if ((PetscRealPart(-rval[j]) > Amax_neg[i-fs]) && i != rcol[j]) Amax_neg[i-fs] = PetscAbsScalar(rval[j]);
289:       if ((PetscRealPart(rval[j])  > Amax_pos[i-fs]) && i != rcol[j]) Amax_pos[i-fs] = PetscAbsScalar(rval[j]);
290:     }
291:     if (ncols > cmax) cmax = ncols;
292:     MatRestoreRow(A,i,&ncols,&rcol,&rval);
293:   }
294:   PetscMalloc2(cmax,&pcols,cmax,&pvals);
295:   VecDestroy(&C);

297:   /* count the on and off processor sparsity patterns for the prolongator */
298:   for (i=0;i<fn;i++) {
299:     /* on */
300:     lsparse[i] = 0;
301:     gsparse[i] = 0;
302:     if (lcid[i] >= 0) {
303:       lsparse[i] = 1;
304:       gsparse[i] = 0;
305:     } else {
306:       MatGetRow(lA,i,&ncols,&rcol,&rval);
307:       for (j = 0;j < ncols;j++) {
308:         col = rcol[j];
309:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
310:           lsparse[i] += 1;
311:         }
312:       }
313:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);
314:       /* off */
315:       if (gA) {
316:         MatGetRow(gA,i,&ncols,&rcol,&rval);
317:         for (j = 0; j < ncols; j++) {
318:           col = rcol[j];
319:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
320:             gsparse[i] += 1;
321:           }
322:         }
323:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
324:       }
325:     }
326:   }

328:   /* preallocate and create the prolongator */
329:   MatCreate(PetscObjectComm((PetscObject)A),P);
330:   MatGetType(G,&mtype);
331:   MatSetType(*P,mtype);
332:   MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);
333:   MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);
334:   MatSeqAIJSetPreallocation(*P,0,lsparse);

336:   /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
337:   for (i = 0;i < fn;i++) {
338:     /* determine on or off */
339:     row_f = i + fs;
340:     row_c = lcid[i];
341:     if (row_c >= 0) {
342:       pij = 1.;
343:       MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);
344:     } else {
345:       g_pos = 0.;
346:       g_neg = 0.;
347:       a_pos = 0.;
348:       a_neg = 0.;
349:       diag = 0.;

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

374:       /* ghosted connections */
375:       if (gA) {
376:         MatGetRow(gA,i,&ncols,&rcol,&rval);
377:         for (j = 0; j < ncols; j++) {
378:           col = rcol[j];
379:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
380:             if (PetscRealPart(rval[j]) > 0.) {
381:               g_pos += rval[j];
382:             } else {
383:               g_neg += rval[j];
384:             }
385:           }
386:           if (PetscRealPart(rval[j]) > 0.) {
387:             a_pos += rval[j];
388:           } else {
389:             a_neg += rval[j];
390:           }
391:         }
392:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
393:       }

395:       if (g_neg == 0.) {
396:         alpha = 0.;
397:       } else {
398:         alpha = -a_neg/g_neg;
399:       }

401:       if (g_pos == 0.) {
402:         diag += a_pos;
403:         beta = 0.;
404:       } else {
405:         beta = -a_pos/g_pos;
406:       }
407:       if (diag == 0.) {
408:         invdiag = 0.;
409:       } else invdiag = 1. / diag;
410:       /* on */
411:       MatGetRow(lA,i,&ncols,&rcol,&rval);
412:       idx = 0;
413:       for (j = 0;j < ncols;j++) {
414:         col = rcol[j];
415:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
416:           row_f = i + fs;
417:           row_c = lcid[col];
418:           /* set the values for on-processor ones */
419:           if (PetscRealPart(rval[j]) < 0.) {
420:             pij = rval[j]*alpha*invdiag;
421:           } else {
422:             pij = rval[j]*beta*invdiag;
423:           }
424:           if (PetscAbsScalar(pij) != 0.) {
425:             pvals[idx] = pij;
426:             pcols[idx] = row_c;
427:             idx++;
428:           }
429:         }
430:       }
431:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);
432:       /* off */
433:       if (gA) {
434:         MatGetRow(gA,i,&ncols,&rcol,&rval);
435:         for (j = 0; j < ncols; j++) {
436:           col = rcol[j];
437:           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
438:             row_f = i + fs;
439:             row_c = gcid[col];
440:             /* set the values for on-processor ones */
441:             if (PetscRealPart(rval[j]) < 0.) {
442:               pij = rval[j]*alpha*invdiag;
443:             } else {
444:               pij = rval[j]*beta*invdiag;
445:             }
446:             if (PetscAbsScalar(pij) != 0.) {
447:               pvals[idx] = pij;
448:               pcols[idx] = row_c;
449:               idx++;
450:             }
451:           }
452:         }
453:         MatRestoreRow(gA,i,&ncols,&rcol,&rval);
454:       }
455:       MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);
456:     }
457:   }

459:   MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
460:   MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);

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

464:   PetscFree2(pcols,pvals);
465:   if (gA) {
466:     PetscSFDestroy(&sf);
467:     PetscFree(gcid);
468:   }

470:   return(0);
471: }

475: PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc,Mat *P)
476: {
477:   PetscInt          j,i,ps,pf,pn,pcs,pcf,pcn,idx,cmax;
478:   PetscErrorCode    ierr;
479:   const PetscScalar *pval;
480:   const PetscInt    *pcol;
481:   PetscScalar       *pnval;
482:   PetscInt          *pncol;
483:   PetscInt          ncols;
484:   Mat               Pnew;
485:   PetscInt          *lsparse,*gsparse;
486:   PetscReal         pmax_pos,pmax_neg,ptot_pos,ptot_neg,pthresh_pos,pthresh_neg;
487:   PC_MG             *mg          = (PC_MG*)pc->data;
488:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
489:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
490:   MatType           mtype;

493:   /* trim and rescale with reallocation */
494:   MatGetOwnershipRange(*P,&ps,&pf);
495:   MatGetOwnershipRangeColumn(*P,&pcs,&pcf);
496:   pn = pf-ps;
497:   pcn = pcf-pcs;
498:   PetscMalloc2(pn,&lsparse,pn,&gsparse);
499:   /* allocate */
500:   cmax = 0;
501:   for (i=ps;i<pf;i++) {
502:     lsparse[i-ps] = 0;
503:     gsparse[i-ps] = 0;
504:     MatGetRow(*P,i,&ncols,&pcol,&pval);
505:     if (ncols > cmax) {
506:       cmax = ncols;
507:     }
508:     pmax_pos = 0.;
509:     pmax_neg = 0.;
510:     for (j=0;j<ncols;j++) {
511:       if (PetscRealPart(pval[j]) > pmax_pos) {
512:         pmax_pos = PetscRealPart(pval[j]);
513:       } else if (PetscRealPart(pval[j]) < pmax_neg) {
514:         pmax_neg = PetscRealPart(pval[j]);
515:       }
516:     }
517:     for (j=0;j<ncols;j++) {
518:       if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
519:         if (pcol[j] >= pcs && pcol[j] < pcf) {
520:           lsparse[i-ps]++;
521:         } else {
522:           gsparse[i-ps]++;
523:         }
524:       }
525:     }
526:     MatRestoreRow(*P,i,&ncols,&pcol,&pval);
527:   }

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

531:   MatGetType(*P,&mtype);
532:   MatCreate(PetscObjectComm((PetscObject)*P),&Pnew);
533:   MatSetType(Pnew, mtype);
534:   MatSetSizes(Pnew,pn,pcn,PETSC_DETERMINE,PETSC_DETERMINE);
535:   MatSeqAIJSetPreallocation(Pnew,0,lsparse);
536:   MatMPIAIJSetPreallocation(Pnew,0,lsparse,0,gsparse);

538:   for (i=ps;i<pf;i++) {
539:     MatGetRow(*P,i,&ncols,&pcol,&pval);
540:     pmax_pos = 0.;
541:     pmax_neg = 0.;
542:     for (j=0;j<ncols;j++) {
543:       if (PetscRealPart(pval[j]) > pmax_pos) {
544:         pmax_pos = PetscRealPart(pval[j]);
545:       } else if (PetscRealPart(pval[j]) < pmax_neg) {
546:         pmax_neg = PetscRealPart(pval[j]);
547:       }
548:     }
549:     pthresh_pos = 0.;
550:     pthresh_neg = 0.;
551:     ptot_pos = 0.;
552:     ptot_neg = 0.;
553:     for (j=0;j<ncols;j++) {
554:       if (PetscRealPart(pval[j]) >= cls->interp_threshold*pmax_pos) {
555:         pthresh_pos += PetscRealPart(pval[j]);
556:       } else if (PetscRealPart(pval[j]) <= cls->interp_threshold*pmax_neg) {
557:         pthresh_neg += PetscRealPart(pval[j]);
558:       }
559:       if (PetscRealPart(pval[j]) > 0.) {
560:         ptot_pos += PetscRealPart(pval[j]);
561:       } else {
562:         ptot_neg += PetscRealPart(pval[j]);
563:       }
564:     }
565:     if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
566:     if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
567:     idx=0;
568:     for (j=0;j<ncols;j++) {
569:       if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold) {
570:         pnval[idx] = ptot_pos*pval[j];
571:         pncol[idx] = pcol[j];
572:         idx++;
573:       } else if (PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
574:         pnval[idx] = ptot_neg*pval[j];
575:         pncol[idx] = pcol[j];
576:         idx++;
577:       }
578:     }
579:     MatRestoreRow(*P,i,&ncols,&pcol,&pval);
580:     MatSetValues(Pnew,1,&i,idx,pncol,pnval,INSERT_VALUES);
581:   }

583:   MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY);
584:   MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY);
585:   MatDestroy(P);

587:   *P = Pnew;
588:   PetscFree2(lsparse,gsparse);
589:   PetscFree2(pnval,pncol);
590:   return(0);
591: }

595: PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
596: {
597:   PetscErrorCode    ierr;
598:   Mat               lA,*lAs;
599:   MatType           mtype;
600:   Vec               cv;
601:   PetscInt          *gcid,*lcid,*lsparse,*gsparse,*picol;
602:   PetscInt          fs,fe,cs,ce,nl,i,j,k,li,lni,ci,ncols,maxcols,fn,cn,cid;
603:   PetscMPIInt       size;
604:   const PetscInt    *lidx,*icol,*gidx;
605:   PetscBool         iscoarse;
606:   PetscScalar       vi,pentry,pjentry;
607:   PetscScalar       *pcontrib,*pvcol;
608:   const PetscScalar *vcol;
609:   PetscReal         diag,jdiag,jwttotal;
610:   PetscInt          pncols;
611:   PetscSF           sf;
612:   PetscLayout       clayout;
613:   IS                lis;

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

851: PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc,Mat A,Mat *P)
852: {

854:   PetscErrorCode    ierr;
855:   PetscInt          f,s,n,cf,cs,i,idx;
856:   PetscInt          *coarserows;
857:   PetscInt          ncols;
858:   const PetscInt    *pcols;
859:   const PetscScalar *pvals;
860:   Mat               Pnew;
861:   Vec               diag;
862:   PC_MG             *mg          = (PC_MG*)pc->data;
863:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
864:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

867:   if (cls->nsmooths == 0) {
868:     PCGAMGTruncateProlongator_Private(pc,P);
869:     return(0);
870:   }
871:   MatGetOwnershipRange(*P,&s,&f);
872:   n = f-s;
873:   MatGetOwnershipRangeColumn(*P,&cs,&cf);
874:   PetscMalloc1(n,&coarserows);
875:   /* identify the rows corresponding to coarse unknowns */
876:   idx = 0;
877:   for (i=s;i<f;i++) {
878:     MatGetRow(*P,i,&ncols,&pcols,&pvals);
879:     /* assume, for now, that it's a coarse unknown if it has a single unit entry */
880:     if (ncols == 1) {
881:       if (pvals[0] == 1.) {
882:         coarserows[idx] = i;
883:         idx++;
884:       }
885:     }
886:     MatRestoreRow(*P,i,&ncols,&pcols,&pvals);
887:   }
888:   MatCreateVecs(A,&diag,0);
889:   MatGetDiagonal(A,diag);
890:   VecReciprocal(diag);
891:   for (i=0;i<cls->nsmooths;i++) {
892:     MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew);
893:     MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL);
894:     MatDiagonalScale(Pnew,diag,0);
895:     MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN);
896:     MatDestroy(P);
897:     *P  = Pnew;
898:     Pnew = NULL;
899:   }
900:   VecDestroy(&diag);
901:   PetscFree(coarserows);
902:   PCGAMGTruncateProlongator_Private(pc,P);
903:   return(0);
904: }

908: PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
909: {
910:   PetscErrorCode    ierr;
911:   PetscErrorCode    (*f)(PC,Mat,Mat,PetscCoarsenData*,Mat*);
912:   PC_MG             *mg          = (PC_MG*)pc->data;
913:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
914:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

917:   PetscFunctionListFind(PCGAMGClassicalProlongatorList,cls->prolongtype,&f);
918:   if (!f)SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot find PCGAMG Classical prolongator type");
919:   (*f)(pc,A,G,agg_lists,P);
920:   return(0);
921: }

925: PetscErrorCode PCGAMGDestroy_Classical(PC pc)
926: {
928:   PC_MG          *mg          = (PC_MG*)pc->data;
929:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;

932:   PetscFree(pc_gamg->subctx);
933:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",NULL);
934:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",NULL);
935:   return(0);
936: }

940: PetscErrorCode PCGAMGSetFromOptions_Classical(PetscOptionItems *PetscOptionsObject,PC pc)
941: {
942:   PC_MG             *mg          = (PC_MG*)pc->data;
943:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
944:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
945:   char              tname[256];
946:   PetscErrorCode    ierr;
947:   PetscBool         flg;

950:   PetscOptionsHead(PetscOptionsObject,"GAMG-Classical options");
951:   PetscOptionsFList("-pc_gamg_classical_type","Type of Classical AMG prolongation","PCGAMGClassicalSetType",PCGAMGClassicalProlongatorList,cls->prolongtype, tname, sizeof(tname), &flg);
952:   if (flg) {
953:     PCGAMGClassicalSetType(pc,tname);
954:   }
955:   PetscOptionsReal("-pc_gamg_classical_interp_threshold","Threshold for classical interpolator entries","",cls->interp_threshold,&cls->interp_threshold,NULL);
956:   PetscOptionsInt("-pc_gamg_classical_nsmooths","Threshold for classical interpolator entries","",cls->nsmooths,&cls->nsmooths,NULL);
957:   PetscOptionsTail();
958:   return(0);
959: }

963: PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
964: {
965:   PC_MG          *mg      = (PC_MG*)pc->data;
966:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;

969:   /* no data for classical AMG */
970:   pc_gamg->data = NULL;
971:   pc_gamg->data_cell_cols = 0;
972:   pc_gamg->data_cell_rows = 0;
973:   pc_gamg->data_sz        = 0;
974:   return(0);
975: }


980: PetscErrorCode PCGAMGClassicalFinalizePackage(void)
981: {

985:   PCGAMGClassicalPackageInitialized = PETSC_FALSE;
986:   PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList);
987:   return(0);
988: }

992: PetscErrorCode PCGAMGClassicalInitializePackage(void)
993: {

997:   if (PCGAMGClassicalPackageInitialized) return(0);
998:   PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALDIRECT,PCGAMGProlongator_Classical_Direct);
999:   PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALSTANDARD,PCGAMGProlongator_Classical_Standard);
1000:   PetscRegisterFinalize(PCGAMGClassicalFinalizePackage);
1001:   return(0);
1002: }

1004: /* -------------------------------------------------------------------------- */
1005: /*
1006:    PCCreateGAMG_Classical

1008: */
1011: PetscErrorCode  PCCreateGAMG_Classical(PC pc)
1012: {
1014:   PC_MG             *mg      = (PC_MG*)pc->data;
1015:   PC_GAMG           *pc_gamg = (PC_GAMG*)mg->innerctx;
1016:   PC_GAMG_Classical *pc_gamg_classical;

1019:   PCGAMGClassicalInitializePackage();
1020:   if (pc_gamg->subctx) {
1021:     /* call base class */
1022:     PCDestroy_GAMG(pc);
1023:   }

1025:   /* create sub context for SA */
1026:   PetscNewLog(pc,&pc_gamg_classical);
1027:   pc_gamg->subctx = pc_gamg_classical;
1028:   pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
1029:   /* reset does not do anything; setup not virtual */

1031:   /* set internal function pointers */
1032:   pc_gamg->ops->destroy        = PCGAMGDestroy_Classical;
1033:   pc_gamg->ops->graph          = PCGAMGGraph_Classical;
1034:   pc_gamg->ops->coarsen        = PCGAMGCoarsen_Classical;
1035:   pc_gamg->ops->prolongator    = PCGAMGProlongator_Classical;
1036:   pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
1037:   pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;

1039:   pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
1040:   pc_gamg_classical->interp_threshold = 0.2;
1041:   pc_gamg_classical->nsmooths         = 0;
1042:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",PCGAMGClassicalSetType_GAMG);
1043:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",PCGAMGClassicalGetType_GAMG);
1044:   PCGAMGClassicalSetType(pc,PCGAMGCLASSICALSTANDARD);
1045:   return(0);
1046: }