Actual source code: agg.c

petsc-3.4.4 2014-03-13
  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:   PetscBool 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: PetscErrorCode PCGAMGSetNSmooths_GAMG(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: PetscErrorCode PCGAMGSetSymGraph_GAMG(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, PetscBool n)
120: {

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

131: PetscErrorCode PCGAMGSetSquareGraph_GAMG(PC pc, PetscBool 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: PetscErrorCode PCSetFromOptions_GAMG_AGG(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;
157:   PetscBool      flag;


161:   PetscOptionsHead("GAMG-AGG options");
162:   {
163:     /* -pc_gamg_agg_nsmooths */
164:     pc_gamg_agg->nsmooths = 1;

166:     PetscOptionsInt("-pc_gamg_agg_nsmooths",
167:                            "smoothing steps for smoothed aggregation, usually 1 (0)",
168:                            "PCGAMGSetNSmooths",
169:                            pc_gamg_agg->nsmooths,
170:                            &pc_gamg_agg->nsmooths,
171:                            &flag);

173:     /* -pc_gamg_sym_graph */
174:     pc_gamg_agg->sym_graph = PETSC_FALSE;

176:     PetscOptionsBool("-pc_gamg_sym_graph",
177:                             "Set for asymmetric matrices",
178:                             "PCGAMGSetSymGraph",
179:                             pc_gamg_agg->sym_graph,
180:                             &pc_gamg_agg->sym_graph,
181:                             &flag);

183:     /* -pc_gamg_square_graph */
184:     pc_gamg_agg->square_graph = PETSC_TRUE;

186:     PetscOptionsBool("-pc_gamg_square_graph",
187:                             "For faster coarsening and lower coarse grid complexity",
188:                             "PCGAMGSetSquareGraph",
189:                             pc_gamg_agg->square_graph,
190:                             &pc_gamg_agg->square_graph,
191:                             &flag);
192:   }
193:   PetscOptionsTail();
194:   return(0);
195: }

197: /* -------------------------------------------------------------------------- */
198: /*
199:    PCDestroy_AGG

201:   Input Parameter:
202:    . pc -
203: */
206: PetscErrorCode PCDestroy_GAMG_AGG(PC pc)
207: {
209:   PC_MG          *mg          = (PC_MG*)pc->data;
210:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;

213:   PetscFree(pc_gamg->subctx);
214:   return(0);
215: }

217: /* -------------------------------------------------------------------------- */
218: /*
219:    PCSetCoordinates_AGG
220:      - collective

222:    Input Parameter:
223:    . pc - the preconditioner context
224:    . ndm - dimesion of data (used for dof/vertex for Stokes)
225:    . a_nloc - number of vertices local
226:    . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
227: */

231: static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords)
232: {
233:   PC_MG          *mg      = (PC_MG*)pc->data;
234:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
236:   PetscInt       arrsz,kk,ii,jj,nloc,ndatarows,ndf;
237:   Mat            mat = pc->pmat;

242:   nloc = a_nloc;

244:   /* SA: null space vectors */
245:   MatGetBlockSize(mat, &ndf); /* this does not work for Stokes */
246:   if (coords && ndf==1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
247:   else if (coords) {
248:     if (ndm > ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"degrees of motion %d > block size %d",ndm,ndf);
249:     pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* displacement elasticity */
250:     if (ndm != ndf) {
251:       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);
252:     }
253:   } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
254:   pc_gamg->data_cell_rows = ndatarows = ndf;
255:   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);
256:   arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols;

258:   /* create data - syntactic sugar that should be refactored at some point */
259:   if (pc_gamg->data==0 || (pc_gamg->data_sz != arrsz)) {
260:     PetscFree(pc_gamg->data);
261:     PetscMalloc((arrsz+1)*sizeof(PetscReal), &pc_gamg->data);
262:   }
263:   /* copy data in - column oriented */
264:   for (kk=0; kk<nloc; kk++) {
265:     const PetscInt M     = nloc*pc_gamg->data_cell_rows; /* stride into data */
266:     PetscReal      *data = &pc_gamg->data[kk*ndatarows]; /* start of cell */
267:     if (pc_gamg->data_cell_cols==1) *data = 1.0;
268:     else {
269:       /* translational modes */
270:       for (ii=0;ii<ndatarows;ii++) {
271:         for (jj=0;jj<ndatarows;jj++) {
272:           if (ii==jj)data[ii*M + jj] = 1.0;
273:           else data[ii*M + jj] = 0.0;
274:         }
275:       }

277:       /* rotational modes */
278:       if (coords) {
279:         if (ndm == 2) {
280:           data   += 2*M;
281:           data[0] = -coords[2*kk+1];
282:           data[1] =  coords[2*kk];
283:         } else {
284:           data   += 3*M;
285:           data[0] = 0.0;             data[M+0] =  coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1];
286:           data[1] = -coords[3*kk+2]; data[M+1] = 0.0;             data[2*M+1] =  coords[3*kk];
287:           data[2] =  coords[3*kk+1]; data[M+2] = -coords[3*kk];   data[2*M+2] = 0.0;
288:         }
289:       }
290:     }
291:   }

293:   pc_gamg->data_sz = arrsz;
294:   return(0);
295: }

297: typedef PetscInt NState;
298: static const NState NOT_DONE=-2;
299: static const NState DELETED =-1;
300: static const NState REMOVED =-3;
301: #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED)

303: /* -------------------------------------------------------------------------- */
304: /*
305:    smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific
306:      - AGG-MG specific: clears singletons out of 'selected_2'

308:    Input Parameter:
309:    . Gmat_2 - glabal matrix of graph (data not defined)
310:    . Gmat_1 - base graph to grab with
311:    Input/Output Parameter:
312:    . aggs_2 - linked list of aggs with gids)
313: */
316: static PetscErrorCode smoothAggs(const Mat Gmat_2, /* base (squared) graph */
317:                                  const Mat Gmat_1,  /* base graph */
318:                                  /* const IS selected_2, [nselected local] selected vertices */
319:                                  PetscCoarsenData *aggs_2)  /* [nselected local] global ID of aggregate */
320: {
322:   PetscBool      isMPI;
323:   Mat_SeqAIJ     *matA_1, *matB_1=0, *matA_2, *matB_2=0;
324:   MPI_Comm       comm;
325:   PetscMPIInt    rank,size;
326:   PetscInt       lid,*ii,*idx,ix,Iend,my0,kk,n,j;
327:   Mat_MPIAIJ     *mpimat_2 = 0, *mpimat_1=0;
328:   const PetscInt nloc      = Gmat_2->rmap->n;
329:   PetscScalar    *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid;
330:   PetscInt       *lid_cprowID_1;
331:   NState         *lid_state;
332:   Vec            ghost_par_orig2;

335:   PetscObjectGetComm((PetscObject)Gmat_2,&comm);
336:   MPI_Comm_rank(comm, &rank);
337:   MPI_Comm_size(comm, &size);
338:   MatGetOwnershipRange(Gmat_1,&my0,&Iend);

340:   if (PETSC_FALSE) {
341:     PetscViewer viewer; char fname[32]; static int llev=0;
342:     sprintf(fname,"Gmat2_%d.m",llev++);
343:     PetscViewerASCIIOpen(comm,fname,&viewer);
344:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
345:     MatView(Gmat_2, viewer);
346:     PetscViewerDestroy(&viewer);
347:   }

349:   /* get submatrices */
350:   PetscObjectTypeCompare((PetscObject)Gmat_1, MATMPIAIJ, &isMPI);
351:   if (isMPI) {
352:     /* grab matrix objects */
353:     mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data;
354:     mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data;
355:     matA_1   = (Mat_SeqAIJ*)mpimat_1->A->data;
356:     matB_1   = (Mat_SeqAIJ*)mpimat_1->B->data;
357:     matA_2   = (Mat_SeqAIJ*)mpimat_2->A->data;
358:     matB_2   = (Mat_SeqAIJ*)mpimat_2->B->data;

360:     /* force compressed row storage for B matrix in AuxMat */
361:     matB_1->compressedrow.check = PETSC_TRUE;

363:     MatCheckCompressedRow(mpimat_1->B,&matB_1->compressedrow,matB_1->i,Gmat_1->rmap->n,-1.0);

365:     PetscMalloc(nloc*sizeof(PetscInt), &lid_cprowID_1);
366:     for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1;
367:     for (ix=0; ix<matB_1->compressedrow.nrows; ix++) {
368:       PetscInt lid = matB_1->compressedrow.rindex[ix];
369:       lid_cprowID_1[lid] = ix;
370:     }
371:   } else {
372:     matA_1        = (Mat_SeqAIJ*)Gmat_1->data;
373:     matA_2        = (Mat_SeqAIJ*)Gmat_2->data;
374:     lid_cprowID_1 = NULL;
375:   }
376:   if (!(matA_1 && !matA_1->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matA_1 && !matA_1->compressedrow.use)");
377:   if (!(matB_1==0 || matB_1->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matB_1==0 || matB_1->compressedrow.use)");
378:   if (!(matA_2 && !matA_2->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matA_2 && !matA_2->compressedrow.use)");
379:   if (!(matB_2==0 || matB_2->compressedrow.use)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"!(matB_2==0 || matB_2->compressedrow.use)");

381:   /* get state of locals and selected gid for deleted */
382:   PetscMalloc(nloc*sizeof(NState), &lid_state);
383:   PetscMalloc(nloc*sizeof(PetscScalar), &lid_parent_gid);
384:   for (lid = 0; lid < nloc; lid++) {
385:     lid_parent_gid[lid] = -1.0;
386:     lid_state[lid]      = DELETED;
387:   }

389:   /* set lid_state */
390:   for (lid = 0; lid < nloc; lid++) {
391:     PetscCDPos pos;
392:     PetscCDGetHeadPos(aggs_2,lid,&pos);
393:     if (pos) {
394:       PetscInt gid1;

396:       PetscLLNGetID(pos, &gid1);
397:       if (gid1 != lid+my0) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"gid1 %D != lid %D + my0 %D",gid1,lid,my0);
398:       lid_state[lid] = gid1;
399:     }
400:   }

402:   /* map local to selected local, DELETED means a ghost owns it */
403:   for (lid=kk=0; lid<nloc; lid++) {
404:     NState state = lid_state[lid];
405:     if (IS_SELECTED(state)) {
406:       PetscCDPos pos;
407:       PetscCDGetHeadPos(aggs_2,lid,&pos);
408:       while (pos) {
409:         PetscInt gid1;
410:         PetscLLNGetID(pos, &gid1);
411:         PetscCDGetNextPos(aggs_2,lid,&pos);

413:         if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1-my0] = (PetscScalar)(lid + my0);
414:       }
415:     }
416:   }
417:   /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
418:   if (isMPI) {
419:     Vec tempVec;
420:     /* get 'cpcol_1_state' */
421:     MatGetVecs(Gmat_1, &tempVec, 0);
422:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
423:       PetscScalar v = (PetscScalar)lid_state[kk];
424:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
425:     }
426:     VecAssemblyBegin(tempVec);
427:     VecAssemblyEnd(tempVec);
428:     VecScatterBegin(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
429:       VecScatterEnd(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
430:     VecGetArray(mpimat_1->lvec, &cpcol_1_state);
431:     /* get 'cpcol_2_state' */
432:     VecScatterBegin(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
433:       VecScatterEnd(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
434:     VecGetArray(mpimat_2->lvec, &cpcol_2_state);
435:     /* get 'cpcol_2_par_orig' */
436:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
437:       PetscScalar v = (PetscScalar)lid_parent_gid[kk];
438:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
439:     }
440:     VecAssemblyBegin(tempVec);
441:     VecAssemblyEnd(tempVec);
442:     VecDuplicate(mpimat_2->lvec, &ghost_par_orig2);
443:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
444:       VecScatterEnd(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
445:     VecGetArray(ghost_par_orig2, &cpcol_2_par_orig);

447:     VecDestroy(&tempVec);
448:   } /* ismpi */

450:   /* doit */
451:   for (lid=0; lid<nloc; lid++) {
452:     NState state = lid_state[lid];
453:     if (IS_SELECTED(state)) {
454:       /* steal locals */
455:       ii  = matA_1->i; n = ii[lid+1] - ii[lid];
456:       idx = matA_1->j + ii[lid];
457:       for (j=0; j<n; j++) {
458:         PetscInt lidj   = idx[j], sgid;
459:         NState   statej = lid_state[lidj];
460:         if (statej==DELETED && (sgid=(PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid+my0) { /* steal local */
461:           lid_parent_gid[lidj] = (PetscScalar)(lid+my0); /* send this if sgid is not local */
462:           if (sgid >= my0 && sgid < Iend) {       /* I'm stealing this local from a local sgid */
463:             PetscInt   hav=0,slid=sgid-my0,gidj=lidj+my0;
464:             PetscCDPos pos,last=NULL;
465:             /* looking for local from local so id_llist_2 works */
466:             PetscCDGetHeadPos(aggs_2,slid,&pos);
467:             while (pos) {
468:               PetscInt gid;
469:               PetscLLNGetID(pos, &gid);
470:               if (gid == gidj) {
471:                 if (!last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null");
472:                 PetscCDRemoveNextNode(aggs_2, slid, last);
473:                 PetscCDAppendNode(aggs_2, lid, pos);
474:                 hav  = 1;
475:                 break;
476:               } else last = pos;

478:               PetscCDGetNextPos(aggs_2,slid,&pos);
479:             }
480:             if (hav!=1) {
481:               if (!hav) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
482:               SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav);
483:             }
484:           } else {            /* I'm stealing this local, owned by a ghost */
485:             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 0.0' if the matrix is structurally symmetric.");
486:             PetscCDAppendID(aggs_2, lid, lidj+my0);
487:           }
488:         }
489:       } /* local neighbors */
490:     } else if (state == DELETED && lid_cprowID_1) {
491:       PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
492:       /* see if I have a selected ghost neighbor that will steal me */
493:       if ((ix=lid_cprowID_1[lid]) != -1) {
494:         ii  = matB_1->compressedrow.i; n = ii[ix+1] - ii[ix];
495:         idx = matB_1->j + ii[ix];
496:         for (j=0; j<n; j++) {
497:           PetscInt cpid   = idx[j];
498:           NState   statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
499:           if (IS_SELECTED(statej) && sgidold != (PetscInt)statej) { /* ghost will steal this, remove from my list */
500:             lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
501:             if (sgidold>=my0 && sgidold<Iend) { /* this was mine */
502:               PetscInt   hav=0,oldslidj=sgidold-my0;
503:               PetscCDPos pos,last=NULL;
504:               /* remove from 'oldslidj' list */
505:               PetscCDGetHeadPos(aggs_2,oldslidj,&pos);
506:               while (pos) {
507:                 PetscInt gid;
508:                 PetscLLNGetID(pos, &gid);
509:                 if (lid+my0 == gid) {
510:                   /* id_llist_2[lastid] = id_llist_2[flid];   /\* remove lid from oldslidj list *\/ */
511:                   if (!last) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"last cannot be null");
512:                   PetscCDRemoveNextNode(aggs_2, oldslidj, last);
513:                   /* ghost (PetscScalar)statej will add this later */
514:                   hav = 1;
515:                   break;
516:                 } else last = pos;

518:                 PetscCDGetNextPos(aggs_2,oldslidj,&pos);
519:               }
520:               if (hav!=1) {
521:                 if (hav==0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
522:                 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"found node %d times???",hav);
523:               }
524:             } else {
525:               /* ghosts remove this later */
526:             }
527:           }
528:         }
529:       }
530:     } /* selected/deleted */
531:   } /* node loop */

533:   if (isMPI) {
534:     PetscScalar   *cpcol_2_parent,*cpcol_2_gid;
535:     Vec           tempVec,ghostgids2,ghostparents2;
536:     PetscInt      cpid,nghost_2;
537:     GAMGHashTable gid_cpid;

539:     VecGetSize(mpimat_2->lvec, &nghost_2);
540:     MatGetVecs(Gmat_2, &tempVec, 0);

542:     /* get 'cpcol_2_parent' */
543:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
544:       VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES);
545:     }
546:     VecAssemblyBegin(tempVec);
547:     VecAssemblyEnd(tempVec);
548:     VecDuplicate(mpimat_2->lvec, &ghostparents2);
549:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
550:     VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
551:     VecGetArray(ghostparents2, &cpcol_2_parent);

553:     /* get 'cpcol_2_gid' */
554:     for (kk=0,j=my0; kk<nloc; kk++,j++) {
555:       PetscScalar v = (PetscScalar)j;
556:       VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES);
557:     }
558:     VecAssemblyBegin(tempVec);
559:     VecAssemblyEnd(tempVec);
560:     VecDuplicate(mpimat_2->lvec, &ghostgids2);
561:     VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
562:       VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
563:     VecGetArray(ghostgids2, &cpcol_2_gid);

565:     VecDestroy(&tempVec);

567:     /* look for deleted ghosts and add to table */
568:     GAMGTableCreate(2*nghost_2+1, &gid_cpid);
569:     for (cpid = 0; cpid < nghost_2; cpid++) {
570:       NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
571:       if (state==DELETED) {
572:         PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
573:         PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
574:         if (sgid_old == -1 && sgid_new != -1) {
575:           PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
576:           GAMGTableAdd(&gid_cpid, gid, cpid);
577:         }
578:       }
579:     }

581:     /* look for deleted ghosts and see if they moved - remove it */
582:     for (lid=0; lid<nloc; lid++) {
583:       NState state = lid_state[lid];
584:       if (IS_SELECTED(state)) {
585:         PetscCDPos pos,last=NULL;
586:         /* look for deleted ghosts and see if they moved */
587:         PetscCDGetHeadPos(aggs_2,lid,&pos);
588:         while (pos) {
589:           PetscInt gid;
590:           PetscLLNGetID(pos, &gid);

592:           if (gid < my0 || gid >= Iend) {
593:             GAMGTableFind(&gid_cpid, gid, &cpid);
594:             if (cpid != -1) {
595:               /* a moved ghost - */
596:               /* id_llist_2[lastid] = id_llist_2[flid];    /\* remove 'flid' from list *\/ */
597:               PetscCDRemoveNextNode(aggs_2, lid, last);
598:             } else last = pos;
599:           } else last = pos;

601:           PetscCDGetNextPos(aggs_2,lid,&pos);
602:         } /* loop over list of deleted */
603:       } /* selected */
604:     }
605:     GAMGTableDestroy(&gid_cpid);

607:     /* look at ghosts, see if they changed - and it */
608:     for (cpid = 0; cpid < nghost_2; cpid++) {
609:       PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
610:       if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */
611:         PetscInt   gid     = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
612:         PetscInt   slid_new=sgid_new-my0,hav=0;
613:         PetscCDPos pos;
614:         /* search for this gid to see if I have it */
615:         PetscCDGetHeadPos(aggs_2,slid_new,&pos);
616:         while (pos) {
617:           PetscInt gidj;
618:           PetscLLNGetID(pos, &gidj);
619:           PetscCDGetNextPos(aggs_2,slid_new,&pos);

621:           if (gidj == gid) { hav = 1; break; }
622:         }
623:         if (hav != 1) {
624:           /* insert 'flidj' into head of llist */
625:           PetscCDAppendID(aggs_2, slid_new, gid);
626:         }
627:       }
628:     }

630:     VecRestoreArray(mpimat_1->lvec, &cpcol_1_state);
631:     VecRestoreArray(mpimat_2->lvec, &cpcol_2_state);
632:     VecRestoreArray(ghostparents2, &cpcol_2_parent);
633:     VecRestoreArray(ghostgids2, &cpcol_2_gid);
634:     PetscFree(lid_cprowID_1);
635:     VecDestroy(&ghostgids2);
636:     VecDestroy(&ghostparents2);
637:     VecDestroy(&ghost_par_orig2);
638:   }

640:   PetscFree(lid_parent_gid);
641:   PetscFree(lid_state);
642:   return(0);
643: }

645: /* -------------------------------------------------------------------------- */
646: /*
647:    PCSetData_AGG - called if data is not set with PCSetCoordinates.
648:       Looks in Mat for near null space.
649:       Does not work for Stokes

651:   Input Parameter:
652:    . pc -
653:    . a_A - matrix to get (near) null space out of.
654: */
657: PetscErrorCode PCSetData_AGG(PC pc, Mat a_A)
658: {
660:   PC_MG          *mg      = (PC_MG*)pc->data;
661:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
662:   MatNullSpace   mnull;

665:   MatGetNearNullSpace(a_A, &mnull);
666:   if (!mnull) {
667:     PetscInt bs,NN,MM;
668:     MatGetBlockSize(a_A, &bs);
669:     MatGetLocalSize(a_A, &MM, &NN);
670:     if (MM % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MM %D must be divisible by bs %D",MM,bs);
671:     PCSetCoordinates_AGG(pc, bs, MM/bs, NULL);
672:   } else {
673:     PetscReal *nullvec;
674:     PetscBool has_const;
675:     PetscInt  i,j,mlocal,nvec,bs;
676:     const Vec *vecs; const PetscScalar *v;
677:     MatGetLocalSize(a_A,&mlocal,NULL);
678:     MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
679:     pc_gamg->data_sz = (nvec+!!has_const)*mlocal;
680:     PetscMalloc((nvec+!!has_const)*mlocal*sizeof(*nullvec),&nullvec);
681:     if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0;
682:     for (i=0; i<nvec; i++) {
683:       VecGetArrayRead(vecs[i],&v);
684:       for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]);
685:       VecRestoreArrayRead(vecs[i],&v);
686:     }
687:     pc_gamg->data           = nullvec;
688:     pc_gamg->data_cell_cols = (nvec+!!has_const);

690:     MatGetBlockSize(a_A, &bs);

692:     pc_gamg->data_cell_rows = bs;
693:   }
694:   return(0);
695: }

697: /* -------------------------------------------------------------------------- */
698: /*
699:  formProl0

701:    Input Parameter:
702:    . agg_llists - list of arrays with aggregates
703:    . bs - block size
704:    . nSAvec - column bs of new P
705:    . my0crs - global index of start of locals
706:    . data_stride - bs*(nloc nodes + ghost nodes)
707:    . data_in[data_stride*nSAvec] - local data on fine grid
708:    . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
709:   Output Parameter:
710:    . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
711:    . a_Prol - prolongation operator
712: */
715: static PetscErrorCode formProl0(const PetscCoarsenData *agg_llists, /* list from selected vertices of aggregate unselected vertices */
716:                                 const PetscInt bs,          /* (row) block size */
717:                                 const PetscInt nSAvec,      /* column bs */
718:                                 const PetscInt my0crs,      /* global index of start of locals */
719:                                 const PetscInt data_stride, /* (nloc+nghost)*bs */
720:                                 PetscReal      data_in[],   /* [data_stride][nSAvec] */
721:                                 const PetscInt flid_fgid[], /* [data_stride/bs] */
722:                                 PetscReal **a_data_out,
723:                                 Mat a_Prol) /* prolongation operator (output)*/
724: {
726:   PetscInt       Istart,my0,Iend,nloc,clid,flid,aggID,kk,jj,ii,mm,ndone,nSelected,minsz,nghosts,out_data_stride;
727:   MPI_Comm       comm;
728:   PetscMPIInt    rank, size;
729:   PetscReal      *out_data;
730:   PetscCDPos     pos;
731:   GAMGHashTable  fgid_flid;

733: /* #define OUT_AGGS */
734: #if defined(OUT_AGGS)
735:   static PetscInt llev = 0; char fname[32]; FILE *file = NULL; PetscInt pM;
736: #endif

739:   PetscObjectGetComm((PetscObject)a_Prol,&comm);
740:   MPI_Comm_rank(comm,&rank);
741:   MPI_Comm_size(comm,&size);
742:   MatGetOwnershipRange(a_Prol, &Istart, &Iend);
743:   nloc = (Iend-Istart)/bs; my0 = Istart/bs;
744:   if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Iend %D - Istart %d must be divisible by bs %D",Iend,Istart,bs);
745:   Iend   /= bs;
746:   nghosts = data_stride/bs - nloc;

748:   GAMGTableCreate(2*nghosts+1, &fgid_flid);
749:   for (kk=0; kk<nghosts; kk++) {
750:     GAMGTableAdd(&fgid_flid, flid_fgid[nloc+kk], nloc+kk);
751:   }

753: #if defined(OUT_AGGS)
754:   sprintf(fname,"aggs_%d_%d.m",llev++,rank);
755:   if (llev==1) file = fopen(fname,"w");
756:   MatGetSize(a_Prol, &pM, &jj);
757: #endif

759:   /* count selected -- same as number of cols of P */
760:   for (nSelected=mm=0; mm<nloc; mm++) {
761:     PetscBool ise;
762:     PetscCDEmptyAt(agg_llists, mm, &ise);
763:     if (!ise) nSelected++;
764:   }
765:   MatGetOwnershipRangeColumn(a_Prol, &ii, &jj);
766:   if ((ii/nSAvec) != my0crs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D /nSAvec %D  != my0crs %D",ii,nSAvec,my0crs);
767:   if (nSelected != (jj-ii)/nSAvec) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"nSelected %D != (jj %D - ii %D)/nSAvec %D",nSelected,jj,ii,nSAvec);

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

772:   PetscMalloc(out_data_stride*nSAvec*sizeof(PetscReal), &out_data);
773:   for (ii=0;ii<out_data_stride*nSAvec;ii++) out_data[ii]=PETSC_MAX_REAL;
774:   *a_data_out = out_data; /* output - stride nSelected*nSAvec */

776:   /* find points and set prolongation */
777:   minsz = 100;
778:   ndone = 0;
779:   for (mm = clid = 0; mm < nloc; mm++) {
780:     PetscCDSizeAt(agg_llists, mm, &jj);
781:     if (jj > 0) {
782:       const PetscInt lid = mm, cgid = my0crs + clid;
783:       PetscInt       cids[100]; /* max bs */
784:       PetscBLASInt   asz  =jj,M=asz*bs,N=nSAvec,INFO;
785:       PetscBLASInt   Mdata=M+((N-M>0) ? N-M : 0),LDA=Mdata,LWORK=N*bs;
786:       PetscScalar    *qqc,*qqr,*TAU,*WORK;
787:       PetscInt       *fids;
788:       PetscReal      *data;
789:       /* count agg */
790:       if (asz<minsz) minsz = asz;

792:       /* get block */
793:       PetscMalloc((Mdata*N)*sizeof(PetscScalar), &qqc);
794:       PetscMalloc((M*N)*sizeof(PetscScalar), &qqr);
795:       PetscMalloc(N*sizeof(PetscScalar), &TAU);
796:       PetscMalloc(LWORK*sizeof(PetscScalar), &WORK);
797:       PetscMalloc(M*sizeof(PetscInt), &fids);

799:       aggID = 0;
800:       PetscCDGetHeadPos(agg_llists,lid,&pos);
801:       while (pos) {
802:         PetscInt gid1;
803:         PetscLLNGetID(pos, &gid1);
804:         PetscCDGetNextPos(agg_llists,lid,&pos);

806:         if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0;
807:         else {
808:           GAMGTableFind(&fgid_flid, gid1, &flid);
809:           if (flid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot find gid1 in table");
810:         }
811:         /* copy in B_i matrix - column oriented */
812:         data = &data_in[flid*bs];
813:         for (kk = ii = 0; ii < bs; ii++) {
814:           for (jj = 0; jj < N; jj++) {
815:             PetscReal d = data[jj*data_stride + ii];
816:             qqc[jj*Mdata + aggID*bs + ii] = d;
817:           }
818:         }
819: #if defined(OUT_AGGS)
820:         if (llev==1) {
821:           char     str[] = "plot(%e,%e,'r*'), hold on,\n", col[] = "rgbkmc", sim[] = "*os+h>d<vx^";
822:           PetscInt MM,pi,pj;
823:           str[12] = col[(clid+17*rank)%6]; str[13] = sim[(clid+17*rank)%11];
824:           MM      = (PetscInt)(PetscSqrtReal((PetscReal)pM));
825:           pj      = gid1/MM; pi = gid1%MM;
826:           fprintf(file,str,(double)pi,(double)pj);
827:           /* fprintf(file,str,data[2*data_stride+1],-data[2*data_stride]); */
828:         }
829: #endif
830:         /* set fine IDs */
831:         for (kk=0; kk<bs; kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk;

833:         aggID++;
834:       }

836:       /* pad with zeros */
837:       for (ii = asz*bs; ii < Mdata; ii++) {
838:         for (jj = 0; jj < N; jj++, kk++) {
839:           qqc[jj*Mdata + ii] = .0;
840:         }
841:       }

843:       ndone += aggID;
844:       /* QR */
845:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
846:       PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO));
847:       PetscFPTrapPop();
848:       if (INFO != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xGEQRF error");
849:       /* get R - column oriented - output B_{i+1} */
850:       {
851:         PetscReal *data = &out_data[clid*nSAvec];
852:         for (jj = 0; jj < nSAvec; jj++) {
853:           for (ii = 0; ii < nSAvec; ii++) {
854:             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);
855:            if (ii <= jj) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]);
856:            else data[jj*out_data_stride + ii] = 0.;
857:           }
858:         }
859:       }

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

865:       for (ii = 0; ii < M; ii++) {
866:         for (jj = 0; jj < N; jj++) {
867:           qqr[N*ii + jj] = qqc[jj*Mdata + ii];
868:         }
869:       }

871:       /* add diagonal block of P0 */
872:       for (kk=0; kk<N; kk++) {
873:         cids[kk] = N*cgid + kk; /* global col IDs in P0 */
874:       }
875:       MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);

877:       PetscFree(qqc);
878:       PetscFree(qqr);
879:       PetscFree(TAU);
880:       PetscFree(WORK);
881:       PetscFree(fids);
882:       clid++;
883:     } /* coarse agg */
884:   } /* for all fine nodes */
885:   MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
886:   MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);

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

893: #if defined(OUT_AGGS)
894:   if (llev==1) fclose(file);
895: #endif
896:   GAMGTableDestroy(&fgid_flid);
897:   return(0);
898: }

900: /* -------------------------------------------------------------------------- */
901: /*
902:    PCGAMGgraph_AGG

904:   Input Parameter:
905:    . pc - this
906:    . Amat - matrix on this fine level
907:   Output Parameter:
908:    . a_Gmat -
909: */
912: PetscErrorCode PCGAMGgraph_AGG(PC pc,const Mat Amat,Mat *a_Gmat)
913: {
914:   PetscErrorCode            ierr;
915:   PC_MG                     *mg          = (PC_MG*)pc->data;
916:   PC_GAMG                   *pc_gamg     = (PC_GAMG*)mg->innerctx;
917:   const PetscInt            verbose      = pc_gamg->verbose;
918:   const PetscReal           vfilter      = pc_gamg->threshold;
919:   PC_GAMG_AGG               *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
920:   PetscMPIInt               rank,size;
921:   Mat                       Gmat;
922:   MPI_Comm                  comm;
923:   PetscBool /* set,flg , */ symm;

926:   PetscObjectGetComm((PetscObject)Amat,&comm);
927: #if defined PETSC_USE_LOG
928:   PetscLogEventBegin(PC_GAMGGgraph_AGG,0,0,0,0);
929: #endif
930:   MPI_Comm_rank(comm, &rank);
931:   MPI_Comm_size(comm, &size);

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

936:   PCGAMGCreateGraph(Amat, &Gmat);
937:   PCGAMGFilterGraph(&Gmat, vfilter, symm, verbose);

939:   *a_Gmat = Gmat;

941: #if defined PETSC_USE_LOG
942:   PetscLogEventEnd(PC_GAMGGgraph_AGG,0,0,0,0);
943: #endif
944:   return(0);
945: }

947: /* -------------------------------------------------------------------------- */
948: /*
949:    PCGAMGCoarsen_AGG

951:   Input Parameter:
952:    . a_pc - this
953:   Input/Output Parameter:
954:    . a_Gmat1 - graph on this fine level - coarsening can change this (squares it)
955:   Output Parameter:
956:    . agg_lists - list of aggregates
957: */
960: PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc,Mat *a_Gmat1,PetscCoarsenData **agg_lists)
961: {
963:   PC_MG          *mg          = (PC_MG*)a_pc->data;
964:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
965:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
966:   Mat            mat,Gmat2, Gmat1 = *a_Gmat1;  /* squared graph */
967:   const PetscInt verbose = pc_gamg->verbose;
968:   IS             perm;
969:   PetscInt       Ii,nloc,bs,n,m;
970:   PetscInt       *permute;
971:   PetscBool      *bIndexSet;
972:   MatCoarsen     crs;
973:   MPI_Comm       comm;
974:   PetscMPIInt    rank,size;

977: #if defined PETSC_USE_LOG
978:   PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);
979: #endif
980:   PetscObjectGetComm((PetscObject)Gmat1,&comm);
981:   MPI_Comm_rank(comm, &rank);
982:   MPI_Comm_size(comm, &size);
983:   MatGetLocalSize(Gmat1, &n, &m);
984:   MatGetBlockSize(Gmat1, &bs);
985:   if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bs %D must be 1",bs);
986:   nloc = n/bs;

988:   if (pc_gamg_agg->square_graph) {
989:     if (verbose > 1) PetscPrintf(comm,"[%d]%s square graph\n",rank,__FUNCT__);
990:     /* MatMatTransposeMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2); */
991:     MatTransposeMatMult(Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2);
992:     if (verbose > 2) {
993:       PetscPrintf(comm,"[%d]%s square graph done\n",rank,__FUNCT__);
994:       /* check for symetry */
995:       if (verbose > 4) {

997:       }
998:     }
999:   } else Gmat2 = Gmat1;

1001:   /* get MIS aggs */
1002:   /* randomize */
1003:   PetscMalloc(nloc*sizeof(PetscInt), &permute);
1004:   PetscMalloc(nloc*sizeof(PetscBool), &bIndexSet);
1005:   for (Ii = 0; Ii < nloc; Ii++) {
1006:     bIndexSet[Ii] = PETSC_FALSE;
1007:     permute[Ii]   = Ii;
1008:   }
1009:   srand(1); /* make deterministic */
1010:   for (Ii = 0; Ii < nloc; Ii++) {
1011:     PetscInt iSwapIndex = rand()%nloc;
1012:     if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1013:       PetscInt iTemp = permute[iSwapIndex];
1014:       permute[iSwapIndex]   = permute[Ii];
1015:       permute[Ii]           = iTemp;
1016:       bIndexSet[iSwapIndex] = PETSC_TRUE;
1017:     }
1018:   }
1019:   PetscFree(bIndexSet);

1021:   if (verbose > 1) PetscPrintf(comm,"[%d]%s coarsen graph\n",rank,__FUNCT__);

1023:   ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);
1024: #if defined PETSC_GAMG_USE_LOG
1025:   PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);
1026: #endif
1027:   MatCoarsenCreate(comm, &crs);
1028:   /* MatCoarsenSetType(crs, MATCOARSENMIS); */
1029:   MatCoarsenSetFromOptions(crs);
1030:   MatCoarsenSetGreedyOrdering(crs, perm);
1031:   MatCoarsenSetAdjacency(crs, Gmat2);
1032:   MatCoarsenSetVerbose(crs, pc_gamg->verbose);
1033:   MatCoarsenSetStrictAggs(crs, PETSC_TRUE);
1034:   MatCoarsenApply(crs);
1035:   MatCoarsenGetData(crs, agg_lists); /* output */
1036:   MatCoarsenDestroy(&crs);

1038:   ISDestroy(&perm);
1039:   PetscFree(permute);
1040: #if defined PETSC_GAMG_USE_LOG
1041:   PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);
1042: #endif
1043:   if (verbose > 2) PetscPrintf(comm,"[%d]%s coarsen graph done\n",rank,__FUNCT__);

1045:   /* smooth aggs */
1046:   if (Gmat2 != Gmat1) {
1047:     const PetscCoarsenData *llist = *agg_lists;
1048:     smoothAggs(Gmat2, Gmat1, *agg_lists);
1049:     MatDestroy(&Gmat1);
1050:     *a_Gmat1 = Gmat2; /* output */
1051:     PetscCDGetMat(llist, &mat);
1052:     if (mat) SETERRQ(comm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????");
1053:   } else {
1054:     const PetscCoarsenData *llist = *agg_lists;
1055:     /* see if we have a matrix that takes precedence (returned from MatCoarsenAppply) */
1056:     PetscCDGetMat(llist, &mat);
1057:     if (mat) {
1058:       MatDestroy(&Gmat1);
1059:       *a_Gmat1 = mat; /* output */
1060:     }
1061:   }
1062: #if defined PETSC_USE_LOG
1063:   PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);
1064: #endif
1065:   if (verbose > 2) PetscPrintf(comm,"[%d]%s PCGAMGCoarsen_AGG done\n",rank,__FUNCT__);
1066:   return(0);
1067: }

1069: /* -------------------------------------------------------------------------- */
1070: /*
1071:  PCGAMGProlongator_AGG

1073:  Input Parameter:
1074:  . pc - this
1075:  . Amat - matrix on this fine level
1076:  . Graph - used to get ghost data for nodes in
1077:  . agg_lists - list of aggregates
1078:  Output Parameter:
1079:  . a_P_out - prolongation operator to the next level
1080:  */
1083: PetscErrorCode PCGAMGProlongator_AGG(PC pc,const Mat Amat,const Mat Gmat,PetscCoarsenData *agg_lists,Mat *a_P_out)
1084: {
1085:   PC_MG          *mg       = (PC_MG*)pc->data;
1086:   PC_GAMG        *pc_gamg  = (PC_GAMG*)mg->innerctx;
1087:   const PetscInt verbose   = pc_gamg->verbose;
1088:   const PetscInt data_cols = pc_gamg->data_cell_cols;
1090:   PetscInt       Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs;
1091:   Mat            Prol;
1092:   PetscMPIInt    rank, size;
1093:   MPI_Comm       comm;
1094:   const PetscInt col_bs = data_cols;
1095:   PetscReal      *data_w_ghost;
1096:   PetscInt       myCrs0, nbnodes=0, *flid_fgid;

1099:   PetscObjectGetComm((PetscObject)Amat,&comm);
1100: #if defined PETSC_USE_LOG
1101:   PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);
1102: #endif
1103:   MPI_Comm_rank(comm, &rank);
1104:   MPI_Comm_size(comm, &size);
1105:   MatGetOwnershipRange(Amat, &Istart, &Iend);
1106:   MatGetBlockSize(Amat, &bs);
1107:   nloc = (Iend-Istart)/bs; my0 = Istart/bs;
1108:   if ((Iend-Istart) % bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"(Iend %D - Istart %D) not divisible by bs %D",Iend,Istart,bs);

1110:   /* get 'nLocalSelected' */
1111:   for (ii=0, nLocalSelected = 0; ii < nloc; ii++) {
1112:     PetscBool ise;
1113:     /* filter out singletons 0 or 1? */
1114:     PetscCDEmptyAt(agg_lists, ii, &ise);
1115:     if (!ise) nLocalSelected++;
1116:   }

1118:   /* create prolongator, create P matrix */
1119:   MatCreate(comm, &Prol);
1120:   MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);
1121:   MatSetBlockSizes(Prol, bs, col_bs);
1122:   MatSetType(Prol, MATAIJ);
1123:   MatSeqAIJSetPreallocation(Prol, data_cols, NULL);
1124:   MatMPIAIJSetPreallocation(Prol,data_cols, NULL,data_cols, NULL);
1125:   /* nloc*bs, nLocalSelected*col_bs, */
1126:   /* PETSC_DETERMINE, PETSC_DETERMINE, */
1127:   /* data_cols, NULL, data_cols, NULL, */
1128:   /* &Prol); */

1130:   /* can get all points "removed" */
1131:    MatGetSize(Prol, &kk, &ii);
1132:   if (ii==0) {
1133:     if (verbose > 0) PetscPrintf(comm,"[%d]%s no selected points on coarse grid\n",rank,__FUNCT__);
1134:     MatDestroy(&Prol);
1135:     *a_P_out = NULL;  /* out */
1136:     return(0);
1137:   }
1138:   if (verbose > 0) PetscPrintf(comm,"\t\t[%d]%s New grid %d nodes\n",rank,__FUNCT__,ii/col_bs);
1139:   MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk);

1141:   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);
1142:   myCrs0 = myCrs0/col_bs;
1143:   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);

1145:   /* create global vector of data in 'data_w_ghost' */
1146: #if defined PETSC_GAMG_USE_LOG
1147:   PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);
1148: #endif
1149:   if (size > 1) { /*  */
1150:     PetscReal *tmp_gdata,*tmp_ldata,*tp2;
1151:     PetscMalloc(nloc*sizeof(PetscReal), &tmp_ldata);
1152:     for (jj = 0; jj < data_cols; jj++) {
1153:       for (kk = 0; kk < bs; kk++) {
1154:         PetscInt        ii,stride;
1155:         const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk;
1156:         for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp;

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

1160:         if (jj==0 && kk==0) { /* now I know how many todal nodes - allocate */
1161:           PetscMalloc(stride*bs*data_cols*sizeof(PetscReal), &data_w_ghost);
1162:           nbnodes = bs*stride;
1163:         }
1164:         tp2 = data_w_ghost + jj*bs*stride + kk;
1165:         for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii];
1166:         PetscFree(tmp_gdata);
1167:       }
1168:     }
1169:     PetscFree(tmp_ldata);
1170:   } else {
1171:     nbnodes      = bs*nloc;
1172:     data_w_ghost = (PetscReal*)pc_gamg->data;
1173:   }

1175:   /* get P0 */
1176:   if (size > 1) {
1177:     PetscReal *fid_glid_loc,*fiddata;
1178:     PetscInt  stride;

1180:     PetscMalloc(nloc*sizeof(PetscReal), &fid_glid_loc);
1181:     for (kk=0; kk<nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk);
1182:     PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata);
1183:     PetscMalloc(stride*sizeof(PetscInt), &flid_fgid);
1184:     for (kk=0; kk<stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1185:     PetscFree(fiddata);

1187:     if (stride != nbnodes/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"stride %D != nbnodes %D/bs %D",stride,nbnodes,bs);
1188:     PetscFree(fid_glid_loc);
1189:   } else {
1190:     PetscMalloc(nloc*sizeof(PetscInt), &flid_fgid);
1191:     for (kk=0; kk<nloc; kk++) flid_fgid[kk] = my0 + kk;
1192:   }
1193: #if defined PETSC_GAMG_USE_LOG
1194:   PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);
1195:   PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);
1196: #endif
1197:   {
1198:     PetscReal *data_out = NULL;
1199:     formProl0(agg_lists, bs, data_cols, myCrs0, nbnodes,data_w_ghost, flid_fgid, &data_out, Prol);
1200:     PetscFree(pc_gamg->data);

1202:     pc_gamg->data           = data_out;
1203:     pc_gamg->data_cell_rows = data_cols;
1204:     pc_gamg->data_sz        = data_cols*data_cols*nLocalSelected;
1205:   }
1206: #if defined PETSC_GAMG_USE_LOG
1207:   PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);
1208: #endif
1209:   if (size > 1) PetscFree(data_w_ghost);
1210:   PetscFree(flid_fgid);

1212:   *a_P_out = Prol;  /* out */
1213: #if defined PETSC_USE_LOG
1214:   PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1215: #endif
1216:   return(0);
1217: }

1219: /* -------------------------------------------------------------------------- */
1220: /*
1221:    PCGAMGOptprol_AGG

1223:   Input Parameter:
1224:    . pc - this
1225:    . Amat - matrix on this fine level
1226:  In/Output Parameter:
1227:    . a_P_out - prolongation operator to the next level
1228: */
1231: PetscErrorCode PCGAMGOptprol_AGG(PC pc,const Mat Amat,Mat *a_P)
1232: {
1234:   PC_MG          *mg          = (PC_MG*)pc->data;
1235:   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
1236:   const PetscInt verbose      = pc_gamg->verbose;
1237:   PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
1238:   PetscInt       jj;
1239:   PetscMPIInt    rank,size;
1240:   Mat            Prol  = *a_P;
1241:   MPI_Comm       comm;

1244:   PetscObjectGetComm((PetscObject)Amat,&comm);
1245: #if defined PETSC_USE_LOG
1246:   PetscLogEventBegin(PC_GAMGOptprol_AGG,0,0,0,0);
1247: #endif
1248:   MPI_Comm_rank(comm, &rank);
1249:   MPI_Comm_size(comm, &size);

1251:   /* smooth P0 */
1252:   for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) {
1253:     Mat       tMat;
1254:     Vec       diag;
1255:     PetscReal alpha, emax, emin;
1256: #if defined PETSC_GAMG_USE_LOG
1257:     PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);
1258: #endif
1259:     if (jj == 0) {
1260:       KSP eksp;
1261:       Vec bb, xx;
1262:       PC  epc;
1263:       MatGetVecs(Amat, &bb, 0);
1264:       MatGetVecs(Amat, &xx, 0);
1265:       {
1266:         PetscRandom rctx;
1267:         PetscRandomCreate(comm,&rctx);
1268:         PetscRandomSetFromOptions(rctx);
1269:         VecSetRandom(bb,rctx);
1270:         PetscRandomDestroy(&rctx);
1271:       }

1273:       /* zeroing out BC rows -- needed for crazy matrices */
1274:       {
1275:         PetscInt    Istart,Iend,ncols,jj,Ii;
1276:         PetscScalar zero = 0.0;
1277:         MatGetOwnershipRange(Amat, &Istart, &Iend);
1278:         for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) {
1279:           MatGetRow(Amat,Ii,&ncols,0,0);
1280:           if (ncols <= 1) {
1281:             VecSetValues(bb, 1, &Ii, &zero, INSERT_VALUES);
1282:           }
1283:           MatRestoreRow(Amat,Ii,&ncols,0,0);
1284:         }
1285:         VecAssemblyBegin(bb);
1286:         VecAssemblyEnd(bb);
1287:       }

1289:       KSPCreate(comm,&eksp);
1290:       KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
1291:       KSPSetNormType(eksp, KSP_NORM_NONE);
1292:       KSPSetOptionsPrefix(eksp,((PetscObject)pc)->prefix);
1293:       KSPAppendOptionsPrefix(eksp, "gamg_est_");
1294:       KSPSetFromOptions(eksp);

1296:       KSPSetInitialGuessNonzero(eksp, PETSC_FALSE);
1297:       KSPSetOperators(eksp, Amat, Amat, SAME_NONZERO_PATTERN);
1298:       KSPSetComputeSingularValues(eksp,PETSC_TRUE);

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

1303:       /* solve - keep stuff out of logging */
1304:       PetscLogEventDeactivate(KSP_Solve);
1305:       PetscLogEventDeactivate(PC_Apply);
1306:       KSPSolve(eksp, bb, xx);
1307:       PetscLogEventActivate(KSP_Solve);
1308:       PetscLogEventActivate(PC_Apply);

1310:       KSPComputeExtremeSingularValues(eksp, &emax, &emin);
1311:       if (verbose > 0) {
1312:         PetscPrintf(comm,"\t\t\t%s smooth P0: max eigen=%e min=%e PC=%s\n",
1313:                     __FUNCT__,emax,emin,PCJACOBI);
1314:       }
1315:       VecDestroy(&xx);
1316:       VecDestroy(&bb);
1317:       KSPDestroy(&eksp);

1319:       if (pc_gamg->emax_id == -1) {
1320:         PetscObjectComposedDataRegister(&pc_gamg->emax_id);
1321:         if (pc_gamg->emax_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"pc_gamg->emax_id == -1");
1322:       }
1323:       PetscObjectComposedDataSetScalar((PetscObject)Amat, pc_gamg->emax_id, emax);
1324:     }

1326:     /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1327:     MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat);
1328:     MatGetVecs(Amat, &diag, 0);
1329:     MatGetDiagonal(Amat, diag); /* effectively PCJACOBI */
1330:     VecReciprocal(diag);
1331:     MatDiagonalScale(tMat, diag, 0);
1332:     VecDestroy(&diag);
1333:     alpha = -1.4/emax;
1334:     MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN);
1335:     MatDestroy(&Prol);
1336:     Prol  = tMat;
1337: #if defined PETSC_GAMG_USE_LOG
1338:     PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);
1339: #endif
1340:   }
1341: #if defined PETSC_USE_LOG
1342:   PetscLogEventEnd(PC_GAMGOptprol_AGG,0,0,0,0);
1343: #endif
1344:   *a_P = Prol;
1345:   return(0);
1346: }

1348: /* -------------------------------------------------------------------------- */
1349: /*
1350:    PCGAMGKKTProl_AGG

1352:   Input Parameter:
1353:    . pc - this
1354:    . Prol11 - matrix on this fine level
1355:    . A21 - matrix on this fine level
1356:  In/Output Parameter:
1357:    . a_P22 - prolongation operator to the next level
1358: */
1361: PetscErrorCode PCGAMGKKTProl_AGG(PC pc,const Mat Prol11,const Mat A21,Mat *a_P22)
1362: {
1363:   PetscErrorCode   ierr;
1364:   PC_MG            *mg      = (PC_MG*)pc->data;
1365:   PC_GAMG          *pc_gamg = (PC_GAMG*)mg->innerctx;
1366:   const PetscInt   verbose  = pc_gamg->verbose;
1367:   /* PC_GAMG_AGG    *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;  */
1368:   PetscMPIInt      rank,size;
1369:   Mat              Prol22,Tmat,Gmat;
1370:   MPI_Comm         comm;
1371:   PetscCoarsenData *agg_lists;

1374:   PetscObjectGetComm((PetscObject)pc,&comm);
1375: #if defined PETSC_USE_LOG
1376:   PetscLogEventBegin(PC_GAMGKKTProl_AGG,0,0,0,0);
1377: #endif
1378:   MPI_Comm_rank(comm, &rank);
1379:   MPI_Comm_size(comm, &size);

1381:   /* form C graph */
1382:   MatMatMult(A21, Prol11, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Tmat);
1383:   MatMatTransposeMult(Tmat, Tmat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat);
1384:   MatDestroy(&Tmat);
1385:   PCGAMGFilterGraph(&Gmat, 0.0, PETSC_FALSE, verbose);

1387:   /* coarsen constraints */
1388:   {
1389:     MatCoarsen crs;
1390:     MatCoarsenCreate(comm, &crs);
1391:     MatCoarsenSetType(crs, MATCOARSENMIS);
1392:     MatCoarsenSetAdjacency(crs, Gmat);
1393:     MatCoarsenSetVerbose(crs, verbose);
1394:     MatCoarsenSetStrictAggs(crs, PETSC_TRUE);
1395:     MatCoarsenApply(crs);
1396:     MatCoarsenGetData(crs, &agg_lists);
1397:     MatCoarsenDestroy(&crs);
1398:   }

1400:   /* form simple prolongation 'Prol22' */
1401:   {
1402:     PetscInt    ii,mm,clid,my0,nloc,nLocalSelected;
1403:     PetscScalar val = 1.0;
1404:     /* get 'nLocalSelected' */
1405:     MatGetLocalSize(Gmat, &nloc, &ii);
1406:     for (ii=0, nLocalSelected = 0; ii < nloc; ii++) {
1407:       PetscBool ise;
1408:       /* filter out singletons 0 or 1? */
1409:       PetscCDEmptyAt(agg_lists, ii, &ise);
1410:       if (!ise) nLocalSelected++;
1411:     }

1413:     MatCreate(comm,&Prol22);
1414:     MatSetSizes(Prol22,nloc, nLocalSelected,PETSC_DETERMINE, PETSC_DETERMINE);
1415:     MatSetType(Prol22, MATAIJ);
1416:     MatSeqAIJSetPreallocation(Prol22,1,NULL);
1417:     MatMPIAIJSetPreallocation(Prol22,1,NULL,1,NULL);
1418:     /* MatCreateAIJ(comm, */
1419:     /*                      nloc, nLocalSelected, */
1420:     /*                      PETSC_DETERMINE, PETSC_DETERMINE, */
1421:     /*                      1, NULL, 1, NULL, */
1422:     /*                      &Prol22); */

1424:     MatGetOwnershipRange(Prol22, &my0, &ii);
1425:     nloc = ii - my0;

1427:     /* make aggregates */
1428:     for (mm = clid = 0; mm < nloc; mm++) {
1429:       PetscCDSizeAt(agg_lists, mm, &ii);
1430:       if (ii > 0) {
1431:         PetscInt   asz=ii,cgid=my0+clid,rids[1000];
1432:         PetscCDPos pos;
1433:         if (asz>1000) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Very large aggregate: %d",asz);
1434:         ii   = 0;
1435:         PetscCDGetHeadPos(agg_lists,mm,&pos);
1436:         while (pos) {
1437:           PetscInt gid1;
1438:           PetscLLNGetID(pos, &gid1);
1439:           PetscCDGetNextPos(agg_lists,mm,&pos);

1441:           rids[ii++] = gid1;
1442:         }
1443:         if (ii != asz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D != asz %D",ii,asz);
1444:         /* add diagonal block of P0 */
1445:         MatSetValues(Prol22,asz,rids,1,&cgid,&val,INSERT_VALUES);

1447:         clid++;
1448:       } /* coarse agg */
1449:     } /* for all fine nodes */
1450:     MatAssemblyBegin(Prol22,MAT_FINAL_ASSEMBLY);
1451:     MatAssemblyEnd(Prol22,MAT_FINAL_ASSEMBLY);
1452:   }

1454:   /* clean up */
1455:   MatDestroy(&Gmat);
1456:   PetscCDDestroy(agg_lists);
1457: #if defined PETSC_USE_LOG
1458:   PetscLogEventEnd(PC_GAMGKKTProl_AGG,0,0,0,0);
1459: #endif
1460:   *a_P22 = Prol22;
1461:   return(0);
1462: }

1464: /* -------------------------------------------------------------------------- */
1465: /*
1466:    PCCreateGAMG_AGG

1468:   Input Parameter:
1469:    . pc -
1470: */
1473: PetscErrorCode  PCCreateGAMG_AGG(PC pc)
1474: {
1476:   PC_MG          *mg      = (PC_MG*)pc->data;
1477:   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1478:   PC_GAMG_AGG    *pc_gamg_agg;

1481:   /* create sub context for SA */
1482:   PetscNewLog(pc, PC_GAMG_AGG, &pc_gamg_agg);
1483:   pc_gamg->subctx = pc_gamg_agg;

1485:   pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1486:   pc_gamg->ops->destroy        = PCDestroy_GAMG_AGG;
1487:   /* reset does not do anything; setup not virtual */

1489:   /* set internal function pointers */
1490:   pc_gamg->ops->graph       = PCGAMGgraph_AGG;
1491:   pc_gamg->ops->coarsen     = PCGAMGCoarsen_AGG;
1492:   pc_gamg->ops->prolongator = PCGAMGProlongator_AGG;
1493:   pc_gamg->ops->optprol     = PCGAMGOptprol_AGG;
1494:   pc_gamg->ops->formkktprol = PCGAMGKKTProl_AGG;

1496:   pc_gamg->ops->createdefaultdata = PCSetData_AGG;

1498:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_AGG);
1499:   return(0);
1500: }