Actual source code: classical.c

petsc-3.6.1 2015-08-06
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;



196:   /* construct the graph if necessary */
197:   if (!G) {
198:     SETERRQ(fcomm,PETSC_ERR_ARG_WRONGSTATE,"Must set Graph in PC in PCGAMG before coarsening");
199:   }

201:   MatCoarsenCreate(fcomm,&crs);
202:   MatCoarsenSetFromOptions(crs);
203:   MatCoarsenSetAdjacency(crs,*G);
204:   MatCoarsenSetStrictAggs(crs,PETSC_TRUE);
205:   MatCoarsenApply(crs);
206:   MatCoarsenGetData(crs,agg_lists);
207:   MatCoarsenDestroy(&crs);

209:   return(0);
210: }

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

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

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

267:    /* create the coarse vector */
268:   VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&C);
269:   VecGetOwnershipRange(C,&cs,&ce);

271:   cn = 0;
272:   for (i=0;i<fn;i++) {
273:     PetscCDEmptyAt(agg_lists,i,&iscoarse);
274:     if (!iscoarse) {
275:       lcid[i] = cs+cn;
276:       cn++;
277:     } else {
278:       lcid[i] = -1;
279:     }
280:   }

282:   if (gA) {
283:     PetscSFBcastBegin(sf,MPIU_INT,lcid,gcid);
284:     PetscSFBcastEnd(sf,MPIU_INT,lcid,gcid);
285:   }

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

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

333:   /* preallocate and create the prolongator */
334:   MatCreate(PetscObjectComm((PetscObject)A),P);
335:   MatGetType(G,&mtype);
336:   MatSetType(*P,mtype);
337:   MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);
338:   MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);
339:   MatSeqAIJSetPreallocation(*P,0,lsparse);

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

356:       /* local connections */
357:       MatGetRow(lA,i,&ncols,&rcol,&rval);
358:       for (j = 0; j < ncols; j++) {
359:         col = rcol[j];
360:         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*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 (col != i) {
368:           if (PetscRealPart(rval[j]) > 0.) {
369:             a_pos += rval[j];
370:           } else {
371:             a_neg += rval[j];
372:           }
373:         } else {
374:           diag = rval[j];
375:         }
376:       }
377:       MatRestoreRow(lA,i,&ncols,&rcol,&rval);

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

400:       if (g_neg == 0.) {
401:         alpha = 0.;
402:       } else {
403:         alpha = -a_neg/g_neg;
404:       }

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

464:   MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
465:   MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);

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

469:   PetscFree2(pcols,pvals);
470:   if (gA) {
471:     PetscSFDestroy(&sf);
472:     PetscFree(gcid);
473:   }

475:   return(0);
476: }

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

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

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

536:   MatGetType(*P,&mtype);
537:   MatCreate(PetscObjectComm((PetscObject)*P),&Pnew);
538:   MatSetType(Pnew, mtype);
539:   MatSetSizes(Pnew,pn,pcn,PETSC_DETERMINE,PETSC_DETERMINE);
540:   MatSeqAIJSetPreallocation(Pnew,0,lsparse);
541:   MatMPIAIJSetPreallocation(Pnew,0,lsparse,0,gsparse);

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

588:   MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY);
589:   MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY);
590:   MatDestroy(P);

592:   *P = Pnew;
593:   PetscFree2(lsparse,gsparse);
594:   PetscFree2(pnval,pncol);
595:   return(0);
596: }

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

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

858: PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc,Mat A,Mat *P)
859: {

861:   PetscErrorCode    ierr;
862:   PetscInt          f,s,n,cf,cs,i,idx;
863:   PetscInt          *coarserows;
864:   PetscInt          ncols;
865:   const PetscInt    *pcols;
866:   const PetscScalar *pvals;
867:   Mat               Pnew;
868:   Vec               diag;
869:   PC_MG             *mg          = (PC_MG*)pc->data;
870:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
871:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

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

915: PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
916: {
917:   PetscErrorCode    ierr;
918:   PetscErrorCode    (*f)(PC,Mat,Mat,PetscCoarsenData*,Mat*);
919:   PC_MG             *mg          = (PC_MG*)pc->data;
920:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
921:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;

924:   PetscFunctionListFind(PCGAMGClassicalProlongatorList,cls->prolongtype,&f);
925:   if (!f)SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot find PCGAMG Classical prolongator type");
926:   (*f)(pc,A,G,agg_lists,P);
927:   return(0);
928: }

932: PetscErrorCode PCGAMGDestroy_Classical(PC pc)
933: {
935:   PC_MG          *mg          = (PC_MG*)pc->data;
936:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;

939:   PetscFree(pc_gamg->subctx);
940:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",NULL);
941:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",NULL);
942:   return(0);
943: }

947: PetscErrorCode PCGAMGSetFromOptions_Classical(PetscOptions *PetscOptionsObject,PC pc)
948: {
949:   PC_MG             *mg          = (PC_MG*)pc->data;
950:   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
951:   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
952:   char              tname[256];
953:   PetscErrorCode    ierr;
954:   PetscBool         flg;

957:   PetscOptionsHead(PetscOptionsObject,"GAMG-Classical options");
958:   PetscOptionsFList("-pc_gamg_classical_type","Type of Classical AMG prolongation","PCGAMGClassicalSetType",PCGAMGClassicalProlongatorList,cls->prolongtype, tname, sizeof(tname), &flg);
959:   if (flg) {
960:     PCGAMGClassicalSetType(pc,tname);
961:   }
962:   PetscOptionsReal("-pc_gamg_classical_interp_threshold","Threshold for classical interpolator entries","",cls->interp_threshold,&cls->interp_threshold,NULL);
963:   PetscOptionsInt("-pc_gamg_classical_nsmooths","Threshold for classical interpolator entries","",cls->nsmooths,&cls->nsmooths,NULL);
964:   PetscOptionsTail();
965:   return(0);
966: }

970: PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
971: {
972:   PC_MG          *mg      = (PC_MG*)pc->data;
973:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;

976:   /* no data for classical AMG */
977:   pc_gamg->data = NULL;
978:   pc_gamg->data_cell_cols = 0;
979:   pc_gamg->data_cell_rows = 0;
980:   pc_gamg->data_sz        = 0;
981:   return(0);
982: }


987: PetscErrorCode PCGAMGClassicalFinalizePackage(void)
988: {

992:   PCGAMGClassicalPackageInitialized = PETSC_FALSE;
993:   PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList);
994:   return(0);
995: }

999: PetscErrorCode PCGAMGClassicalInitializePackage(void)
1000: {

1004:   if (PCGAMGClassicalPackageInitialized) return(0);
1005:   PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALDIRECT,PCGAMGProlongator_Classical_Direct);
1006:   PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALSTANDARD,PCGAMGProlongator_Classical_Standard);
1007:   PetscRegisterFinalize(PCGAMGClassicalFinalizePackage);
1008:   return(0);
1009: }

1011: /* -------------------------------------------------------------------------- */
1012: /*
1013:    PCCreateGAMG_Classical

1015: */
1018: PetscErrorCode  PCCreateGAMG_Classical(PC pc)
1019: {
1021:   PC_MG             *mg      = (PC_MG*)pc->data;
1022:   PC_GAMG           *pc_gamg = (PC_GAMG*)mg->innerctx;
1023:   PC_GAMG_Classical *pc_gamg_classical;

1026:   PCGAMGClassicalInitializePackage();
1027:   if (pc_gamg->subctx) {
1028:     /* call base class */
1029:     PCDestroy_GAMG(pc);
1030:   }

1032:   /* create sub context for SA */
1033:   PetscNewLog(pc,&pc_gamg_classical);
1034:   pc_gamg->subctx = pc_gamg_classical;
1035:   pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
1036:   /* reset does not do anything; setup not virtual */

1038:   /* set internal function pointers */
1039:   pc_gamg->ops->destroy        = PCGAMGDestroy_Classical;
1040:   pc_gamg->ops->graph          = PCGAMGGraph_Classical;
1041:   pc_gamg->ops->coarsen        = PCGAMGCoarsen_Classical;
1042:   pc_gamg->ops->prolongator    = PCGAMGProlongator_Classical;
1043:   pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
1044:   pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;

1046:   pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
1047:   pc_gamg_classical->interp_threshold = 0.2;
1048:   pc_gamg_classical->nsmooths         = 0;
1049:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",PCGAMGClassicalSetType_GAMG);
1050:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",PCGAMGClassicalGetType_GAMG);
1051:   PCGAMGClassicalSetType(pc,PCGAMGCLASSICALSTANDARD);
1052:   return(0);
1053: }