Actual source code: agg.c

petsc-3.7.3 2016-07-24
Report Typos and Errors
  1: /*
  2:  GAMG geometric-algebric multiogrid PC - Mark Adams 2011
  3:  */

  5: #include <../src/ksp/pc/impls/gamg/gamg.h>        /*I "petscpc.h" I*/
  6: #include <petsc/private/kspimpl.h>

  8: #include <petscblaslapack.h>

 10: typedef struct {
 11:   PetscInt  nsmooths;
 12:   PetscBool sym_graph;
 13:   PetscInt square_graph;
 14: } PC_GAMG_AGG;

 18: /*@
 19:    PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical)

 21:    Not Collective on PC

 23:    Input Parameters:
 24: .  pc - the preconditioner context

 26:    Options Database Key:
 27: .  -pc_gamg_agg_nsmooths

 29:    Level: intermediate

 31:    Concepts: Aggregation AMG preconditioner

 33: .seealso: ()
 34: @*/
 35: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
 36: {

 41:   PetscTryMethod(pc,"PCGAMGSetNSmooths_C",(PC,PetscInt),(pc,n));
 42:   return(0);
 43: }

 47: static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n)
 48: {
 49:   PC_MG       *mg          = (PC_MG*)pc->data;
 50:   PC_GAMG     *pc_gamg     = (PC_GAMG*)mg->innerctx;
 51:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

 54:   pc_gamg_agg->nsmooths = n;
 55:   return(0);
 56: }

 60: /*@
 61:    PCGAMGSetSymGraph -

 63:    Not Collective on PC

 65:    Input Parameters:
 66: .  pc - the preconditioner context

 68:    Options Database Key:
 69: .  -pc_gamg_sym_graph

 71:    Level: intermediate

 73:    Concepts: Aggregation AMG preconditioner

 75: .seealso: ()
 76: @*/
 77: PetscErrorCode PCGAMGSetSymGraph(PC pc, PetscBool n)
 78: {

 83:   PetscTryMethod(pc,"PCGAMGSetSymGraph_C",(PC,PetscBool),(pc,n));
 84:   return(0);
 85: }

 89: static PetscErrorCode PCGAMGSetSymGraph_AGG(PC pc, PetscBool n)
 90: {
 91:   PC_MG       *mg          = (PC_MG*)pc->data;
 92:   PC_GAMG     *pc_gamg     = (PC_GAMG*)mg->innerctx;
 93:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

 96:   pc_gamg_agg->sym_graph = n;
 97:   return(0);
 98: }

102: /*@
103:    PCGAMGSetSquareGraph -

105:    Not Collective on PC

107:    Input Parameters:
108: .  pc - the preconditioner context

110:    Options Database Key:
111: .  -pc_gamg_square_graph

113:    Level: intermediate

115:    Concepts: Aggregation AMG preconditioner

117: .seealso: ()
118: @*/
119: PetscErrorCode PCGAMGSetSquareGraph(PC pc, PetscInt n)
120: {

125:   PetscTryMethod(pc,"PCGAMGSetSquareGraph_C",(PC,PetscInt),(pc,n));
126:   return(0);
127: }

131: static PetscErrorCode PCGAMGSetSquareGraph_AGG(PC pc, PetscInt n)
132: {
133:   PC_MG       *mg          = (PC_MG*)pc->data;
134:   PC_GAMG     *pc_gamg     = (PC_GAMG*)mg->innerctx;
135:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

138:   pc_gamg_agg->square_graph = n;
139:   return(0);
140: }

142: /* -------------------------------------------------------------------------- */
143: /*
144:    PCSetFromOptions_GAMG_AGG

146:   Input Parameter:
147:    . pc -
148: */
151: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PetscOptionItems *PetscOptionsObject,PC pc)
152: {
154:   PC_MG          *mg          = (PC_MG*)pc->data;
155:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
156:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

159:   PetscOptionsHead(PetscOptionsObject,"GAMG-AGG options");
160:   {
161:     /* -pc_gamg_agg_nsmooths */
162:     pc_gamg_agg->nsmooths = 1;

164:     PetscOptionsInt("-pc_gamg_agg_nsmooths","smoothing steps for smoothed aggregation, usually 1","PCGAMGSetNSmooths",pc_gamg_agg->nsmooths,&pc_gamg_agg->nsmooths,NULL);

166:     /* -pc_gamg_sym_graph */
167:     pc_gamg_agg->sym_graph = PETSC_FALSE;

169:     PetscOptionsBool("-pc_gamg_sym_graph","Set for asymmetric matrices","PCGAMGSetSymGraph",pc_gamg_agg->sym_graph,&pc_gamg_agg->sym_graph,NULL);

171:     /* -pc_gamg_square_graph */
172:     pc_gamg_agg->square_graph = 1;

174:     PetscOptionsInt("-pc_gamg_square_graph","Number of levels to square graph for faster coarsening and lower coarse grid complexity","PCGAMGSetSquareGraph",pc_gamg_agg->square_graph,&pc_gamg_agg->square_graph,NULL);
175:   }
176:   PetscOptionsTail();
177:   return(0);
178: }

180: /* -------------------------------------------------------------------------- */
181: /*
182:    PCDestroy_AGG

184:   Input Parameter:
185:    . pc -
186: */
189: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
190: {
192:   PC_MG          *mg          = (PC_MG*)pc->data;
193:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;

196:   PetscFree(pc_gamg->subctx);
197:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);
198:   return(0);
199: }

201: /* -------------------------------------------------------------------------- */
202: /*
203:    PCSetCoordinates_AGG
204:      - collective

206:    Input Parameter:
207:    . pc - the preconditioner context
208:    . ndm - dimesion of data (used for dof/vertex for Stokes)
209:    . a_nloc - number of vertices local
210:    . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
211: */

215: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
216: {
217:   PC_MG          *mg      = (PC_MG*)pc->data;
218:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
220:   PetscInt       arrsz,kk,ii,jj,nloc,ndatarows,ndf;
221:   Mat            mat = pc->pmat;

226:   nloc = a_nloc;

228:   /* SA: null space vectors */
229:   MatGetBlockSize(mat, &ndf); /* this does not work for Stokes */
230:   if (coords && ndf==1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
231:   else if (coords) {
232:     if (ndm > ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"degrees of motion %d > block size %d",ndm,ndf);
233:     pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* displacement elasticity */
234:     if (ndm != ndf) {
235:       if (pc_gamg->data_cell_cols != ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Don't know how to create null space for ndm=%d, ndf=%d.  Use MatSetNearNullSpace.",ndm,ndf);
236:     }
237:   } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
238:   pc_gamg->data_cell_rows = ndatarows = ndf;
239:   if (pc_gamg->data_cell_cols <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->data_cell_cols %D <= 0",pc_gamg->data_cell_cols);
240:   arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols;

242:   /* create data - syntactic sugar that should be refactored at some point */
243:   if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
244:     PetscFree(pc_gamg->data);
245:     PetscMalloc1(arrsz+1, &pc_gamg->data);
246:   }
247:   /* copy data in - column oriented */
248:   for (kk=0; kk<nloc; kk++) {
249:     const PetscInt M     = nloc*pc_gamg->data_cell_rows; /* stride into data */
250:     PetscReal      *data = &pc_gamg->data[kk*ndatarows]; /* start of cell */
251:     if (pc_gamg->data_cell_cols==1) *data = 1.0;
252:     else {
253:       /* translational modes */
254:       for (ii=0;ii<ndatarows;ii++) {
255:         for (jj=0;jj<ndatarows;jj++) {
256:           if (ii==jj)data[ii*M + jj] = 1.0;
257:           else data[ii*M + jj] = 0.0;
258:         }
259:       }

261:       /* rotational modes */
262:       if (coords) {
263:         if (ndm == 2) {
264:           data   += 2*M;
265:           data[0] = -coords[2*kk+1];
266:           data[1] =  coords[2*kk];
267:         } else {
268:           data   += 3*M;
269:           data[0] = 0.0;             data[M+0] =  coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1];
270:           data[1] = -coords[3*kk+2]; data[M+1] = 0.0;             data[2*M+1] =  coords[3*kk];
271:           data[2] =  coords[3*kk+1]; data[M+2] = -coords[3*kk];   data[2*M+2] = 0.0;
272:         }
273:       }
274:     }
275:   }

277:   pc_gamg->data_sz = arrsz;
278:   return(0);
279: }

281: typedef PetscInt NState;
282: static const NState NOT_DONE=-2;
283: static const NState DELETED =-1;
284: static const NState REMOVED =-3;
285: #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED)

287: /* -------------------------------------------------------------------------- */
288: /*
289:    smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific
290:      - AGG-MG specific: clears singletons out of 'selected_2'

292:    Input Parameter:
293:    . Gmat_2 - glabal matrix of graph (data not defined)   base (squared) graph
294:    . Gmat_1 - base graph to grab with                 base graph
295:    Input/Output Parameter:
296:    . aggs_2 - linked list of aggs with gids)
297: */
300: static PetscErrorCode smoothAggs(Mat Gmat_2, Mat Gmat_1,PetscCoarsenData *aggs_2)
301: {
303:   PetscBool      isMPI;
304:   Mat_SeqAIJ     *matA_1, *matB_1=0, *matA_2, *matB_2=0;
305:   MPI_Comm       comm;
306:   PetscInt       lid,*ii,*idx,ix,Iend,my0,kk,n,j;
307:   Mat_MPIAIJ     *mpimat_2 = 0, *mpimat_1=0;
308:   const PetscInt nloc      = Gmat_2->rmap->n;
309:   PetscScalar    *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid;
310:   PetscInt       *lid_cprowID_1;
311:   NState         *lid_state;
312:   Vec            ghost_par_orig2;

315:   PetscObjectGetComm((PetscObject)Gmat_2,&comm);
316:   MatGetOwnershipRange(Gmat_1,&my0,&Iend);

318:   if (PETSC_FALSE) {
319:     PetscViewer viewer; char fname[32]; static int llev=0;
320:     sprintf(fname,"Gmat2_%d.m",llev++);
321:     PetscViewerASCIIOpen(comm,fname,&viewer);
322:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
323:     MatView(Gmat_2, viewer);
324:     PetscViewerPopFormat(viewer);
325:     PetscViewerDestroy(&viewer);
326:   }

328:   /* get submatrices */
329:   PetscObjectTypeCompare((PetscObject)Gmat_1, MATMPIAIJ, &isMPI);
330:   if (isMPI) {
331:     /* grab matrix objects */
332:     mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data;
333:     mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data;
334:     matA_1   = (Mat_SeqAIJ*)mpimat_1->A->data;
335:     matB_1   = (Mat_SeqAIJ*)mpimat_1->B->data;
336:     matA_2   = (Mat_SeqAIJ*)mpimat_2->A->data;
337:     matB_2   = (Mat_SeqAIJ*)mpimat_2->B->data;

339:     /* force compressed row storage for B matrix in AuxMat */
340:     MatCheckCompressedRow(mpimat_1->B,matB_1->nonzerorowcnt,&matB_1->compressedrow,matB_1->i,Gmat_1->rmap->n,-1.0);

342:     PetscMalloc1(nloc, &lid_cprowID_1);
343:     for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
344:     for (ix=0; ix<matB_1->compressedrow.nrows; ix++) {
345:       PetscInt lid = matB_1->compressedrow.rindex[ix];
346:       lid_cprowID_1[lid] = ix;
347:     }
348:   } else {
349:     matA_1        = (Mat_SeqAIJ*)Gmat_1->data;
350:     matA_2        = (Mat_SeqAIJ*)Gmat_2->data;
351:     lid_cprowID_1 = NULL;
352:   }
353:   if (nloc>0) {
354:     if (!(matA_1 && !matA_1->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matA_1 && !matA_1->compressedrow.use)");
355:     if (!(!matB_1 || matB_1->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matB_1==0 || matB_1->compressedrow.use)");
356:     if (!(matA_2 && !matA_2->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matA_2 && !matA_2->compressedrow.use)");
357:     if (!(!matB_2 || matB_2->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matB_2==0 || matB_2->compressedrow.use)");
358:   }
359:   /* get state of locals and selected gid for deleted */
360:   PetscMalloc2(nloc, &lid_state,nloc, &lid_parent_gid);
361:   for (lid = 0; lid < nloc; lid++) {
362:     lid_parent_gid[lid] = -1.0;
363:     lid_state[lid]      = DELETED;
364:   }

366:   /* set lid_state */
367:   for (lid = 0; lid < nloc; lid++) {
368:     PetscCDPos pos;
369:     PetscCDGetHeadPos(aggs_2,lid,&pos);
370:     if (pos) {
371:       PetscInt gid1;

373:       PetscLLNGetID(pos, &gid1);
374:       if (gid1 != lid+my0) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"gid1 %D != lid %D + my0 %D",gid1,lid,my0);
375:       lid_state[lid] = gid1;
376:     }
377:   }

379:   /* map local to selected local, DELETED means a ghost owns it */
380:   for (lid=kk=0; lid<nloc; lid++) {
381:     NState state = lid_state[lid];
382:     if (IS_SELECTED(state)) {
383:       PetscCDPos pos;
384:       PetscCDGetHeadPos(aggs_2,lid,&pos);
385:       while (pos) {
386:         PetscInt gid1;
387:         PetscLLNGetID(pos, &gid1);
388:         PetscCDGetNextPos(aggs_2,lid,&pos);

390:         if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1-my0] = (PetscScalar)(lid + my0);
391:       }
392:     }
393:   }
394:   /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
395:   if (isMPI) {
396:     Vec tempVec;
397:     /* get 'cpcol_1_state' */
398:     MatCreateVecs(Gmat_1, &tempVec, 0);
399:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
400:       PetscScalar v = (PetscScalar)lid_state[kk];
401:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
402:     }
403:     VecAssemblyBegin(tempVec);
404:     VecAssemblyEnd(tempVec);
405:     VecScatterBegin(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
406:     VecScatterEnd(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
407:     VecGetArray(mpimat_1->lvec, &cpcol_1_state);
408:     /* get 'cpcol_2_state' */
409:     VecScatterBegin(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
410:       VecScatterEnd(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
411:     VecGetArray(mpimat_2->lvec, &cpcol_2_state);
412:     /* get 'cpcol_2_par_orig' */
413:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
414:       PetscScalar v = (PetscScalar)lid_parent_gid[kk];
415:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
416:     }
417:     VecAssemblyBegin(tempVec);
418:     VecAssemblyEnd(tempVec);
419:     VecDuplicate(mpimat_2->lvec, &ghost_par_orig2);
420:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
421:     VecScatterEnd(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
422:     VecGetArray(ghost_par_orig2, &cpcol_2_par_orig);

424:     VecDestroy(&tempVec);
425:   } /* ismpi */

427:   /* doit */
428:   for (lid=0; lid<nloc; lid++) {
429:     NState state = lid_state[lid];
430:     if (IS_SELECTED(state)) {
431:       /* steal locals */
432:       ii  = matA_1->i; n = ii[lid+1] - ii[lid];
433:       idx = matA_1->j + ii[lid];
434:       for (j=0; j<n; j++) {
435:         PetscInt lidj   = idx[j], sgid;
436:         NState   statej = lid_state[lidj];
437:         if (statej==DELETED && (sgid=(PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid+my0) { /* steal local */
438:           lid_parent_gid[lidj] = (PetscScalar)(lid+my0); /* send this if sgid is not local */
439:           if (sgid >= my0 && sgid < Iend) {       /* I'm stealing this local from a local sgid */
440:             PetscInt   hav=0,slid=sgid-my0,gidj=lidj+my0;
441:             PetscCDPos pos,last=NULL;
442:             /* looking for local from local so id_llist_2 works */
443:             PetscCDGetHeadPos(aggs_2,slid,&pos);
444:             while (pos) {
445:               PetscInt gid;
446:               PetscLLNGetID(pos, &gid);
447:               if (gid == gidj) {
448:                 if (!last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null");
449:                 PetscCDRemoveNextNode(aggs_2, slid, last);
450:                 PetscCDAppendNode(aggs_2, lid, pos);
451:                 hav  = 1;
452:                 break;
453:               } else last = pos;

455:               PetscCDGetNextPos(aggs_2,slid,&pos);
456:             }
457:             if (hav!=1) {
458:               if (!hav) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
459:               SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav);
460:             }
461:           } else {            /* I'm stealing this local, owned by a ghost */
462:             if (sgid != -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Have un-symmetric graph (apparently). Use '-pc_gamg_sym_graph true' to symetrize the graph or '-pc_gamg_threshold -1.0' if the matrix is structurally symmetric.");
463:             PetscCDAppendID(aggs_2, lid, lidj+my0);
464:           }
465:         }
466:       } /* local neighbors */
467:     } else if (state == DELETED && lid_cprowID_1) {
468:       PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
469:       /* see if I have a selected ghost neighbor that will steal me */
470:       if ((ix=lid_cprowID_1[lid]) != -1) {
471:         ii  = matB_1->compressedrow.i; n = ii[ix+1] - ii[ix];
472:         idx = matB_1->j + ii[ix];
473:         for (j=0; j<n; j++) {
474:           PetscInt cpid   = idx[j];
475:           NState   statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
476:           if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */
477:             lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
478:             if (sgidold>=my0 && sgidold<Iend) { /* this was mine */
479:               PetscInt   hav=0,oldslidj=sgidold-my0;
480:               PetscCDPos pos,last=NULL;
481:               /* remove from 'oldslidj' list */
482:               PetscCDGetHeadPos(aggs_2,oldslidj,&pos);
483:               while (pos) {
484:                 PetscInt gid;
485:                 PetscLLNGetID(pos, &gid);
486:                 if (lid+my0 == gid) {
487:                   /* id_llist_2[lastid] = id_llist_2[flid];   /\* remove lid from oldslidj list *\/ */
488:                   if (!last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null");
489:                   PetscCDRemoveNextNode(aggs_2, oldslidj, last);
490:                   /* ghost (PetscScalar)statej will add this later */
491:                   hav = 1;
492:                   break;
493:                 } else last = pos;

495:                 PetscCDGetNextPos(aggs_2,oldslidj,&pos);
496:               }
497:               if (hav!=1) {
498:                 if (!hav) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
499:                 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav);
500:               }
501:             } else {
502:               /* ghosts remove this later */
503:             }
504:           }
505:         }
506:       }
507:     } /* selected/deleted */
508:   } /* node loop */

510:   if (isMPI) {
511:     PetscScalar   *cpcol_2_parent,*cpcol_2_gid;
512:     Vec           tempVec,ghostgids2,ghostparents2;
513:     PetscInt      cpid,nghost_2;
514:     GAMGHashTable gid_cpid;

516:     VecGetSize(mpimat_2->lvec, &nghost_2);
517:     MatCreateVecs(Gmat_2, &tempVec, 0);

519:     /* get 'cpcol_2_parent' */
520:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
521:       VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES);
522:     }
523:     VecAssemblyBegin(tempVec);
524:     VecAssemblyEnd(tempVec);
525:     VecDuplicate(mpimat_2->lvec, &ghostparents2);
526:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
527:     VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
528:     VecGetArray(ghostparents2, &cpcol_2_parent);

530:     /* get 'cpcol_2_gid' */
531:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
532:       PetscScalar v = (PetscScalar)j;
533:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
534:     }
535:     VecAssemblyBegin(tempVec);
536:     VecAssemblyEnd(tempVec);
537:     VecDuplicate(mpimat_2->lvec, &ghostgids2);
538:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
539:     VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
540:     VecGetArray(ghostgids2, &cpcol_2_gid);

542:     VecDestroy(&tempVec);

544:     /* look for deleted ghosts and add to table */
545:     GAMGTableCreate(2*nghost_2+1, &gid_cpid);
546:     for (cpid = 0; cpid < nghost_2; cpid++) {
547:       NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
548:       if (state==DELETED) {
549:         PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
550:         PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
551:         if (sgid_old == -1 && sgid_new != -1) {
552:           PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
553:           GAMGTableAdd(&gid_cpid, gid, cpid);
554:         }
555:       }
556:     }

558:     /* look for deleted ghosts and see if they moved - remove it */
559:     for (lid=0; lid<nloc; lid++) {
560:       NState state = lid_state[lid];
561:       if (IS_SELECTED(state)) {
562:         PetscCDPos pos,last=NULL;
563:         /* look for deleted ghosts and see if they moved */
564:         PetscCDGetHeadPos(aggs_2,lid,&pos);
565:         while (pos) {
566:           PetscInt gid;
567:           PetscLLNGetID(pos, &gid);

569:           if (gid < my0 || gid >= Iend) {
570:             GAMGTableFind(&gid_cpid, gid, &cpid);
571:             if (cpid != -1) {
572:               /* a moved ghost - */
573:               /* id_llist_2[lastid] = id_llist_2[flid];    /\* remove 'flid' from list *\/ */
574:               PetscCDRemoveNextNode(aggs_2, lid, last);
575:             } else last = pos;
576:           } else last = pos;

578:           PetscCDGetNextPos(aggs_2,lid,&pos);
579:         } /* loop over list of deleted */
580:       } /* selected */
581:     }
582:     GAMGTableDestroy(&gid_cpid);

584:     /* look at ghosts, see if they changed - and it */
585:     for (cpid = 0; cpid < nghost_2; cpid++) {
586:       PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
587:       if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
588:         PetscInt   gid     = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
589:         PetscInt   slid_new=sgid_new-my0,hav=0;
590:         PetscCDPos pos;
591:         /* search for this gid to see if I have it */
592:         PetscCDGetHeadPos(aggs_2,slid_new,&pos);
593:         while (pos) {
594:           PetscInt gidj;
595:           PetscLLNGetID(pos, &gidj);
596:           PetscCDGetNextPos(aggs_2,slid_new,&pos);

598:           if (gidj == gid) { hav = 1; break; }
599:         }
600:         if (hav != 1) {
601:           /* insert 'flidj' into head of llist */
602:           PetscCDAppendID(aggs_2, slid_new, gid);
603:         }
604:       }
605:     }

607:     VecRestoreArray(mpimat_1->lvec, &cpcol_1_state);
608:     VecRestoreArray(mpimat_2->lvec, &cpcol_2_state);
609:     VecRestoreArray(ghostparents2, &cpcol_2_parent);
610:     VecRestoreArray(ghostgids2, &cpcol_2_gid);
611:     PetscFree(lid_cprowID_1);
612:     VecDestroy(&ghostgids2);
613:     VecDestroy(&ghostparents2);
614:     VecDestroy(&ghost_par_orig2);
615:   }

617:   PetscFree2(lid_state,lid_parent_gid);
618:   return(0);
619: }

621: /* -------------------------------------------------------------------------- */
622: /*
623:    PCSetData_AGG - called if data is not set with PCSetCoordinates.
624:       Looks in Mat for near null space.
625:       Does not work for Stokes

627:   Input Parameter:
628:    . pc -
629:    . a_A - matrix to get (near) null space out of.
630: */
633: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
634: {
636:   PC_MG          *mg      = (PC_MG*)pc->data;
637:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
638:   MatNullSpace   mnull;

641:   MatGetNearNullSpace(a_A, &mnull);
642:   if (!mnull) {
643:     PetscInt bs,NN,MM;
644:     MatGetBlockSize(a_A, &bs);
645:     MatGetLocalSize(a_A, &MM, &NN);
646:     if (MM % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MM %D must be divisible by bs %D",MM,bs);
647:     PCSetCoordinates_AGG(pc, bs, MM/bs, NULL);
648:   } else {
649:     PetscReal *nullvec;
650:     PetscBool has_const;
651:     PetscInt  i,j,mlocal,nvec,bs;
652:     const Vec *vecs; const PetscScalar *v;
653:     MatGetLocalSize(a_A,&mlocal,NULL);
654:     MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
655:     pc_gamg->data_sz = (nvec+!!has_const)*mlocal;
656:     PetscMalloc1((nvec+!!has_const)*mlocal,&nullvec);
657:     if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0;
658:     for (i=0; i<nvec; i++) {
659:       VecGetArrayRead(vecs[i],&v);
660:       for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]);
661:       VecRestoreArrayRead(vecs[i],&v);
662:     }
663:     pc_gamg->data           = nullvec;
664:     pc_gamg->data_cell_cols = (nvec+!!has_const);

666:     MatGetBlockSize(a_A, &bs);

668:     pc_gamg->data_cell_rows = bs;
669:   }
670:   return(0);
671: }

673: /* -------------------------------------------------------------------------- */
674: /*
675:  formProl0

677:    Input Parameter:
678:    . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices 
679:    . bs - row block size
680:    . nSAvec - column bs of new P
681:    . my0crs - global index of start of locals
682:    . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
683:    . data_in[data_stride*nSAvec] - local data on fine grid
684:    . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
685:   Output Parameter:
686:    . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
687:    . a_Prol - prolongation operator
688: */
691: static PetscErrorCode formProl0(PetscCoarsenData *agg_llists,PetscInt bs,PetscInt nSAvec,PetscInt my0crs,PetscInt data_stride,PetscReal data_in[],const PetscInt flid_fgid[],PetscReal **a_data_out,Mat a_Prol)
692: {
694:   PetscInt       Istart,my0,Iend,nloc,clid,flid = 0,aggID,kk,jj,ii,mm,ndone,nSelected,minsz,nghosts,out_data_stride;
695:   MPI_Comm       comm;
696:   PetscMPIInt    rank;
697:   PetscReal      *out_data;
698:   PetscCDPos     pos;
699:   GAMGHashTable  fgid_flid;

701: /* #define OUT_AGGS */
702: #if defined(OUT_AGGS)
703:   static PetscInt llev = 0; char fname[32]; FILE *file = NULL; PetscInt pM;
704: #endif

707:   PetscObjectGetComm((PetscObject)a_Prol,&comm);
708:   MPI_Comm_rank(comm,&rank);
709:   MatGetOwnershipRange(a_Prol, &Istart, &Iend);
710:   nloc = (Iend-Istart)/bs; my0 = Istart/bs;
711:   if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Iend %D - Istart %d must be divisible by bs %D",Iend,Istart,bs);
712:   Iend   /= bs;
713:   nghosts = data_stride/bs - nloc;

715:   GAMGTableCreate(2*nghosts+1, &fgid_flid);
716:   for (kk=0; kk<nghosts; kk++) {
717:     GAMGTableAdd(&fgid_flid, flid_fgid[nloc+kk], nloc+kk);
718:   }

720: #if defined(OUT_AGGS)
721:   sprintf(fname,"aggs_%d_%d.m",llev++,rank);
722:   if (llev==1) file = fopen(fname,"w");
723:   MatGetSize(a_Prol, &pM, &jj);
724: #endif

726:   /* count selected -- same as number of cols of P */
727:   for (nSelected=mm=0; mm<nloc; mm++) {
728:     PetscBool ise;
729:     PetscCDEmptyAt(agg_llists, mm, &ise);
730:     if (!ise) nSelected++;
731:   }
732:   MatGetOwnershipRangeColumn(a_Prol, &ii, &jj);
733:   if ((ii/nSAvec) != my0crs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D /nSAvec %D  != my0crs %D",ii,nSAvec,my0crs);
734:   if (nSelected != (jj-ii)/nSAvec) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nSelected %D != (jj %D - ii %D)/nSAvec %D",nSelected,jj,ii,nSAvec);

736:   /* aloc space for coarse point data (output) */
737:   out_data_stride = nSelected*nSAvec;

739:   PetscMalloc1(out_data_stride*nSAvec, &out_data);
740:   for (ii=0;ii<out_data_stride*nSAvec;ii++) out_data[ii]=PETSC_MAX_REAL;
741:   *a_data_out = out_data; /* output - stride nSelected*nSAvec */

743:   /* find points and set prolongation */
744:   minsz = 100;
745:   ndone = 0;
746:   for (mm = clid = 0; mm < nloc; mm++) {
747:     PetscCDSizeAt(agg_llists, mm, &jj);
748:     if (jj > 0) {
749:       const PetscInt lid = mm, cgid = my0crs + clid;
750:       PetscInt       cids[100]; /* max bs */
751:       PetscBLASInt   asz  =jj,M=asz*bs,N=nSAvec,INFO;
752:       PetscBLASInt   Mdata=M+((N-M>0) ? N-M : 0),LDA=Mdata,LWORK=N*bs;
753:       PetscScalar    *qqc,*qqr,*TAU,*WORK;
754:       PetscInt       *fids;
755:       PetscReal      *data;
756:       /* count agg */
757:       if (asz<minsz) minsz = asz;

759:       /* get block */
760:       PetscMalloc5(Mdata*N, &qqc,M*N, &qqr,N, &TAU,LWORK, &WORK,M, &fids);

762:       aggID = 0;
763:       PetscCDGetHeadPos(agg_llists,lid,&pos);
764:       while (pos) {
765:         PetscInt gid1;
766:         PetscLLNGetID(pos, &gid1);
767:         PetscCDGetNextPos(agg_llists,lid,&pos);

769:         if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
770:         else {
771:           GAMGTableFind(&fgid_flid, gid1, &flid);
772:           if (flid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot find gid1 in table");
773:         }
774:         /* copy in B_i matrix - column oriented */
775:         data = &data_in[flid*bs];
776:         for (ii = 0; ii < bs; ii++) {
777:           for (jj = 0; jj < N; jj++) {
778:             PetscReal d = data[jj*data_stride + ii];
779:             qqc[jj*Mdata + aggID*bs + ii] = d;
780:           }
781:         }
782: #if defined(OUT_AGGS)
783:         if (llev==1) {
784:           char     str[] = "plot(%e,%e,'r*'), hold on,\n", col[] = "rgbkmc", sim[] = "*os+h>d<vx^";
785:           PetscInt MM,pi,pj;
786:           str[12] = col[(clid+17*rank)%6]; str[13] = sim[(clid+17*rank)%11];
787:           MM      = (PetscInt)(PetscSqrtReal((PetscReal)pM));
788:           pj      = gid1/MM; pi = gid1%MM;
789:           fprintf(file,str,(double)pi,(double)pj);
790:           /* fprintf(file,str,data[2*data_stride+1],-data[2*data_stride]); */
791:         }
792: #endif
793:         /* set fine IDs */
794:         for (kk=0; kk<bs; kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk;

796:         aggID++;
797:       }

799:       /* pad with zeros */
800:       for (ii = asz*bs; ii < Mdata; ii++) {
801:         for (jj = 0; jj < N; jj++, kk++) {
802:           qqc[jj*Mdata + ii] = .0;
803:         }
804:       }

806:       ndone += aggID;
807:       /* QR */
808:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
809:       PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
810:       PetscFPTrapPop();
811:       if (INFO != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xGEQRF error");
812:       /* get R - column oriented - output B_{i+1} */
813:       {
814:         PetscReal *data = &out_data[clid*nSAvec];
815:         for (jj = 0; jj < nSAvec; jj++) {
816:           for (ii = 0; ii < nSAvec; ii++) {
817:             if (data[jj*out_data_stride + ii] != PETSC_MAX_REAL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"data[jj*out_data_stride + ii] != %e",PETSC_MAX_REAL);
818:            if (ii <= jj) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]);
819:            else data[jj*out_data_stride + ii] = 0.;
820:           }
821:         }
822:       }

824:       /* get Q - row oriented */
825:       PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
826:       if (INFO != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xORGQR error arg %d",-INFO);

828:       for (ii = 0; ii < M; ii++) {
829:         for (jj = 0; jj < N; jj++) {
830:           qqr[N*ii + jj] = qqc[jj*Mdata + ii];
831:         }
832:       }

834:       /* add diagonal block of P0 */
835:       for (kk=0; kk<N; kk++) {
836:         cids[kk] = N*cgid + kk; /* global col IDs in P0 */
837:       }
838:       MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);

840:       PetscFree5(qqc,qqr,TAU,WORK,fids);
841:       clid++;
842:     } /* coarse agg */
843:   } /* for all fine nodes */
844:   MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
845:   MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);

847: /* MPIU_Allreduce(&ndone, &ii, 1, MPIU_INT, MPIU_SUM, comm); */
848: /* MatGetSize(a_Prol, &kk, &jj); */
849: /* MPIU_Allreduce(&minsz, &jj, 1, MPIU_INT, MPIU_MIN, comm); */
850: /* PetscPrintf(comm," **** [%d]%s %d total done, %d nodes (%d local done), min agg. size = %d\n",rank,__FUNCT__,ii,kk/bs,ndone,jj); */

852: #if defined(OUT_AGGS)
853:   if (llev==1) fclose(file);
854: #endif
855:   GAMGTableDestroy(&fgid_flid);
856:   return(0);
857: }

861: static PetscErrorCode PCView_GAMG_AGG(PC pc,PetscViewer viewer)
862: {
864:   PC_MG          *mg      = (PC_MG*)pc->data;
865:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
866:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

869:   PetscViewerASCIIPrintf(viewer,"      AGG specific options\n");
870:   PetscViewerASCIIPrintf(viewer,"        Symmetric graph %s\n",pc_gamg_agg->sym_graph ? "true" : "false");
871:   return(0);
872: }

874: /* -------------------------------------------------------------------------- */
875: /*
876:    PCGAMGGraph_AGG

878:   Input Parameter:
879:    . pc - this
880:    . Amat - matrix on this fine level
881:   Output Parameter:
882:    . a_Gmat -
883: */
886: static PetscErrorCode PCGAMGGraph_AGG(PC pc,Mat Amat,Mat *a_Gmat)
887: {
888:   PetscErrorCode            ierr;
889:   PC_MG                     *mg          = (PC_MG*)pc->data;
890:   PC_GAMG                   *pc_gamg     = (PC_GAMG*)mg->innerctx;
891:   const PetscReal           vfilter      = pc_gamg->threshold;
892:   PC_GAMG_AGG               *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
893:   Mat                       Gmat;
894:   MPI_Comm                  comm;
895:   PetscBool /* set,flg , */ symm;

898:   PetscObjectGetComm((PetscObject)Amat,&comm);
899:   PetscLogEventBegin(PC_GAMGGraph_AGG,0,0,0,0);

901:   /* MatIsSymmetricKnown(Amat, &set, &flg); || !(set && flg) -- this causes lot of symm calls */
902:   symm = (PetscBool)(pc_gamg_agg->sym_graph); /* && !pc_gamg_agg->square_graph; */

904:   PCGAMGCreateGraph(Amat, &Gmat);
905:   PCGAMGFilterGraph(&Gmat, vfilter, symm);

907:   *a_Gmat = Gmat;

909:   PetscLogEventEnd(PC_GAMGGraph_AGG,0,0,0,0);
910:   return(0);
911: }

913: /* -------------------------------------------------------------------------- */
914: /*
915:    PCGAMGCoarsen_AGG

917:   Input Parameter:
918:    . a_pc - this
919:   Input/Output Parameter:
920:    . a_Gmat1 - graph on this fine level - coarsening can change this (squares it)
921:   Output Parameter:
922:    . agg_lists - list of aggregates
923: */
926: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc,Mat *a_Gmat1,PetscCoarsenData **agg_lists)
927: {
929:   PC_MG          *mg          = (PC_MG*)a_pc->data;
930:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
931:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
932:   Mat            mat,Gmat2, Gmat1 = *a_Gmat1;  /* squared graph */
933:   IS             perm;
934:   PetscInt       Istart,Iend,Ii,nloc,bs,n,m;
935:   PetscInt       *permute;
936:   PetscBool      *bIndexSet;
937:   MatCoarsen     crs;
938:   MPI_Comm       comm;
939:   PetscMPIInt    rank;
940:   PetscReal      rr;
941:   PetscInt       iSwapIndex;

944:   PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);
945:   PetscObjectGetComm((PetscObject)Gmat1,&comm);
946:   MPI_Comm_rank(comm, &rank);
947:   MatGetLocalSize(Gmat1, &n, &m);
948:   MatGetBlockSize(Gmat1, &bs);
949:   if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bs %D must be 1",bs);
950:   nloc = n/bs;

952:   if (pc_gamg->current_level < pc_gamg_agg->square_graph) {
953:     PetscInfo2(a_pc,"Square Graph on level %d of %d to square\n",pc_gamg->current_level+1,pc_gamg_agg->square_graph);
954:     MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);
955:   } else Gmat2 = Gmat1;

957:   /* get MIS aggs - randomize */
958:   PetscMalloc1(nloc, &permute);
959:   PetscCalloc1(nloc, &bIndexSet);
960:   for (Ii = 0; Ii < nloc; Ii++) {
961:     permute[Ii]   = Ii;
962:   }
963:   MatGetOwnershipRange(Gmat1, &Istart, &Iend);
964:   for (Ii = 0; Ii < nloc; Ii++) {
965:     PetscRandomGetValueReal(pc_gamg->random,&rr);
966:     iSwapIndex = (PetscInt) (rr*nloc);
967:     if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
968:       PetscInt iTemp = permute[iSwapIndex];
969:       permute[iSwapIndex]   = permute[Ii];
970:       permute[Ii]           = iTemp;
971:       bIndexSet[iSwapIndex] = PETSC_TRUE;
972:     }
973:   }
974:   PetscFree(bIndexSet);

976:   ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);
977: #if defined PETSC_GAMG_USE_LOG
978:   PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);
979: #endif
980:   MatCoarsenCreate(comm, &crs);
981:   MatCoarsenSetFromOptions(crs);
982:   MatCoarsenSetGreedyOrdering(crs, perm);
983:   MatCoarsenSetAdjacency(crs, Gmat2);
984:   MatCoarsenSetStrictAggs(crs, PETSC_TRUE);
985:   MatCoarsenApply(crs);
986:   MatCoarsenGetData(crs, agg_lists); /* output */
987:   MatCoarsenDestroy(&crs);

989:   ISDestroy(&perm);
990:   PetscFree(permute);
991: #if defined PETSC_GAMG_USE_LOG
992:   PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);
993: #endif

995:   /* smooth aggs */
996:   if (Gmat2 != Gmat1) {
997:     const PetscCoarsenData *llist = *agg_lists;
998:     smoothAggs(Gmat2, Gmat1, *agg_lists);
999:     MatDestroy(&Gmat1);
1000:     *a_Gmat1 = Gmat2; /* output */
1001:     PetscCDGetMat(llist, &mat);
1002:     if (mat) SETERRQ(comm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????");
1003:   } else {
1004:     const PetscCoarsenData *llist = *agg_lists;
1005:     /* see if we have a matrix that takes precedence (returned from MatCoarsenApply) */
1006:     PetscCDGetMat(llist, &mat);
1007:     if (mat) {
1008:       MatDestroy(&Gmat1);
1009:       *a_Gmat1 = mat; /* output */
1010:     }
1011:   }
1012:   PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);
1013:   return(0);
1014: }

1016: /* -------------------------------------------------------------------------- */
1017: /*
1018:  PCGAMGProlongator_AGG

1020:  Input Parameter:
1021:  . pc - this
1022:  . Amat - matrix on this fine level
1023:  . Graph - used to get ghost data for nodes in
1024:  . agg_lists - list of aggregates
1025:  Output Parameter:
1026:  . a_P_out - prolongation operator to the next level
1027:  */
1030: static PetscErrorCode PCGAMGProlongator_AGG(PC pc,Mat Amat,Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out)
1031: {
1032:   PC_MG          *mg       = (PC_MG*)pc->data;
1033:   PC_GAMG        *pc_gamg  = (PC_GAMG*)mg->innerctx;
1034:   const PetscInt data_cols = pc_gamg->data_cell_cols;
1036:   PetscInt       Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs;
1037:   Mat            Prol;
1038:   PetscMPIInt    rank, size;
1039:   MPI_Comm       comm;
1040:   const PetscInt col_bs = data_cols;
1041:   PetscReal      *data_w_ghost;
1042:   PetscInt       myCrs0, nbnodes=0, *flid_fgid;
1043:   MatType        mtype;

1046:   PetscObjectGetComm((PetscObject)Amat,&comm);
1047:   PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);
1048:   MPI_Comm_rank(comm, &rank);
1049:   MPI_Comm_size(comm, &size);
1050:   MatGetOwnershipRange(Amat, &Istart, &Iend);
1051:   MatGetBlockSize(Amat, &bs);
1052:   nloc = (Iend-Istart)/bs; my0 = Istart/bs;
1053:   if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) not divisible by bs %D",Iend,Istart,bs);

1055:   /* get 'nLocalSelected' */
1056:   for (ii=0, nLocalSelected = 0; ii < nloc; ii++) {
1057:     PetscBool ise;
1058:     /* filter out singletons 0 or 1? */
1059:     PetscCDEmptyAt(agg_lists, ii, &ise);
1060:     if (!ise) nLocalSelected++;
1061:   }

1063:   /* create prolongator, create P matrix */
1064:   MatGetType(Amat,&mtype);
1065:   MatCreate(comm, &Prol);
1066:   MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);
1067:   MatSetBlockSizes(Prol, bs, col_bs);
1068:   MatSetType(Prol, mtype);
1069:   MatSeqAIJSetPreallocation(Prol, data_cols, NULL);
1070:   MatMPIAIJSetPreallocation(Prol,data_cols, NULL,data_cols, NULL);
1071:   /* nloc*bs, nLocalSelected*col_bs, */
1072:   /* PETSC_DETERMINE, PETSC_DETERMINE, */
1073:   /* data_cols, NULL, data_cols, NULL, */
1074:   /* &Prol); */

1076:   /* can get all points "removed" */
1077:    MatGetSize(Prol, &kk, &ii);
1078:   if (!ii) {
1079:     PetscInfo(pc,"No selected points on coarse grid\n");
1080:     MatDestroy(&Prol);
1081:     *a_P_out = NULL;  /* out */
1082:     PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1083:     return(0);
1084:   }
1085:   PetscInfo1(pc,"New grid %D nodes\n",ii/col_bs);
1086:   MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);

1088:   if ((kk-myCrs0) % col_bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(kk %D -myCrs0 %D) not divisible by col_bs %D",kk,myCrs0,col_bs);
1089:   myCrs0 = myCrs0/col_bs;
1090:   if ((kk/col_bs-myCrs0) != nLocalSelected) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(kk %D/col_bs %D - myCrs0 %D) != nLocalSelected %D)",kk,col_bs,myCrs0,nLocalSelected);

1092:   /* create global vector of data in 'data_w_ghost' */
1093: #if defined PETSC_GAMG_USE_LOG
1094:   PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);
1095: #endif
1096:   if (size > 1) { /*  */
1097:     PetscReal *tmp_gdata,*tmp_ldata,*tp2;
1098:     PetscMalloc1(nloc, &tmp_ldata);
1099:     for (jj = 0; jj < data_cols; jj++) {
1100:       for (kk = 0; kk < bs; kk++) {
1101:         PetscInt        ii,stride;
1102:         const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk;
1103:         for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;

1105:         PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata);

1107:         if (!jj && !kk) { /* now I know how many todal nodes - allocate */
1108:           PetscMalloc1(stride*bs*data_cols, &data_w_ghost);
1109:           nbnodes = bs*stride;
1110:         }
1111:         tp2 = data_w_ghost + jj*bs*stride + kk;
1112:         for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1113:         PetscFree(tmp_gdata);
1114:       }
1115:     }
1116:     PetscFree(tmp_ldata);
1117:   } else {
1118:     nbnodes      = bs*nloc;
1119:     data_w_ghost = (PetscReal*)pc_gamg->data;
1120:   }

1122:   /* get P0 */
1123:   if (size > 1) {
1124:     PetscReal *fid_glid_loc,*fiddata;
1125:     PetscInt  stride;

1127:     PetscMalloc1(nloc, &fid_glid_loc);
1128:     for (kk=0; kk<nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk);
1129:     PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);
1130:     PetscMalloc1(stride, &flid_fgid);
1131:     for (kk=0; kk<stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1132:     PetscFree(fiddata);

1134:     if (stride != nbnodes/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"stride %D != nbnodes %D/bs %D",stride,nbnodes,bs);
1135:     PetscFree(fid_glid_loc);
1136:   } else {
1137:     PetscMalloc1(nloc, &flid_fgid);
1138:     for (kk=0; kk<nloc; kk++) flid_fgid[kk] = my0 + kk;
1139:   }
1140: #if defined PETSC_GAMG_USE_LOG
1141:   PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);
1142:   PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);
1143: #endif
1144:   {
1145:     PetscReal *data_out = NULL;
1146:     formProl0(agg_lists, bs, data_cols, myCrs0, nbnodes,data_w_ghost, flid_fgid, &data_out, Prol);
1147:     PetscFree(pc_gamg->data);

1149:     pc_gamg->data           = data_out;
1150:     pc_gamg->data_cell_rows = data_cols;
1151:     pc_gamg->data_sz        = data_cols*data_cols*nLocalSelected;
1152:   }
1153: #if defined PETSC_GAMG_USE_LOG
1154:   PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);
1155: #endif
1156:   if (size > 1) {PetscFree(data_w_ghost);}
1157:   PetscFree(flid_fgid);

1159:   *a_P_out = Prol;  /* out */

1161:   PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1162:   return(0);
1163: }

1165: /* -------------------------------------------------------------------------- */
1166: /*
1167:    PCGAMGOptProlongator_AGG

1169:   Input Parameter:
1170:    . pc - this
1171:    . Amat - matrix on this fine level
1172:  In/Output Parameter:
1173:    . a_P_out - prolongation operator to the next level
1174: */
1177: static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc,Mat Amat,Mat *a_P)
1178: {
1180:   PC_MG          *mg          = (PC_MG*)pc->data;
1181:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
1182:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
1183:   PetscInt       jj;
1184:   Mat            Prol  = *a_P;
1185:   MPI_Comm       comm;

1188:   PetscObjectGetComm((PetscObject)Amat,&comm);
1189:   PetscLogEventBegin(PC_GAMGOptProlongator_AGG,0,0,0,0);

1191:   /* smooth P0 */
1192:   for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1193:     Mat       tMat;
1194:     Vec       diag;
1195:     PetscReal alpha, emax, emin;
1196: #if defined PETSC_GAMG_USE_LOG
1197:     PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);
1198: #endif
1199:     if (!jj) {
1200:       KSP         eksp;
1201:       Vec         bb, xx;
1202:       PC          epc;

1204:       MatCreateVecs(Amat, &bb, 0);
1205:       MatCreateVecs(Amat, &xx, 0);

1207:       VecSetRandom(bb,pc_gamg->random);

1209:       KSPCreate(comm,&eksp);
1210:       KSPSetErrorIfNotConverged(eksp,pc->erroriffailure);
1211:       KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
1212:       KSPSetNormType(eksp, KSP_NORM_NONE);
1213:       KSPSetOptionsPrefix(eksp,((PetscObject)pc)->prefix);
1214:       KSPAppendOptionsPrefix(eksp, "gamg_est_");
1215:       KSPSetFromOptions(eksp);

1217:       KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);
1218:       KSPSetOperators(eksp, Amat, Amat);
1219:       KSPSetComputeSingularValues(eksp,PETSC_TRUE);

1221:       KSPGetPC(eksp, &epc);
1222:       PCSetType(epc, PCJACOBI);  /* smoother in smoothed agg. */

1224:       /* solve - keep stuff out of logging */
1225:       PetscLogEventDeactivate(KSP_Solve);
1226:       PetscLogEventDeactivate(PC_Apply);
1227:       KSPSolve(eksp, bb, xx);
1228:       PetscLogEventActivate(KSP_Solve);
1229:       PetscLogEventActivate(PC_Apply);

1231:       KSPComputeExtremeSingularValues(eksp, &emax, &emin);
1232:       PetscInfo3(pc,"Smooth P0: max eigen=%e min=%e PC=%s\n",emax,emin,PCJACOBI);
1233:       VecDestroy(&xx);
1234:       VecDestroy(&bb);
1235:       KSPDestroy(&eksp);
1236:     }

1238:     /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1239:     MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);
1240:     MatCreateVecs(Amat, &diag, 0);
1241:     MatGetDiagonal(Amat, diag); /* effectively PCJACOBI */
1242:     VecReciprocal(diag);
1243:     MatDiagonalScale(tMat, diag, 0);
1244:     VecDestroy(&diag);
1245:     alpha = -1.4/emax;
1246:     MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);
1247:     MatDestroy(&Prol);
1248:     Prol  = tMat;
1249: #if defined PETSC_GAMG_USE_LOG
1250:     PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);
1251: #endif
1252:   }
1253:   PetscLogEventEnd(PC_GAMGOptProlongator_AGG,0,0,0,0);
1254:   *a_P = Prol;
1255:   return(0);
1256: }

1258: /* -------------------------------------------------------------------------- */
1259: /*
1260:    PCCreateGAMG_AGG

1262:   Input Parameter:
1263:    . pc -
1264: */
1267: PetscErrorCode  PCCreateGAMG_AGG(PC pc)
1268: {
1270:   PC_MG          *mg      = (PC_MG*)pc->data;
1271:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1272:   PC_GAMG_AGG    *pc_gamg_agg;

1275:   /* create sub context for SA */
1276:   PetscNewLog(pc,&pc_gamg_agg);
1277:   pc_gamg->subctx = pc_gamg_agg;

1279:   pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1280:   pc_gamg->ops->destroy        = PCDestroy_GAMG_AGG;
1281:   /* reset does not do anything; setup not virtual */

1283:   /* set internal function pointers */
1284:   pc_gamg->ops->graph             = PCGAMGGraph_AGG;
1285:   pc_gamg->ops->coarsen           = PCGAMGCoarsen_AGG;
1286:   pc_gamg->ops->prolongator       = PCGAMGProlongator_AGG;
1287:   pc_gamg->ops->optprolongator    = PCGAMGOptProlongator_AGG;
1288:   pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1289:   pc_gamg->ops->view              = PCView_GAMG_AGG;

1291:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNSmooths_C",PCGAMGSetNSmooths_AGG);
1292:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetSymGraph_C",PCGAMGSetSymGraph_AGG);
1293:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetSquareGraph_C",PCGAMGSetSquareGraph_AGG);
1294:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_AGG);
1295:   return(0);
1296: }