Actual source code: agg.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  1: /*
  2:  GAMG geometric-algebric multiogrid PC - Mark Adams 2011
  3:  */

  5:  #include <../src/ksp/pc/impls/gamg/gamg.h>
  6: /*  Next line needed to deactivate KSP_Solve logging */
  7:  #include <petsc/private/kspimpl.h>
  8:  #include <petscblaslapack.h>
  9:  #include <petscdm.h>

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

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

 20:    Not Collective on PC

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

 25:    Options Database Key:
 26: .  -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation

 28:    Level: intermediate

 30:    Concepts: Aggregation AMG preconditioner

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

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

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

 51:   pc_gamg_agg->nsmooths = n;
 52:   return(0);
 53: }

 55: /*@
 56:    PCGAMGSetSymGraph - Symmetrize the graph before computing the aggregation. Some algorithms require the graph be symmetric

 58:    Not Collective on PC

 60:    Input Parameters:
 61: +  pc - the preconditioner context
 62: .  n - PETSC_TRUE or PETSC_FALSE

 64:    Options Database Key:
 65: .  -pc_gamg_sym_graph <true,default=false> - symmetrize the graph before computing the aggregation

 67:    Level: intermediate

 69:    Concepts: Aggregation AMG preconditioner

 71: .seealso: PCGAMGSetSquareGraph()
 72: @*/
 73: PetscErrorCode PCGAMGSetSymGraph(PC pc, PetscBool n)
 74: {

 79:   PetscTryMethod(pc,"PCGAMGSetSymGraph_C",(PC,PetscBool),(pc,n));
 80:   return(0);
 81: }

 83: static PetscErrorCode PCGAMGSetSymGraph_AGG(PC pc, PetscBool n)
 84: {
 85:   PC_MG       *mg          = (PC_MG*)pc->data;
 86:   PC_GAMG     *pc_gamg     = (PC_GAMG*)mg->innerctx;
 87:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

 90:   pc_gamg_agg->sym_graph = n;
 91:   return(0);
 92: }

 94: /*@
 95:    PCGAMGSetSquareGraph -  Square the graph, ie. compute A'*A before aggregating it

 97:    Not Collective on PC

 99:    Input Parameters:
100: +  pc - the preconditioner context
101: -  n - PETSC_TRUE or PETSC_FALSE

103:    Options Database Key:
104: .  -pc_gamg_square_graph <n,default = 1> - number of levels to square the graph on before aggregating it

106:    Level: intermediate

108:    Concepts: Aggregation AMG preconditioner

110: .seealso: PCGAMGSetSymGraph()
111: @*/
112: PetscErrorCode PCGAMGSetSquareGraph(PC pc, PetscInt n)
113: {

118:   PetscTryMethod(pc,"PCGAMGSetSquareGraph_C",(PC,PetscInt),(pc,n));
119:   return(0);
120: }

122: static PetscErrorCode PCGAMGSetSquareGraph_AGG(PC pc, PetscInt n)
123: {
124:   PC_MG       *mg          = (PC_MG*)pc->data;
125:   PC_GAMG     *pc_gamg     = (PC_GAMG*)mg->innerctx;
126:   PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

129:   pc_gamg_agg->square_graph = n;
130:   return(0);
131: }

133: static PetscErrorCode PCSetFromOptions_GAMG_AGG(PetscOptionItems *PetscOptionsObject,PC pc)
134: {
136:   PC_MG          *mg          = (PC_MG*)pc->data;
137:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
138:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

141:   PetscOptionsHead(PetscOptionsObject,"GAMG-AGG options");
142:   {
143:     PetscOptionsInt("-pc_gamg_agg_nsmooths","smoothing steps for smoothed aggregation, usually 1","PCGAMGSetNSmooths",pc_gamg_agg->nsmooths,&pc_gamg_agg->nsmooths,NULL);
144:     PetscOptionsBool("-pc_gamg_sym_graph","Set for asymmetric matrices","PCGAMGSetSymGraph",pc_gamg_agg->sym_graph,&pc_gamg_agg->sym_graph,NULL);
145:     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);
146:   }
147:   PetscOptionsTail();
148:   return(0);
149: }

151: /* -------------------------------------------------------------------------- */
152: static PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
153: {
155:   PC_MG          *mg          = (PC_MG*)pc->data;
156:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;

159:   PetscFree(pc_gamg->subctx);
160:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);
161:   return(0);
162: }

164: /* -------------------------------------------------------------------------- */
165: /*
166:    PCSetCoordinates_AGG
167:      - collective

169:    Input Parameter:
170:    . pc - the preconditioner context
171:    . ndm - dimesion of data (used for dof/vertex for Stokes)
172:    . a_nloc - number of vertices local
173:    . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
174: */

176: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
177: {
178:   PC_MG          *mg      = (PC_MG*)pc->data;
179:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
181:   PetscInt       arrsz,kk,ii,jj,nloc,ndatarows,ndf;
182:   Mat            mat = pc->pmat;

187:   nloc = a_nloc;

189:   /* SA: null space vectors */
190:   MatGetBlockSize(mat, &ndf); /* this does not work for Stokes */
191:   if (coords && ndf==1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
192:   else if (coords) {
193:     if (ndm > ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"degrees of motion %d > block size %d",ndm,ndf);
194:     pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* displacement elasticity */
195:     if (ndm != ndf) {
196:       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);
197:     }
198:   } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
199:   pc_gamg->data_cell_rows = ndatarows = ndf;
200:   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);
201:   arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols;

203:   /* create data - syntactic sugar that should be refactored at some point */
204:   if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) {
205:     PetscFree(pc_gamg->data);
206:     PetscMalloc1(arrsz+1, &pc_gamg->data);
207:   }
208:   /* copy data in - column oriented */
209:   for (kk=0; kk<nloc; kk++) {
210:     const PetscInt M     = nloc*pc_gamg->data_cell_rows; /* stride into data */
211:     PetscReal      *data = &pc_gamg->data[kk*ndatarows]; /* start of cell */
212:     if (pc_gamg->data_cell_cols==1) *data = 1.0;
213:     else {
214:       /* translational modes */
215:       for (ii=0;ii<ndatarows;ii++) {
216:         for (jj=0;jj<ndatarows;jj++) {
217:           if (ii==jj)data[ii*M + jj] = 1.0;
218:           else data[ii*M + jj] = 0.0;
219:         }
220:       }

222:       /* rotational modes */
223:       if (coords) {
224:         if (ndm == 2) {
225:           data   += 2*M;
226:           data[0] = -coords[2*kk+1];
227:           data[1] =  coords[2*kk];
228:         } else {
229:           data   += 3*M;
230:           data[0] = 0.0;             data[M+0] =  coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1];
231:           data[1] = -coords[3*kk+2]; data[M+1] = 0.0;             data[2*M+1] =  coords[3*kk];
232:           data[2] =  coords[3*kk+1]; data[M+2] = -coords[3*kk];   data[2*M+2] = 0.0;
233:         }
234:       }
235:     }
236:   }

238:   pc_gamg->data_sz = arrsz;
239:   return(0);
240: }

242: typedef PetscInt NState;
243: static const NState NOT_DONE=-2;
244: static const NState DELETED =-1;
245: static const NState REMOVED =-3;
246: #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED)

248: /* -------------------------------------------------------------------------- */
249: /*
250:    smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific
251:      - AGG-MG specific: clears singletons out of 'selected_2'

253:    Input Parameter:
254:    . Gmat_2 - glabal matrix of graph (data not defined)   base (squared) graph
255:    . Gmat_1 - base graph to grab with                 base graph
256:    Input/Output Parameter:
257:    . aggs_2 - linked list of aggs with gids)
258: */
259: static PetscErrorCode smoothAggs(PC pc,Mat Gmat_2, Mat Gmat_1,PetscCoarsenData *aggs_2)
260: {
262:   PetscBool      isMPI;
263:   Mat_SeqAIJ     *matA_1, *matB_1=0;
264:   MPI_Comm       comm;
265:   PetscInt       lid,*ii,*idx,ix,Iend,my0,kk,n,j;
266:   Mat_MPIAIJ     *mpimat_2 = 0, *mpimat_1=0;
267:   const PetscInt nloc      = Gmat_2->rmap->n;
268:   PetscScalar    *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid;
269:   PetscInt       *lid_cprowID_1;
270:   NState         *lid_state;
271:   Vec            ghost_par_orig2;

274:   PetscObjectGetComm((PetscObject)Gmat_2,&comm);
275:   MatGetOwnershipRange(Gmat_1,&my0,&Iend);

277:   if (PETSC_FALSE) {
278:     PetscViewer viewer; char fname[32]; static int llev=0;
279:     sprintf(fname,"Gmat2_%d.m",llev++);
280:     PetscViewerASCIIOpen(comm,fname,&viewer);
281:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
282:     MatView(Gmat_2, viewer);
283:     PetscViewerPopFormat(viewer);
284:     PetscViewerDestroy(&viewer);
285:   }

287:   /* get submatrices */
288:   PetscObjectTypeCompare((PetscObject)Gmat_1, MATMPIAIJ, &isMPI);
289:   if (isMPI) {
290:     /* grab matrix objects */
291:     mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data;
292:     mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data;
293:     matA_1   = (Mat_SeqAIJ*)mpimat_1->A->data;
294:     matB_1   = (Mat_SeqAIJ*)mpimat_1->B->data;

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

299:     PetscMalloc1(nloc, &lid_cprowID_1);
300:     for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
301:     for (ix=0; ix<matB_1->compressedrow.nrows; ix++) {
302:       PetscInt lid = matB_1->compressedrow.rindex[ix];
303:       lid_cprowID_1[lid] = ix;
304:     }
305:   } else {
306:     PetscBool        isAIJ;
307:     PetscObjectBaseTypeCompare((PetscObject)Gmat_1,MATSEQAIJ,&isAIJ);
308:     if (!isAIJ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Require AIJ matrix.");
309:     matA_1        = (Mat_SeqAIJ*)Gmat_1->data;
310:     lid_cprowID_1 = NULL;
311:   }
312:   if (nloc>0) {
313:     if (matB_1 && !matB_1->compressedrow.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"matB_1 && !matB_1->compressedrow.use: PETSc bug???");
314:   }
315:   /* get state of locals and selected gid for deleted */
316:   PetscMalloc2(nloc, &lid_state,nloc, &lid_parent_gid);
317:   for (lid = 0; lid < nloc; lid++) {
318:     lid_parent_gid[lid] = -1.0;
319:     lid_state[lid]      = DELETED;
320:   }

322:   /* set lid_state */
323:   for (lid = 0; lid < nloc; lid++) {
324:     PetscCDIntNd *pos;
325:     PetscCDGetHeadPos(aggs_2,lid,&pos);
326:     if (pos) {
327:       PetscInt gid1;

329:       PetscCDIntNdGetID(pos, &gid1);
330:       if (gid1 != lid+my0) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"gid1 %D != lid %D + my0 %D",gid1,lid,my0);
331:       lid_state[lid] = gid1;
332:     }
333:   }

335:   /* map local to selected local, DELETED means a ghost owns it */
336:   for (lid=kk=0; lid<nloc; lid++) {
337:     NState state = lid_state[lid];
338:     if (IS_SELECTED(state)) {
339:       PetscCDIntNd *pos;
340:       PetscCDGetHeadPos(aggs_2,lid,&pos);
341:       while (pos) {
342:         PetscInt gid1;
343:         PetscCDIntNdGetID(pos, &gid1);
344:         PetscCDGetNextPos(aggs_2,lid,&pos);

346:         if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1-my0] = (PetscScalar)(lid + my0);
347:       }
348:     }
349:   }
350:   /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
351:   if (isMPI) {
352:     Vec tempVec;
353:     /* get 'cpcol_1_state' */
354:     MatCreateVecs(Gmat_1, &tempVec, 0);
355:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
356:       PetscScalar v = (PetscScalar)lid_state[kk];
357:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
358:     }
359:     VecAssemblyBegin(tempVec);
360:     VecAssemblyEnd(tempVec);
361:     VecScatterBegin(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
362:     VecScatterEnd(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
363:     VecGetArray(mpimat_1->lvec, &cpcol_1_state);
364:     /* get 'cpcol_2_state' */
365:     VecScatterBegin(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
366:       VecScatterEnd(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
367:     VecGetArray(mpimat_2->lvec, &cpcol_2_state);
368:     /* get 'cpcol_2_par_orig' */
369:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
370:       PetscScalar v = (PetscScalar)lid_parent_gid[kk];
371:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
372:     }
373:     VecAssemblyBegin(tempVec);
374:     VecAssemblyEnd(tempVec);
375:     VecDuplicate(mpimat_2->lvec, &ghost_par_orig2);
376:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
377:     VecScatterEnd(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
378:     VecGetArray(ghost_par_orig2, &cpcol_2_par_orig);

380:     VecDestroy(&tempVec);
381:   } /* ismpi */

383:   /* doit */
384:   for (lid=0; lid<nloc; lid++) {
385:     NState state = lid_state[lid];
386:     if (IS_SELECTED(state)) {
387:       /* steal locals */
388:       ii  = matA_1->i; n = ii[lid+1] - ii[lid];
389:       idx = matA_1->j + ii[lid];
390:       for (j=0; j<n; j++) {
391:         PetscInt lidj   = idx[j], sgid;
392:         NState   statej = lid_state[lidj];
393:         if (statej==DELETED && (sgid=(PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid+my0) { /* steal local */
394:           lid_parent_gid[lidj] = (PetscScalar)(lid+my0); /* send this if sgid is not local */
395:           if (sgid >= my0 && sgid < Iend) {       /* I'm stealing this local from a local sgid */
396:             PetscInt   hav=0,slid=sgid-my0,gidj=lidj+my0;
397:             PetscCDIntNd *pos,*last=NULL;
398:             /* looking for local from local so id_llist_2 works */
399:             PetscCDGetHeadPos(aggs_2,slid,&pos);
400:             while (pos) {
401:               PetscInt gid;
402:               PetscCDIntNdGetID(pos, &gid);
403:               if (gid == gidj) {
404:                 if (!last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null");
405:                 PetscCDRemoveNextNode(aggs_2, slid, last);
406:                 PetscCDAppendNode(aggs_2, lid, pos);
407:                 hav  = 1;
408:                 break;
409:               } else last = pos;

411:               PetscCDGetNextPos(aggs_2,slid,&pos);
412:             }
413:             if (hav!=1) {
414:               if (!hav) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
415:               SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav);
416:             }
417:           } else {            /* I'm stealing this local, owned by a ghost */
418:             if (sgid != -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Have un-symmetric graph (apparently). Use '-%spc_gamg_sym_graph true' to symetrize the graph or '-%spc_gamg_threshold -1' if the matrix is structurally symmetric.",((PetscObject)pc)->prefix,((PetscObject)pc)->prefix);
419:             PetscCDAppendID(aggs_2, lid, lidj+my0);
420:           }
421:         }
422:       } /* local neighbors */
423:     } else if (state == DELETED && lid_cprowID_1) {
424:       PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
425:       /* see if I have a selected ghost neighbor that will steal me */
426:       if ((ix=lid_cprowID_1[lid]) != -1) {
427:         ii  = matB_1->compressedrow.i; n = ii[ix+1] - ii[ix];
428:         idx = matB_1->j + ii[ix];
429:         for (j=0; j<n; j++) {
430:           PetscInt cpid   = idx[j];
431:           NState   statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
432:           if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */
433:             lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
434:             if (sgidold>=my0 && sgidold<Iend) { /* this was mine */
435:               PetscInt   hav=0,oldslidj=sgidold-my0;
436:               PetscCDIntNd *pos,*last=NULL;
437:               /* remove from 'oldslidj' list */
438:               PetscCDGetHeadPos(aggs_2,oldslidj,&pos);
439:               while (pos) {
440:                 PetscInt gid;
441:                 PetscCDIntNdGetID(pos, &gid);
442:                 if (lid+my0 == gid) {
443:                   /* id_llist_2[lastid] = id_llist_2[flid];   /\* remove lid from oldslidj list *\/ */
444:                   if (!last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null");
445:                   PetscCDRemoveNextNode(aggs_2, oldslidj, last);
446:                   /* ghost (PetscScalar)statej will add this later */
447:                   hav = 1;
448:                   break;
449:                 } else last = pos;

451:                 PetscCDGetNextPos(aggs_2,oldslidj,&pos);
452:               }
453:               if (hav!=1) {
454:                 if (!hav) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
455:                 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav);
456:               }
457:             } else {
458:               /* ghosts remove this later */
459:             }
460:           }
461:         }
462:       }
463:     } /* selected/deleted */
464:   } /* node loop */

466:   if (isMPI) {
467:     PetscScalar     *cpcol_2_parent,*cpcol_2_gid;
468:     Vec             tempVec,ghostgids2,ghostparents2;
469:     PetscInt        cpid,nghost_2;
470:     PCGAMGHashTable gid_cpid;

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

475:     /* get 'cpcol_2_parent' */
476:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
477:       VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES);
478:     }
479:     VecAssemblyBegin(tempVec);
480:     VecAssemblyEnd(tempVec);
481:     VecDuplicate(mpimat_2->lvec, &ghostparents2);
482:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
483:     VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
484:     VecGetArray(ghostparents2, &cpcol_2_parent);

486:     /* get 'cpcol_2_gid' */
487:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
488:       PetscScalar v = (PetscScalar)j;
489:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
490:     }
491:     VecAssemblyBegin(tempVec);
492:     VecAssemblyEnd(tempVec);
493:     VecDuplicate(mpimat_2->lvec, &ghostgids2);
494:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
495:     VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
496:     VecGetArray(ghostgids2, &cpcol_2_gid);

498:     VecDestroy(&tempVec);

500:     /* look for deleted ghosts and add to table */
501:     PCGAMGHashTableCreate(2*nghost_2+1, &gid_cpid);
502:     for (cpid = 0; cpid < nghost_2; cpid++) {
503:       NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
504:       if (state==DELETED) {
505:         PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
506:         PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
507:         if (sgid_old == -1 && sgid_new != -1) {
508:           PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
509:           PCGAMGHashTableAdd(&gid_cpid, gid, cpid);
510:         }
511:       }
512:     }

514:     /* look for deleted ghosts and see if they moved - remove it */
515:     for (lid=0; lid<nloc; lid++) {
516:       NState state = lid_state[lid];
517:       if (IS_SELECTED(state)) {
518:         PetscCDIntNd *pos,*last=NULL;
519:         /* look for deleted ghosts and see if they moved */
520:         PetscCDGetHeadPos(aggs_2,lid,&pos);
521:         while (pos) {
522:           PetscInt gid;
523:           PetscCDIntNdGetID(pos, &gid);

525:           if (gid < my0 || gid >= Iend) {
526:             PCGAMGHashTableFind(&gid_cpid, gid, &cpid);
527:             if (cpid != -1) {
528:               /* a moved ghost - */
529:               /* id_llist_2[lastid] = id_llist_2[flid];    /\* remove 'flid' from list *\/ */
530:               PetscCDRemoveNextNode(aggs_2, lid, last);
531:             } else last = pos;
532:           } else last = pos;

534:           PetscCDGetNextPos(aggs_2,lid,&pos);
535:         } /* loop over list of deleted */
536:       } /* selected */
537:     }
538:     PCGAMGHashTableDestroy(&gid_cpid);

540:     /* look at ghosts, see if they changed - and it */
541:     for (cpid = 0; cpid < nghost_2; cpid++) {
542:       PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
543:       if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
544:         PetscInt     gid     = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
545:         PetscInt     slid_new=sgid_new-my0,hav=0;
546:         PetscCDIntNd *pos;

548:         /* search for this gid to see if I have it */
549:         PetscCDGetHeadPos(aggs_2,slid_new,&pos);
550:         while (pos) {
551:           PetscInt gidj;
552:           PetscCDIntNdGetID(pos, &gidj);
553:           PetscCDGetNextPos(aggs_2,slid_new,&pos);

555:           if (gidj == gid) { hav = 1; break; }
556:         }
557:         if (hav != 1) {
558:           /* insert 'flidj' into head of llist */
559:           PetscCDAppendID(aggs_2, slid_new, gid);
560:         }
561:       }
562:     }

564:     VecRestoreArray(mpimat_1->lvec, &cpcol_1_state);
565:     VecRestoreArray(mpimat_2->lvec, &cpcol_2_state);
566:     VecRestoreArray(ghostparents2, &cpcol_2_parent);
567:     VecRestoreArray(ghostgids2, &cpcol_2_gid);
568:     PetscFree(lid_cprowID_1);
569:     VecDestroy(&ghostgids2);
570:     VecDestroy(&ghostparents2);
571:     VecDestroy(&ghost_par_orig2);
572:   }

574:   PetscFree2(lid_state,lid_parent_gid);
575:   return(0);
576: }

578: /* -------------------------------------------------------------------------- */
579: /*
580:    PCSetData_AGG - called if data is not set with PCSetCoordinates.
581:       Looks in Mat for near null space.
582:       Does not work for Stokes

584:   Input Parameter:
585:    . pc -
586:    . a_A - matrix to get (near) null space out of.
587: */
588: static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
589: {
591:   PC_MG          *mg      = (PC_MG*)pc->data;
592:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
593:   MatNullSpace   mnull;

596:   MatGetNearNullSpace(a_A, &mnull);
597:   if (!mnull) {
598:     DM dm;
599:     PCGetDM(pc, &dm);
600:     if (!dm) {
601:       MatGetDM(a_A, &dm);
602:     }
603:     if (dm) {
604:       PetscObject deformation;
605:       PetscInt    Nf;

607:       DMGetNumFields(dm, &Nf);
608:       if (Nf) {
609:         DMGetField(dm, 0, &deformation);
610:         PetscObjectQuery((PetscObject)deformation,"nearnullspace",(PetscObject*)&mnull);
611:         if (!mnull) {
612:           PetscObjectQuery((PetscObject)deformation,"nullspace",(PetscObject*)&mnull);
613:         }
614:       }
615:     }
616:   }

618:   if (!mnull) {
619:     PetscInt bs,NN,MM;
620:     MatGetBlockSize(a_A, &bs);
621:     MatGetLocalSize(a_A, &MM, &NN);
622:     if (MM % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MM %D must be divisible by bs %D",MM,bs);
623:     PCSetCoordinates_AGG(pc, bs, MM/bs, NULL);
624:   } else {
625:     PetscReal *nullvec;
626:     PetscBool has_const;
627:     PetscInt  i,j,mlocal,nvec,bs;
628:     const Vec *vecs; const PetscScalar *v;
629:     MatGetLocalSize(a_A,&mlocal,NULL);
630:     MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
631:     pc_gamg->data_sz = (nvec+!!has_const)*mlocal;
632:     PetscMalloc1((nvec+!!has_const)*mlocal,&nullvec);
633:     if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0;
634:     for (i=0; i<nvec; i++) {
635:       VecGetArrayRead(vecs[i],&v);
636:       for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]);
637:       VecRestoreArrayRead(vecs[i],&v);
638:     }
639:     pc_gamg->data           = nullvec;
640:     pc_gamg->data_cell_cols = (nvec+!!has_const);

642:     MatGetBlockSize(a_A, &bs);

644:     pc_gamg->data_cell_rows = bs;
645:   }
646:   return(0);
647: }

649: /* -------------------------------------------------------------------------- */
650: /*
651:  formProl0

653:    Input Parameter:
654:    . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices 
655:    . bs - row block size
656:    . nSAvec - column bs of new P
657:    . my0crs - global index of start of locals
658:    . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec]
659:    . data_in[data_stride*nSAvec] - local data on fine grid
660:    . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
661:   Output Parameter:
662:    . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
663:    . a_Prol - prolongation operator
664: */
665: 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)
666: {
667:   PetscErrorCode  ierr;
668:   PetscInt        Istart,my0,Iend,nloc,clid,flid = 0,aggID,kk,jj,ii,mm,ndone,nSelected,minsz,nghosts,out_data_stride;
669:   MPI_Comm        comm;
670:   PetscMPIInt     rank;
671:   PetscReal       *out_data;
672:   PetscCDIntNd    *pos;
673:   PCGAMGHashTable fgid_flid;

675: /* #define OUT_AGGS */
676: #if defined(OUT_AGGS)
677:   static PetscInt llev = 0; char fname[32]; FILE *file = NULL; PetscInt pM;
678: #endif

681:   PetscObjectGetComm((PetscObject)a_Prol,&comm);
682:   MPI_Comm_rank(comm,&rank);
683:   MatGetOwnershipRange(a_Prol, &Istart, &Iend);
684:   nloc = (Iend-Istart)/bs; my0 = Istart/bs;
685:   if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Iend %D - Istart %d must be divisible by bs %D",Iend,Istart,bs);
686:   Iend   /= bs;
687:   nghosts = data_stride/bs - nloc;

689:   PCGAMGHashTableCreate(2*nghosts+1, &fgid_flid);
690:   for (kk=0; kk<nghosts; kk++) {
691:     PCGAMGHashTableAdd(&fgid_flid, flid_fgid[nloc+kk], nloc+kk);
692:   }

694: #if defined(OUT_AGGS)
695:   sprintf(fname,"aggs_%d_%d.m",llev++,rank);
696:   if (llev==1) file = fopen(fname,"w");
697:   MatGetSize(a_Prol, &pM, &jj);
698: #endif

700:   /* count selected -- same as number of cols of P */
701:   for (nSelected=mm=0; mm<nloc; mm++) {
702:     PetscBool ise;
703:     PetscCDEmptyAt(agg_llists, mm, &ise);
704:     if (!ise) nSelected++;
705:   }
706:   MatGetOwnershipRangeColumn(a_Prol, &ii, &jj);
707:   if ((ii/nSAvec) != my0crs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D /nSAvec %D  != my0crs %D",ii,nSAvec,my0crs);
708:   if (nSelected != (jj-ii)/nSAvec) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nSelected %D != (jj %D - ii %D)/nSAvec %D",nSelected,jj,ii,nSAvec);

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

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

717:   /* find points and set prolongation */
718:   minsz = 100;
719:   ndone = 0;
720:   for (mm = clid = 0; mm < nloc; mm++) {
721:     PetscCDSizeAt(agg_llists, mm, &jj);
722:     if (jj > 0) {
723:       const PetscInt lid = mm, cgid = my0crs + clid;
724:       PetscInt       cids[100]; /* max bs */
725:       PetscBLASInt   asz  =jj,M=asz*bs,N=nSAvec,INFO;
726:       PetscBLASInt   Mdata=M+((N-M>0) ? N-M : 0),LDA=Mdata,LWORK=N*bs;
727:       PetscScalar    *qqc,*qqr,*TAU,*WORK;
728:       PetscInt       *fids;
729:       PetscReal      *data;
730:       /* count agg */
731:       if (asz<minsz) minsz = asz;

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

736:       aggID = 0;
737:       PetscCDGetHeadPos(agg_llists,lid,&pos);
738:       while (pos) {
739:         PetscInt gid1;
740:         PetscCDIntNdGetID(pos, &gid1);
741:         PetscCDGetNextPos(agg_llists,lid,&pos);

743:         if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
744:         else {
745:           PCGAMGHashTableFind(&fgid_flid, gid1, &flid);
746:           if (flid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot find gid1 in table");
747:         }
748:         /* copy in B_i matrix - column oriented */
749:         data = &data_in[flid*bs];
750:         for (ii = 0; ii < bs; ii++) {
751:           for (jj = 0; jj < N; jj++) {
752:             PetscReal d = data[jj*data_stride + ii];
753:             qqc[jj*Mdata + aggID*bs + ii] = d;
754:           }
755:         }
756: #if defined(OUT_AGGS)
757:         if (llev==1) {
758:           char     str[] = "plot(%e,%e,'r*'), hold on,\n", col[] = "rgbkmc", sim[] = "*os+h>d<vx^";
759:           PetscInt MM,pi,pj;
760:           str[12] = col[(clid+17*rank)%6]; str[13] = sim[(clid+17*rank)%11];
761:           MM      = (PetscInt)(PetscSqrtReal((PetscReal)pM));
762:           pj      = gid1/MM; pi = gid1%MM;
763:           fprintf(file,str,(double)pi,(double)pj);
764:           /* fprintf(file,str,data[2*data_stride+1],-data[2*data_stride]); */
765:         }
766: #endif
767:         /* set fine IDs */
768:         for (kk=0; kk<bs; kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk;

770:         aggID++;
771:       }

773:       /* pad with zeros */
774:       for (ii = asz*bs; ii < Mdata; ii++) {
775:         for (jj = 0; jj < N; jj++, kk++) {
776:           qqc[jj*Mdata + ii] = .0;
777:         }
778:       }

780:       ndone += aggID;
781:       /* QR */
782:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
783:       PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
784:       PetscFPTrapPop();
785:       if (INFO != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xGEQRF error");
786:       /* get R - column oriented - output B_{i+1} */
787:       {
788:         PetscReal *data = &out_data[clid*nSAvec];
789:         for (jj = 0; jj < nSAvec; jj++) {
790:           for (ii = 0; ii < nSAvec; ii++) {
791:             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);
792:            if (ii <= jj) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]);
793:            else data[jj*out_data_stride + ii] = 0.;
794:           }
795:         }
796:       }

798:       /* get Q - row oriented */
799:       PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
800:       if (INFO != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xORGQR error arg %d",-INFO);

802:       for (ii = 0; ii < M; ii++) {
803:         for (jj = 0; jj < N; jj++) {
804:           qqr[N*ii + jj] = qqc[jj*Mdata + ii];
805:         }
806:       }

808:       /* add diagonal block of P0 */
809:       for (kk=0; kk<N; kk++) {
810:         cids[kk] = N*cgid + kk; /* global col IDs in P0 */
811:       }
812:       MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);

814:       PetscFree5(qqc,qqr,TAU,WORK,fids);
815:       clid++;
816:     } /* coarse agg */
817:   } /* for all fine nodes */
818:   MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
819:   MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);

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

826: #if defined(OUT_AGGS)
827:   if (llev==1) fclose(file);
828: #endif
829:   PCGAMGHashTableDestroy(&fgid_flid);
830:   return(0);
831: }

833: static PetscErrorCode PCView_GAMG_AGG(PC pc,PetscViewer viewer)
834: {
836:   PC_MG          *mg      = (PC_MG*)pc->data;
837:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
838:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;

841:   PetscViewerASCIIPrintf(viewer,"      AGG specific options\n");
842:   PetscViewerASCIIPrintf(viewer,"        Symmetric graph %s\n",pc_gamg_agg->sym_graph ? "true" : "false");
843:   PetscViewerASCIIPrintf(viewer,"        Number of levels to square graph %D\n",pc_gamg_agg->square_graph);
844:   PetscViewerASCIIPrintf(viewer,"        Number smoothing steps %D\n",pc_gamg_agg->nsmooths);
845:   return(0);
846: }

848: /* -------------------------------------------------------------------------- */
849: /*
850:    PCGAMGGraph_AGG

852:   Input Parameter:
853:    . pc - this
854:    . Amat - matrix on this fine level
855:   Output Parameter:
856:    . a_Gmat -
857: */
858: static PetscErrorCode PCGAMGGraph_AGG(PC pc,Mat Amat,Mat *a_Gmat)
859: {
860:   PetscErrorCode            ierr;
861:   PC_MG                     *mg          = (PC_MG*)pc->data;
862:   PC_GAMG                   *pc_gamg     = (PC_GAMG*)mg->innerctx;
863:   const PetscReal           vfilter      = pc_gamg->threshold[pc_gamg->current_level];
864:   PC_GAMG_AGG               *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
865:   Mat                       Gmat;
866:   MPI_Comm                  comm;
867:   PetscBool /* set,flg , */ symm;

870:   PetscObjectGetComm((PetscObject)Amat,&comm);
871:   PetscLogEventBegin(PC_GAMGGraph_AGG,0,0,0,0);

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

876:   PCGAMGCreateGraph(Amat, &Gmat);
877:   PCGAMGFilterGraph(&Gmat, vfilter, symm);

879:   *a_Gmat = Gmat;

881:   PetscLogEventEnd(PC_GAMGGraph_AGG,0,0,0,0);
882:   return(0);
883: }

885: /* -------------------------------------------------------------------------- */
886: /*
887:    PCGAMGCoarsen_AGG

889:   Input Parameter:
890:    . a_pc - this
891:   Input/Output Parameter:
892:    . a_Gmat1 - graph on this fine level - coarsening can change this (squares it)
893:   Output Parameter:
894:    . agg_lists - list of aggregates
895: */
896: static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc,Mat *a_Gmat1,PetscCoarsenData **agg_lists)
897: {
899:   PC_MG          *mg          = (PC_MG*)a_pc->data;
900:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
901:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
902:   Mat            mat,Gmat2, Gmat1 = *a_Gmat1;  /* squared graph */
903:   IS             perm;
904:   PetscInt       Istart,Iend,Ii,nloc,bs,n,m;
905:   PetscInt       *permute;
906:   PetscBool      *bIndexSet;
907:   MatCoarsen     crs;
908:   MPI_Comm       comm;
909:   PetscMPIInt    rank;
910:   PetscReal      hashfact;
911:   PetscInt       iSwapIndex;
912:   PetscRandom    random;

915:   PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);
916:   PetscObjectGetComm((PetscObject)Gmat1,&comm);
917:   MPI_Comm_rank(comm, &rank);
918:   MatGetLocalSize(Gmat1, &n, &m);
919:   MatGetBlockSize(Gmat1, &bs);
920:   if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bs %D must be 1",bs);
921:   nloc = n/bs;

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

928:   /* get MIS aggs - randomize */
929:   PetscMalloc1(nloc, &permute);
930:   PetscCalloc1(nloc, &bIndexSet);
931:   for (Ii = 0; Ii < nloc; Ii++) {
932:     permute[Ii]   = Ii;
933:   }
934:   PetscRandomCreate(PETSC_COMM_SELF,&random);
935:   MatGetOwnershipRange(Gmat1, &Istart, &Iend);
936:   for (Ii = 0; Ii < nloc; Ii++) {
937:     PetscRandomGetValueReal(random,&hashfact);
938:     iSwapIndex = (PetscInt) (hashfact*nloc)%nloc;
939:     if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
940:       PetscInt iTemp = permute[iSwapIndex];
941:       permute[iSwapIndex]   = permute[Ii];
942:       permute[Ii]           = iTemp;
943:       bIndexSet[iSwapIndex] = PETSC_TRUE;
944:     }
945:   }
946:   PetscFree(bIndexSet);
947:   PetscRandomDestroy(&random);
948:   ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);
949: #if defined PETSC_GAMG_USE_LOG
950:   PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);
951: #endif
952:   MatCoarsenCreate(comm, &crs);
953:   MatCoarsenSetFromOptions(crs);
954:   MatCoarsenSetGreedyOrdering(crs, perm);
955:   MatCoarsenSetAdjacency(crs, Gmat2);
956:   MatCoarsenSetStrictAggs(crs, PETSC_TRUE);
957:   MatCoarsenApply(crs);
958:   MatCoarsenGetData(crs, agg_lists); /* output */
959:   MatCoarsenDestroy(&crs);

961:   ISDestroy(&perm);
962:   PetscFree(permute);
963: #if defined PETSC_GAMG_USE_LOG
964:   PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);
965: #endif

967:   /* smooth aggs */
968:   if (Gmat2 != Gmat1) {
969:     const PetscCoarsenData *llist = *agg_lists;
970:     smoothAggs(a_pc,Gmat2, Gmat1, *agg_lists);
971:     MatDestroy(&Gmat1);
972:     *a_Gmat1 = Gmat2; /* output */
973:     PetscCDGetMat(llist, &mat);
974:     if (mat) SETERRQ(comm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????");
975:   } else {
976:     const PetscCoarsenData *llist = *agg_lists;
977:     /* see if we have a matrix that takes precedence (returned from MatCoarsenApply) */
978:     PetscCDGetMat(llist, &mat);
979:     if (mat) {
980:       MatDestroy(&Gmat1);
981:       *a_Gmat1 = mat; /* output */
982:     }
983:   }
984:   PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);
985:   return(0);
986: }

988: /* -------------------------------------------------------------------------- */
989: /*
990:  PCGAMGProlongator_AGG

992:  Input Parameter:
993:  . pc - this
994:  . Amat - matrix on this fine level
995:  . Graph - used to get ghost data for nodes in
996:  . agg_lists - list of aggregates
997:  Output Parameter:
998:  . a_P_out - prolongation operator to the next level
999:  */
1000: static PetscErrorCode PCGAMGProlongator_AGG(PC pc,Mat Amat,Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out)
1001: {
1002:   PC_MG          *mg       = (PC_MG*)pc->data;
1003:   PC_GAMG        *pc_gamg  = (PC_GAMG*)mg->innerctx;
1004:   const PetscInt col_bs = pc_gamg->data_cell_cols;
1006:   PetscInt       Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs;
1007:   Mat            Prol;
1008:   PetscMPIInt    rank, size;
1009:   MPI_Comm       comm;
1010:   PetscReal      *data_w_ghost;
1011:   PetscInt       myCrs0, nbnodes=0, *flid_fgid;
1012:   MatType        mtype;

1015:   PetscObjectGetComm((PetscObject)Amat,&comm);
1016:   if (col_bs < 1) SETERRQ(comm,PETSC_ERR_PLIB,"Column bs cannot be less than 1");
1017:   PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);
1018:   MPI_Comm_rank(comm, &rank);
1019:   MPI_Comm_size(comm, &size);
1020:   MatGetOwnershipRange(Amat, &Istart, &Iend);
1021:   MatGetBlockSize(Amat, &bs);
1022:   nloc = (Iend-Istart)/bs; my0 = Istart/bs;
1023:   if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) not divisible by bs %D",Iend,Istart,bs);

1025:   /* get 'nLocalSelected' */
1026:   for (ii=0, nLocalSelected = 0; ii < nloc; ii++) {
1027:     PetscBool ise;
1028:     /* filter out singletons 0 or 1? */
1029:     PetscCDEmptyAt(agg_lists, ii, &ise);
1030:     if (!ise) nLocalSelected++;
1031:   }

1033:   /* create prolongator, create P matrix */
1034:   MatGetType(Amat,&mtype);
1035:   MatCreate(comm, &Prol);
1036:   MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);
1037:   MatSetBlockSizes(Prol, bs, col_bs);
1038:   MatSetType(Prol, mtype);
1039:   MatSeqAIJSetPreallocation(Prol, col_bs, NULL);
1040:   MatMPIAIJSetPreallocation(Prol,col_bs, NULL,col_bs, NULL);

1042:   /* can get all points "removed" */
1043:    MatGetSize(Prol, &kk, &ii);
1044:   if (!ii) {
1045:     PetscInfo(pc,"No selected points on coarse grid\n");
1046:     MatDestroy(&Prol);
1047:     *a_P_out = NULL;  /* out */
1048:     PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1049:     return(0);
1050:   }
1051:   PetscInfo1(pc,"New grid %D nodes\n",ii/col_bs);
1052:   MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);

1054:   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);
1055:   myCrs0 = myCrs0/col_bs;
1056:   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);

1058:   /* create global vector of data in 'data_w_ghost' */
1059: #if defined PETSC_GAMG_USE_LOG
1060:   PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);
1061: #endif
1062:   if (size > 1) { /*  */
1063:     PetscReal *tmp_gdata,*tmp_ldata,*tp2;
1064:     PetscMalloc1(nloc, &tmp_ldata);
1065:     for (jj = 0; jj < col_bs; jj++) {
1066:       for (kk = 0; kk < bs; kk++) {
1067:         PetscInt        ii,stride;
1068:         const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk;
1069:         for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;

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

1073:         if (!jj && !kk) { /* now I know how many todal nodes - allocate */
1074:           PetscMalloc1(stride*bs*col_bs, &data_w_ghost);
1075:           nbnodes = bs*stride;
1076:         }
1077:         tp2 = data_w_ghost + jj*bs*stride + kk;
1078:         for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1079:         PetscFree(tmp_gdata);
1080:       }
1081:     }
1082:     PetscFree(tmp_ldata);
1083:   } else {
1084:     nbnodes      = bs*nloc;
1085:     data_w_ghost = (PetscReal*)pc_gamg->data;
1086:   }

1088:   /* get P0 */
1089:   if (size > 1) {
1090:     PetscReal *fid_glid_loc,*fiddata;
1091:     PetscInt  stride;

1093:     PetscMalloc1(nloc, &fid_glid_loc);
1094:     for (kk=0; kk<nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk);
1095:     PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);
1096:     PetscMalloc1(stride, &flid_fgid);
1097:     for (kk=0; kk<stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1098:     PetscFree(fiddata);

1100:     if (stride != nbnodes/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"stride %D != nbnodes %D/bs %D",stride,nbnodes,bs);
1101:     PetscFree(fid_glid_loc);
1102:   } else {
1103:     PetscMalloc1(nloc, &flid_fgid);
1104:     for (kk=0; kk<nloc; kk++) flid_fgid[kk] = my0 + kk;
1105:   }
1106: #if defined PETSC_GAMG_USE_LOG
1107:   PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);
1108:   PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);
1109: #endif
1110:   {
1111:     PetscReal *data_out = NULL;
1112:     formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes,data_w_ghost, flid_fgid, &data_out, Prol);
1113:     PetscFree(pc_gamg->data);

1115:     pc_gamg->data           = data_out;
1116:     pc_gamg->data_cell_rows = col_bs;
1117:     pc_gamg->data_sz        = col_bs*col_bs*nLocalSelected;
1118:   }
1119: #if defined PETSC_GAMG_USE_LOG
1120:   PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);
1121: #endif
1122:   if (size > 1) {PetscFree(data_w_ghost);}
1123:   PetscFree(flid_fgid);

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

1127:   PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1128:   return(0);
1129: }

1131: /* -------------------------------------------------------------------------- */
1132: /*
1133:    PCGAMGOptProlongator_AGG

1135:   Input Parameter:
1136:    . pc - this
1137:    . Amat - matrix on this fine level
1138:  In/Output Parameter:
1139:    . a_P - prolongation operator to the next level
1140: */
1141: static PetscErrorCode PCGAMGOptProlongator_AGG(PC pc,Mat Amat,Mat *a_P)
1142: {
1144:   PC_MG          *mg          = (PC_MG*)pc->data;
1145:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
1146:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
1147:   PetscInt       jj;
1148:   Mat            Prol  = *a_P;
1149:   MPI_Comm       comm;
1150:   KSP            eksp;
1151:   Vec            bb, xx;
1152:   PC             epc;
1153:   PetscReal      alpha, emax, emin;
1154:   PetscRandom    random;

1157:   PetscObjectGetComm((PetscObject)Amat,&comm);
1158:   PetscLogEventBegin(PC_GAMGOptProlongator_AGG,0,0,0,0);

1160:   /* compute maximum value of operator to be used in smoother */
1161:   if (0 < pc_gamg_agg->nsmooths) {
1162:     MatCreateVecs(Amat, &bb, 0);
1163:     MatCreateVecs(Amat, &xx, 0);
1164:     PetscRandomCreate(PETSC_COMM_SELF,&random);
1165:     VecSetRandom(bb,random);
1166:     PetscRandomDestroy(&random);

1168:     KSPCreate(comm,&eksp);
1169:     KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
1170:     KSPSetNormType(eksp, KSP_NORM_NONE);
1171:     KSPSetErrorIfNotConverged(eksp,PETSC_FALSE);

1173:     KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);
1174:     KSPSetOperators(eksp, Amat, Amat);
1175:     KSPSetComputeSingularValues(eksp,PETSC_TRUE);

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

1180:     KSPSetOptionsPrefix(eksp,((PetscObject)pc)->prefix);
1181:     KSPAppendOptionsPrefix(eksp, "gamg_est_");
1182:     KSPSetFromOptions(eksp);

1184:     /* solve - keep stuff out of logging */
1185:     PetscLogEventDeactivate(KSP_Solve);
1186:     PetscLogEventDeactivate(PC_Apply);
1187:     KSPSolve(eksp, bb, xx);
1188:     PetscLogEventActivate(KSP_Solve);
1189:     PetscLogEventActivate(PC_Apply);

1191:     KSPComputeExtremeSingularValues(eksp, &emax, &emin);
1192:     PetscInfo3(pc,"Smooth P0: max eigen=%e min=%e PC=%s\n",emax,emin,PCJACOBI);
1193:     VecDestroy(&xx);
1194:     VecDestroy(&bb);
1195:     KSPDestroy(&eksp);
1196:   }

1198:   /* smooth P0 */
1199:   for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1200:     Mat       tMat;
1201:     Vec       diag;

1203: #if defined PETSC_GAMG_USE_LOG
1204:     PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);
1205: #endif

1207:     /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1208:     MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);
1209:     MatCreateVecs(Amat, &diag, 0);
1210:     MatGetDiagonal(Amat, diag); /* effectively PCJACOBI */
1211:     VecReciprocal(diag);
1212:     MatDiagonalScale(tMat, diag, 0);
1213:     VecDestroy(&diag);
1214:     alpha = -1.4/emax;
1215:     MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);
1216:     MatDestroy(&Prol);
1217:     Prol  = tMat;
1218: #if defined PETSC_GAMG_USE_LOG
1219:     PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);
1220: #endif
1221:   }
1222:   PetscLogEventEnd(PC_GAMGOptProlongator_AGG,0,0,0,0);
1223:   *a_P = Prol;
1224:   return(0);
1225: }

1227: /* -------------------------------------------------------------------------- */
1228: /*
1229:    PCCreateGAMG_AGG

1231:   Input Parameter:
1232:    . pc -
1233: */
1234: PetscErrorCode  PCCreateGAMG_AGG(PC pc)
1235: {
1237:   PC_MG          *mg      = (PC_MG*)pc->data;
1238:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1239:   PC_GAMG_AGG    *pc_gamg_agg;

1242:   /* create sub context for SA */
1243:   PetscNewLog(pc,&pc_gamg_agg);
1244:   pc_gamg->subctx = pc_gamg_agg;

1246:   pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1247:   pc_gamg->ops->destroy        = PCDestroy_GAMG_AGG;
1248:   /* reset does not do anything; setup not virtual */

1250:   /* set internal function pointers */
1251:   pc_gamg->ops->graph             = PCGAMGGraph_AGG;
1252:   pc_gamg->ops->coarsen           = PCGAMGCoarsen_AGG;
1253:   pc_gamg->ops->prolongator       = PCGAMGProlongator_AGG;
1254:   pc_gamg->ops->optprolongator    = PCGAMGOptProlongator_AGG;
1255:   pc_gamg->ops->createdefaultdata = PCSetData_AGG;
1256:   pc_gamg->ops->view              = PCView_GAMG_AGG;

1258:   pc_gamg_agg->square_graph = 1;
1259:   pc_gamg_agg->sym_graph    = PETSC_FALSE;
1260:   pc_gamg_agg->nsmooths     = 1;


1263:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNSmooths_C",PCGAMGSetNSmooths_AGG);
1264:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetSymGraph_C",PCGAMGSetSymGraph_AGG);
1265:   PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetSquareGraph_C",PCGAMGSetSquareGraph_AGG);
1266:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_AGG);
1267:   return(0);
1268: }