Actual source code: agg.c
petsc-3.3-p7 2013-05-11
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 <assert.h>
9: #include <petscblaslapack.h>
11: typedef struct {
12: PetscInt nsmooths;
13: PetscBool sym_graph;
14: PetscBool square_graph;
15: }PC_GAMG_AGG;
19: /*@
20: PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical)
22: Not Collective on PC
24: Input Parameters:
25: . pc - the preconditioner context
27: Options Database Key:
28: . -pc_gamg_agg_nsmooths
30: Level: intermediate
32: Concepts: Aggregation AMG preconditioner
34: .seealso: ()
35: @*/
36: PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n)
37: {
42: PetscTryMethod(pc,"PCGAMGSetNSmooths_C",(PC,PetscInt),(pc,n));
43: return(0);
44: }
46: EXTERN_C_BEGIN
49: PetscErrorCode PCGAMGSetNSmooths_GAMG(PC pc, PetscInt n)
50: {
51: PC_MG *mg = (PC_MG*)pc->data;
52: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
53: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
56: pc_gamg_agg->nsmooths = n;
57: return(0);
58: }
59: EXTERN_C_END
63: /*@
64: PCGAMGSetSymGraph -
66: Not Collective on PC
68: Input Parameters:
69: . pc - the preconditioner context
71: Options Database Key:
72: . -pc_gamg_sym_graph
74: Level: intermediate
76: Concepts: Aggregation AMG preconditioner
78: .seealso: ()
79: @*/
80: PetscErrorCode PCGAMGSetSymGraph(PC pc, PetscBool n)
81: {
86: PetscTryMethod(pc,"PCGAMGSetSymGraph_C",(PC,PetscBool),(pc,n));
87: return(0);
88: }
90: EXTERN_C_BEGIN
93: PetscErrorCode PCGAMGSetSymGraph_GAMG(PC pc, PetscBool n)
94: {
95: PC_MG *mg = (PC_MG*)pc->data;
96: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
97: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
100: pc_gamg_agg->sym_graph = n;
101: return(0);
102: }
103: EXTERN_C_END
107: /*@
108: PCGAMGSetSquareGraph -
110: Not Collective on PC
112: Input Parameters:
113: . pc - the preconditioner context
115: Options Database Key:
116: . -pc_gamg_square_graph
118: Level: intermediate
120: Concepts: Aggregation AMG preconditioner
122: .seealso: ()
123: @*/
124: PetscErrorCode PCGAMGSetSquareGraph(PC pc, PetscBool n)
125: {
130: PetscTryMethod(pc,"PCGAMGSetSquareGraph_C",(PC,PetscBool),(pc,n));
131: return(0);
132: }
134: EXTERN_C_BEGIN
137: PetscErrorCode PCGAMGSetSquareGraph_GAMG(PC pc, PetscBool n)
138: {
139: PC_MG *mg = (PC_MG*)pc->data;
140: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
141: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
144: pc_gamg_agg->square_graph = n;
145: return(0);
146: }
147: EXTERN_C_END
149: /* -------------------------------------------------------------------------- */
150: /*
151: PCSetFromOptions_GAMG_AGG
153: Input Parameter:
154: . pc -
155: */
158: PetscErrorCode PCSetFromOptions_GAMG_AGG( PC pc )
159: {
160: PetscErrorCode ierr;
161: PC_MG *mg = (PC_MG*)pc->data;
162: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
163: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
164: PetscBool flag;
165:
167: /* call base class */
168: PCSetFromOptions_GAMG( pc );
170: PetscOptionsHead("GAMG-AGG options");
171: {
172: /* -pc_gamg_agg_nsmooths */
173: pc_gamg_agg->nsmooths = 0;
174: PetscOptionsInt("-pc_gamg_agg_nsmooths",
175: "smoothing steps for smoothed aggregation, usually 1 (0)",
176: "PCGAMGSetNSmooths",
177: pc_gamg_agg->nsmooths,
178: &pc_gamg_agg->nsmooths,
179: &flag);
180:
182: /* -pc_gamg_sym_graph */
183: pc_gamg_agg->sym_graph = PETSC_FALSE;
184: PetscOptionsBool("-pc_gamg_sym_graph",
185: "Set for asymmetric matrices",
186: "PCGAMGSetSymGraph",
187: pc_gamg_agg->sym_graph,
188: &pc_gamg_agg->sym_graph,
189: &flag);
190:
192: /* -pc_gamg_square_graph */
193: pc_gamg_agg->square_graph = PETSC_TRUE;
194: PetscOptionsBool("-pc_gamg_square_graph",
195: "For faster coarsening and lower coarse grid complexity",
196: "PCGAMGSetSquareGraph",
197: pc_gamg_agg->square_graph,
198: &pc_gamg_agg->square_graph,
199: &flag);
200:
201: }
202: PetscOptionsTail();
203:
204: return(0);
205: }
207: /* -------------------------------------------------------------------------- */
208: /*
209: PCDestroy_AGG
211: Input Parameter:
212: . pc -
213: */
216: PetscErrorCode PCDestroy_AGG( PC pc )
217: {
218: PetscErrorCode ierr;
219: PC_MG *mg = (PC_MG*)pc->data;
220: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
221: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
222:
224: if( pc_gamg_agg ) {
225: PetscFree(pc_gamg_agg);
226: pc_gamg_agg = 0;
227: }
229: /* call base class */
230: PCDestroy_GAMG( pc );
231:
232: return(0);
233: }
235: /* -------------------------------------------------------------------------- */
236: /*
237: PCSetCoordinates_AGG
238: - collective
240: Input Parameter:
241: . pc - the preconditioner context
242: . ndm - dimesion of data (used for dof/vertex for Stokes)
243: . a_nloc - number of vertices local
244: . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...}
245: */
246: EXTERN_C_BEGIN
249: PetscErrorCode PCSetCoordinates_AGG( PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords )
250: {
251: PC_MG *mg = (PC_MG*)pc->data;
252: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
253: PetscErrorCode ierr;
254: PetscInt arrsz,kk,ii,jj,nloc,ndatarows,ndf;
255: Mat mat = pc->pmat;
256: /* MPI_Comm wcomm = ((PetscObject)pc)->comm; */
261: nloc = a_nloc;
263: /* SA: null space vectors */
264: MatGetBlockSize( mat, &ndf ); CHKERRQ( ierr ); /* this does not work for Stokes */
265: if( coords && ndf==1 ) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */
266: else if( coords ) {
267: if(ndm > ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"degrees of motion %d > block size %d",ndm,ndf);
268: pc_gamg->data_cell_cols = (ndm==2 ? 3 : 6); /* displacement elasticity */
269: if(ndm != ndf) {
270: if(pc_gamg->data_cell_cols != ndf) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Don't know how to create null space for ndm=%d, ndf=%d. Use MatSetNearNullSpace.",ndm,ndf);
271: }
272: }
273: else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */
274: pc_gamg->data_cell_rows = ndatarows = ndf; assert(pc_gamg->data_cell_cols>0);
275: arrsz = nloc*pc_gamg->data_cell_rows*pc_gamg->data_cell_cols;
277: /* create data - syntactic sugar that should be refactored at some point */
278: if (pc_gamg->data==0 || (pc_gamg->data_sz != arrsz)) {
279: PetscFree( pc_gamg->data );
280: PetscMalloc((arrsz+1)*sizeof(PetscReal), &pc_gamg->data );
281: }
282: /* copy data in - column oriented */
283: for(kk=0;kk<nloc;kk++){
284: const PetscInt M = nloc*pc_gamg->data_cell_rows; /* stride into data */
285: PetscReal *data = &pc_gamg->data[kk*ndatarows]; /* start of cell */
286: if( pc_gamg->data_cell_cols==1 ) *data = 1.0;
287: else {
288: /* translational modes */
289: for(ii=0;ii<ndatarows;ii++)
290: for(jj=0;jj<ndatarows;jj++)
291: if(ii==jj)data[ii*M + jj] = 1.0;
292: else data[ii*M + jj] = 0.0;
294: /* rotational modes */
295: if( coords ) {
296: if( ndm == 2 ){
297: data += 2*M;
298: data[0] = -coords[2*kk+1];
299: data[1] = coords[2*kk];
300: }
301: else {
302: data += 3*M;
303: data[0] = 0.0; data[M+0] = coords[3*kk+2]; data[2*M+0] = -coords[3*kk+1];
304: data[1] = -coords[3*kk+2]; data[M+1] = 0.0; data[2*M+1] = coords[3*kk];
305: data[2] = coords[3*kk+1]; data[M+2] = -coords[3*kk]; data[2*M+2] = 0.0;
306: }
307: }
308: }
309: }
311: pc_gamg->data_sz = arrsz;
313: return(0);
314: }
315: EXTERN_C_END
317: typedef PetscInt NState;
318: static const NState NOT_DONE=-2;
319: static const NState DELETED=-1;
320: static const NState REMOVED=-3;
321: #define IS_SELECTED(s) (s!=DELETED && s!=NOT_DONE && s!=REMOVED)
323: /* -------------------------------------------------------------------------- */
324: /*
325: smoothAggs - greedy grab of with G1 (unsquared graph) -- AIJ specific
326: - AGG-MG specific: clears singletons out of 'selected_2'
328: Input Parameter:
329: . Gmat_2 - glabal matrix of graph (data not defined)
330: . Gmat_1 - base graph to grab with
331: Input/Output Parameter:
332: . aggs_2 - linked list of aggs with gids )
333: */
336: static PetscErrorCode smoothAggs( const Mat Gmat_2, /* base (squared) graph */
337: const Mat Gmat_1, /* base graph */
338: /* const IS selected_2, [nselected local] selected vertices */
339: PetscCoarsenData *aggs_2 /* [nselected local] global ID of aggregate */
340: )
341: {
343: PetscBool isMPI;
344: Mat_SeqAIJ *matA_1, *matB_1=0, *matA_2, *matB_2=0;
345: MPI_Comm wcomm = ((PetscObject)Gmat_2)->comm;
346: PetscMPIInt mype,npe;
347: PetscInt lid,*ii,*idx,ix,Iend,my0,kk,n,j;
348: Mat_MPIAIJ *mpimat_2 = 0, *mpimat_1=0;
349: const PetscInt nloc = Gmat_2->rmap->n;
350: PetscScalar *cpcol_1_state,*cpcol_2_state,*cpcol_2_par_orig,*lid_parent_gid;
351: PetscInt *lid_cprowID_1;
352: NState *lid_state;
353: Vec ghost_par_orig2;
356: MPI_Comm_rank( wcomm, &mype );
357: MPI_Comm_size( wcomm, &npe );
358: MatGetOwnershipRange(Gmat_1,&my0,&Iend);
360: if( PETSC_FALSE ) {
361: PetscViewer viewer; char fname[32]; static int llev=0;
362: sprintf(fname,"Gmat2_%d.m",llev++);
363: PetscViewerASCIIOpen(wcomm,fname,&viewer);
364: PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB);
365: MatView(Gmat_2, viewer );
366: PetscViewerDestroy( &viewer );
367: }
369: /* get submatrices */
370: PetscObjectTypeCompare( (PetscObject)Gmat_1, MATMPIAIJ, &isMPI );
371: if(isMPI) {
372: /* grab matrix objects */
373: mpimat_2 = (Mat_MPIAIJ*)Gmat_2->data;
374: mpimat_1 = (Mat_MPIAIJ*)Gmat_1->data;
375: matA_1 = (Mat_SeqAIJ*)mpimat_1->A->data;
376: matB_1 = (Mat_SeqAIJ*)mpimat_1->B->data;
377: matA_2 = (Mat_SeqAIJ*)mpimat_2->A->data;
378: matB_2 = (Mat_SeqAIJ*)mpimat_2->B->data;
380: /* force compressed row storage for B matrix in AuxMat */
381: matB_1->compressedrow.check = PETSC_TRUE;
382: MatCheckCompressedRow(mpimat_1->B,&matB_1->compressedrow,matB_1->i,Gmat_1->rmap->n,-1.0);
383:
385: PetscMalloc( nloc*sizeof(PetscInt), &lid_cprowID_1 );
386: for( lid = 0 ; lid < nloc ; lid++ ) lid_cprowID_1[lid] = -1;
387: for (ix=0; ix<matB_1->compressedrow.nrows; ix++) {
388: PetscInt lid = matB_1->compressedrow.rindex[ix];
389: lid_cprowID_1[lid] = ix;
390: }
391: }
392: else {
393: matA_1 = (Mat_SeqAIJ*)Gmat_1->data;
394: matA_2 = (Mat_SeqAIJ*)Gmat_2->data;
395: lid_cprowID_1 = PETSC_NULL;
396: }
397: assert( matA_1 && !matA_1->compressedrow.use );
398: assert( matB_1==0 || matB_1->compressedrow.use );
399: assert( matA_2 && !matA_2->compressedrow.use );
400: assert( matB_2==0 || matB_2->compressedrow.use );
402: /* get state of locals and selected gid for deleted */
403: PetscMalloc( nloc*sizeof(NState), &lid_state );
404: PetscMalloc( nloc*sizeof(PetscScalar), &lid_parent_gid );
405: for( lid = 0 ; lid < nloc ; lid++ ) {
406: lid_parent_gid[lid] = -1.0;
407: lid_state[lid] = DELETED;
408: }
409:
410: /* set lid_state */
411: for( lid = 0 ; lid < nloc ; lid++ ) {
412: PetscCDPos pos;
413: PetscCDGetHeadPos(aggs_2,lid,&pos);
414: if( pos ) {
415: PetscInt gid1;
416: PetscLLNGetID( pos, &gid1 ); assert(gid1==lid+my0);
417: lid_state[lid] = gid1;
418: }
419: }
421: /* map local to selected local, DELETED means a ghost owns it */
422: for(lid=kk=0;lid<nloc;lid++){
423: NState state = lid_state[lid];
424: if( IS_SELECTED(state) ){
425: PetscCDPos pos;
426: PetscCDGetHeadPos(aggs_2,lid,&pos);
427: while(pos){
428: PetscInt gid1;
429: PetscLLNGetID( pos, &gid1 );
430: PetscCDGetNextPos(aggs_2,lid,&pos);
431:
432: if( gid1 >= my0 && gid1 < Iend ){
433: lid_parent_gid[gid1-my0] = (PetscScalar)(lid + my0);
434: }
435: }
436: }
437: }
438: /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */
439: if (isMPI) {
440: Vec tempVec;
441: /* get 'cpcol_1_state' */
442: MatGetVecs( Gmat_1, &tempVec, 0 );
443: for(kk=0,j=my0;kk<nloc;kk++,j++){
444: PetscScalar v = (PetscScalar)lid_state[kk];
445: VecSetValues( tempVec, 1, &j, &v, INSERT_VALUES );
446: }
447: VecAssemblyBegin( tempVec );
448: VecAssemblyEnd( tempVec );
449: VecScatterBegin(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
450:
451: VecScatterEnd(mpimat_1->Mvctx,tempVec, mpimat_1->lvec,INSERT_VALUES,SCATTER_FORWARD);
452:
453: VecGetArray( mpimat_1->lvec, &cpcol_1_state );
454: /* get 'cpcol_2_state' */
455: VecScatterBegin(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
456:
457: VecScatterEnd(mpimat_2->Mvctx,tempVec, mpimat_2->lvec,INSERT_VALUES,SCATTER_FORWARD);
458:
459: VecGetArray( mpimat_2->lvec, &cpcol_2_state );
460: /* get 'cpcol_2_par_orig' */
461: for(kk=0,j=my0;kk<nloc;kk++,j++){
462: PetscScalar v = (PetscScalar)lid_parent_gid[kk];
463: VecSetValues( tempVec, 1, &j, &v, INSERT_VALUES );
464: }
465: VecAssemblyBegin( tempVec );
466: VecAssemblyEnd( tempVec );
467: VecDuplicate( mpimat_2->lvec, &ghost_par_orig2 );
468: VecScatterBegin(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
469:
470: VecScatterEnd(mpimat_2->Mvctx,tempVec, ghost_par_orig2,INSERT_VALUES,SCATTER_FORWARD);
471:
472: VecGetArray( ghost_par_orig2, &cpcol_2_par_orig );
474: VecDestroy( &tempVec );
475: } /* ismpi */
477: /* doit */
478: for(lid=0;lid<nloc;lid++){
479: NState state = lid_state[lid];
480: if( IS_SELECTED(state) ) {
481: /* steal locals */
482: ii = matA_1->i; n = ii[lid+1] - ii[lid];
483: idx = matA_1->j + ii[lid];
484: for (j=0; j<n; j++) {
485: PetscInt lidj = idx[j], sgid;
486: NState statej = lid_state[lidj];
487: if (statej==DELETED && (sgid=(PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid+my0) { /* steal local */
488: lid_parent_gid[lidj] = (PetscScalar)(lid+my0); /* send this if sgid is not local */
489: if( sgid >= my0 && sgid < Iend ){ /* I'm stealing this local from a local sgid */
490: PetscInt hav=0,slid=sgid-my0,gidj=lidj+my0;
491: PetscCDPos pos,last=PETSC_NULL;
492: /* looking for local from local so id_llist_2 works */
493: PetscCDGetHeadPos(aggs_2,slid,&pos);
494: while(pos){
495: PetscInt gid;
496: PetscLLNGetID( pos, &gid );
497: if( gid == gidj ) {
498: assert(last);
499: PetscCDRemoveNextNode( aggs_2, slid, last );
500: PetscCDAppendNode( aggs_2, lid, pos );
501: hav = 1;
502: break;
503: }
504: else last = pos;
506: PetscCDGetNextPos(aggs_2,slid,&pos);
507: }
508: if(hav!=1){
509: if(hav==0)SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
510: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"found node %d times???",hav);
511: }
512: }
513: else{ /* I'm stealing this local, owned by a ghost */
514: assert(sgid==-1);
515: PetscCDAppendID( aggs_2, lid, lidj+my0 );
516: }
517: }
518: } /* local neighbors */
519: }
520: else if( state == DELETED && lid_cprowID_1 ) {
521: PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]);
522: /* see if I have a selected ghost neighbor that will steal me */
523: if( (ix=lid_cprowID_1[lid]) != -1 ){
524: ii = matB_1->compressedrow.i; n = ii[ix+1] - ii[ix];
525: idx = matB_1->j + ii[ix];
526: for( j=0 ; j<n ; j++ ) {
527: PetscInt cpid = idx[j];
528: NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]);
529: if( IS_SELECTED(statej) && sgidold != (PetscInt)statej ) { /* ghost will steal this, remove from my list */
530: lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */
531: if( sgidold>=my0 && sgidold<Iend ) { /* this was mine */
532: PetscInt hav=0,oldslidj=sgidold-my0;
533: PetscCDPos pos,last=PETSC_NULL;
534: /* remove from 'oldslidj' list */
535: PetscCDGetHeadPos(aggs_2,oldslidj,&pos);
536: while( pos ) {
537: PetscInt gid;
538: PetscLLNGetID( pos, &gid );
539: if( lid+my0 == gid ) {
540: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */
541: assert(last);
542: PetscCDRemoveNextNode( aggs_2, oldslidj, last );
543: /* ghost (PetscScalar)statej will add this later */
544: hav = 1;
545: break;
546: }
547: else last = pos;
549: PetscCDGetNextPos(aggs_2,oldslidj,&pos);
550: }
551: if(hav!=1){
552: if(hav==0)SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"failed to find adj in 'selected' lists - structurally unsymmetric matrix");
553: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"found node %d times???",hav);
554: }
555: }
556: else {
557: /* ghosts remove this later */
558: }
559: }
560: }
561: }
562: } /* selected/deleted */
563: } /* node loop */
565: if( isMPI ) {
566: PetscScalar *cpcol_2_parent,*cpcol_2_gid;
567: Vec tempVec,ghostgids2,ghostparents2;
568: PetscInt cpid,nghost_2;
569: GAMGHashTable gid_cpid;
571: VecGetSize( mpimat_2->lvec, &nghost_2 );
572: MatGetVecs( Gmat_2, &tempVec, 0 );
574: /* get 'cpcol_2_parent' */
575: for(kk=0,j=my0;kk<nloc;kk++,j++){
576: VecSetValues( tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES );
577: }
578: VecAssemblyBegin( tempVec );
579: VecAssemblyEnd( tempVec );
580: VecDuplicate( mpimat_2->lvec, &ghostparents2 );
581: VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
582:
583: VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostparents2,INSERT_VALUES,SCATTER_FORWARD);
584:
585: VecGetArray( ghostparents2, &cpcol_2_parent );
587: /* get 'cpcol_2_gid' */
588: for(kk=0,j=my0;kk<nloc;kk++,j++){
589: PetscScalar v = (PetscScalar)j;
590: VecSetValues( tempVec, 1, &j, &v, INSERT_VALUES );
591: }
592: VecAssemblyBegin( tempVec );
593: VecAssemblyEnd( tempVec );
594: VecDuplicate( mpimat_2->lvec, &ghostgids2 );
595: VecScatterBegin(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
596:
597: VecScatterEnd(mpimat_2->Mvctx,tempVec, ghostgids2,INSERT_VALUES,SCATTER_FORWARD);
598:
599: VecGetArray( ghostgids2, &cpcol_2_gid );
601: VecDestroy( &tempVec );
603: /* look for deleted ghosts and add to table */
604: GAMGTableCreate( 2*nghost_2, &gid_cpid );
605: for( cpid = 0 ; cpid < nghost_2 ; cpid++ ) {
606: NState state = (NState)PetscRealPart(cpcol_2_state[cpid]);
607: if( state==DELETED ) {
608: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
609: PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]);
610: if( sgid_old == -1 && sgid_new != -1 ) {
611: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
612: GAMGTableAdd( &gid_cpid, gid, cpid );
613: }
614: }
615: }
617: /* look for deleted ghosts and see if they moved - remove it */
618: for(lid=0;lid<nloc;lid++){
619: NState state = lid_state[lid];
620: if( IS_SELECTED(state) ){
621: PetscCDPos pos,last=PETSC_NULL;
622: /* look for deleted ghosts and see if they moved */
623: PetscCDGetHeadPos(aggs_2,lid,&pos);
624: while(pos){
625: PetscInt gid;
626: PetscLLNGetID( pos, &gid );
628: if( gid < my0 || gid >= Iend ) {
629: GAMGTableFind( &gid_cpid, gid, &cpid );
630: if( cpid != -1 ) {
631: /* a moved ghost - */
632: /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */
633: PetscCDRemoveNextNode( aggs_2, lid, last );
634: }
635: else last = pos;
636: }
637: else last = pos;
639: PetscCDGetNextPos(aggs_2,lid,&pos);
640: } /* loop over list of deleted */
641: } /* selected */
642: }
643: GAMGTableDestroy( &gid_cpid );
645: /* look at ghosts, see if they changed - and it */
646: for( cpid = 0 ; cpid < nghost_2 ; cpid++ ){
647: PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]);
648: if( sgid_new >= my0 && sgid_new < Iend ) { /* this is mine */
649: PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]);
650: PetscInt slid_new=sgid_new-my0,hav=0;
651: PetscCDPos pos;
652: /* search for this gid to see if I have it */
653: PetscCDGetHeadPos(aggs_2,slid_new,&pos);
654: while(pos){
655: PetscInt gidj;
656: PetscLLNGetID( pos, &gidj );
657: PetscCDGetNextPos(aggs_2,slid_new,&pos);
658:
659: if( gidj == gid ) { hav = 1; break; }
660: }
661: if( hav != 1 ){
662: /* insert 'flidj' into head of llist */
663: PetscCDAppendID( aggs_2, slid_new, gid );
664: }
665: }
666: }
668: VecRestoreArray( mpimat_1->lvec, &cpcol_1_state );
669: VecRestoreArray( mpimat_2->lvec, &cpcol_2_state );
670: VecRestoreArray( ghostparents2, &cpcol_2_parent );
671: VecRestoreArray( ghostgids2, &cpcol_2_gid );
672: PetscFree( lid_cprowID_1 );
673: VecDestroy( &ghostgids2 );
674: VecDestroy( &ghostparents2 );
675: VecDestroy( &ghost_par_orig2 );
676: }
678: PetscFree( lid_parent_gid );
679: PetscFree( lid_state );
681: return(0);
682: }
684: /* -------------------------------------------------------------------------- */
685: /*
686: PCSetData_AGG - called if data is not set with PCSetCoordinates.
687: Looks in Mat for near null space.
688: Does not work for Stokes
690: Input Parameter:
691: . pc -
692: . a_A - matrix to get (near) null space out of.
693: */
696: PetscErrorCode PCSetData_AGG( PC pc, Mat a_A )
697: {
698: PetscErrorCode ierr;
699: PC_MG *mg = (PC_MG*)pc->data;
700: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
701: MatNullSpace mnull;
704: MatGetNearNullSpace( a_A, &mnull );
705: if( !mnull ) {
706: PetscInt bs,NN,MM;
707: MatGetBlockSize( a_A, &bs ); CHKERRQ( ierr );
708: MatGetLocalSize( a_A, &MM, &NN ); CHKERRQ( ierr ); assert(MM%bs==0);
709: PetscPrintf(((PetscObject)a_A)->comm,"[%d]%s bs=%d MM=%d\n",0,__FUNCT__,bs,MM);
710: PCSetCoordinates_AGG( pc, bs, MM/bs, PETSC_NULL );
711: }
712: else {
713: PetscReal *nullvec;
714: PetscBool has_const;
715: PetscInt i,j,mlocal,nvec,bs;
716: const Vec *vecs; const PetscScalar *v;
717: MatGetLocalSize(a_A,&mlocal,PETSC_NULL);
718: MatNullSpaceGetVecs( mnull, &has_const, &nvec, &vecs );
719: PetscMalloc((nvec+!!has_const)*mlocal*sizeof *nullvec,&nullvec);
720: if (has_const) for (i=0; i<mlocal; i++) nullvec[i] = 1.0;
721: for (i=0; i<nvec; i++) {
722: VecGetArrayRead(vecs[i],&v);
723: for (j=0; j<mlocal; j++) nullvec[(i+!!has_const)*mlocal + j] = PetscRealPart(v[j]);
724: VecRestoreArrayRead(vecs[i],&v);
725: }
726: pc_gamg->data = nullvec;
727: pc_gamg->data_cell_cols = (nvec+!!has_const);
728: MatGetBlockSize( a_A, &bs ); CHKERRQ( ierr ); /* this does not work for Stokes */
729: pc_gamg->data_cell_rows = bs;
730: }
731: return(0);
732: }
734: /* -------------------------------------------------------------------------- */
735: /*
736: formProl0
738: Input Parameter:
739: . agg_llists - list of arrays with aggregates
740: . bs - block size
741: . nSAvec - column bs of new P
742: . my0crs - global index of start of locals
743: . data_stride - bs*(nloc nodes + ghost nodes)
744: . data_in[data_stride*nSAvec] - local data on fine grid
745: . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist'
746: Output Parameter:
747: . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data
748: . a_Prol - prolongation operator
749: */
752: static PetscErrorCode formProl0(const PetscCoarsenData *agg_llists,/* list from selected vertices of aggregate unselected vertices */
753: const PetscInt bs, /* (row) block size */
754: const PetscInt nSAvec, /* column bs */
755: const PetscInt my0crs, /* global index of start of locals */
756: const PetscInt data_stride, /* (nloc+nghost)*bs */
757: PetscReal data_in[], /* [data_stride][nSAvec] */
758: const PetscInt flid_fgid[], /* [data_stride/bs] */
759: PetscReal **a_data_out,
760: Mat a_Prol /* prolongation operator (output)*/
761: )
762: {
764: PetscInt Istart,my0,Iend,nloc,clid,flid,aggID,kk,jj,ii,mm,ndone,nSelected,minsz,nghosts,out_data_stride;
765: MPI_Comm wcomm = ((PetscObject)a_Prol)->comm;
766: PetscMPIInt mype, npe;
767: PetscReal *out_data;
768: PetscCDPos pos;
769: GAMGHashTable fgid_flid;
771: /* #define OUT_AGGS */
772: #ifdef OUT_AGGS
773: static PetscInt llev = 0; char fname[32]; FILE *file = PETSC_NULL; PetscInt pM;
774: #endif
777: MPI_Comm_rank(wcomm,&mype);
778: MPI_Comm_size(wcomm,&npe);
779: MatGetOwnershipRange( a_Prol, &Istart, &Iend );
780: nloc = (Iend-Istart)/bs; my0 = Istart/bs; assert((Iend-Istart)%bs==0);
781: Iend /= bs;
782: nghosts = data_stride/bs - nloc;
784: GAMGTableCreate( 2*nghosts, &fgid_flid );
785: for(kk=0;kk<nghosts;kk++) {
786: GAMGTableAdd( &fgid_flid, flid_fgid[nloc+kk], nloc+kk );
787: }
789: #ifdef OUT_AGGS
790: sprintf(fname,"aggs_%d_%d.m",llev++,mype);
791: if(llev==1) {
792: file = fopen(fname,"w");
793: }
794: MatGetSize( a_Prol, &pM, &jj );
795: #endif
797: /* count selected -- same as number of cols of P */
798: for(nSelected=mm=0;mm<nloc;mm++) {
799: PetscBool ise;
800: PetscCDEmptyAt( agg_llists, mm, &ise );
801: if( !ise ) nSelected++;
802: }
803: MatGetOwnershipRangeColumn( a_Prol, &ii, &jj );
804: assert((ii/nSAvec)==my0crs); assert(nSelected==(jj-ii)/nSAvec);
806: /* aloc space for coarse point data (output) */
807: out_data_stride = nSelected*nSAvec;
808: PetscMalloc( out_data_stride*nSAvec*sizeof(PetscReal), &out_data );
809: for(ii=0;ii<out_data_stride*nSAvec;ii++) {
810: out_data[ii]=1.e300;
811: }
812: *a_data_out = out_data; /* output - stride nSelected*nSAvec */
814: /* find points and set prolongation */
815: minsz = 100;
816: ndone = 0;
817: for( mm = clid = 0 ; mm < nloc ; mm++ ){
818: PetscCDSizeAt( agg_llists, mm, &jj );
819: if( jj > 0 ) {
820: const PetscInt lid = mm, cgid = my0crs + clid;
821: PetscInt cids[100]; /* max bs */
822: PetscBLASInt asz=jj,M=asz*bs,N=nSAvec,INFO;
823: PetscBLASInt Mdata=M+((N-M>0)?N-M:0),LDA=Mdata,LWORK=N*bs;
824: PetscScalar *qqc,*qqr,*TAU,*WORK;
825: PetscInt *fids;
826: PetscReal *data;
827: /* count agg */
828: if( asz<minsz ) minsz = asz;
830: /* get block */
831: PetscMalloc( (Mdata*N)*sizeof(PetscScalar), &qqc );
832: PetscMalloc( (M*N)*sizeof(PetscScalar), &qqr );
833: PetscMalloc( N*sizeof(PetscScalar), &TAU );
834: PetscMalloc( LWORK*sizeof(PetscScalar), &WORK );
835: PetscMalloc( M*sizeof(PetscInt), &fids );
837: aggID = 0;
838: PetscCDGetHeadPos(agg_llists,lid,&pos);
839: while(pos){
840: PetscInt gid1;
841: PetscLLNGetID( pos, &gid1 );
842: PetscCDGetNextPos(agg_llists,lid,&pos);
844: if( gid1 >= my0 && gid1 < Iend ) flid = gid1 - my0;
845: else {
846: GAMGTableFind( &fgid_flid, gid1, &flid );
847: assert(flid>=0);
848: }
849: /* copy in B_i matrix - column oriented */
850: data = &data_in[flid*bs];
851: for( kk = ii = 0; ii < bs ; ii++ ) {
852: for( jj = 0; jj < N ; jj++ ) {
853: PetscReal d = data[jj*data_stride + ii];
854: qqc[jj*Mdata + aggID*bs + ii] = d;
855: }
856: }
857: #ifdef OUT_AGGS
858: if(llev==1) {
859: char str[] = "plot(%e,%e,'r*'), hold on,\n", col[] = "rgbkmc", sim[] = "*os+h>d<vx^";
860: PetscInt MM,pi,pj;
861: str[12] = col[(clid+17*mype)%6]; str[13] = sim[(clid+17*mype)%11];
862: MM = (PetscInt)(PetscSqrtReal((PetscReal)pM));
863: pj = gid1/MM; pi = gid1%MM;
864: fprintf(file,str,(double)pi,(double)pj);
865: /* fprintf(file,str,data[2*data_stride+1],-data[2*data_stride]); */
866: }
867: #endif
868: /* set fine IDs */
869: for(kk=0;kk<bs;kk++) fids[aggID*bs + kk] = flid_fgid[flid]*bs + kk;
870:
871: aggID++;
872: }
874: /* pad with zeros */
875: for( ii = asz*bs; ii < Mdata ; ii++ ) {
876: for( jj = 0; jj < N ; jj++, kk++ ) {
877: qqc[jj*Mdata + ii] = .0;
878: }
879: }
881: ndone += aggID;
882: /* QR */
883: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
884: LAPACKgeqrf_( &Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO );
885: PetscFPTrapPop();
886: if( INFO != 0 ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRS error");
887: /* get R - column oriented - output B_{i+1} */
888: {
889: PetscReal *data = &out_data[clid*nSAvec];
890: for( jj = 0; jj < nSAvec ; jj++ ) {
891: for( ii = 0; ii < nSAvec ; ii++ ) {
892: assert(data[jj*out_data_stride + ii] == 1.e300);
893: if( ii <= jj ) data[jj*out_data_stride + ii] = PetscRealPart(qqc[jj*Mdata + ii]);
894: else data[jj*out_data_stride + ii] = 0.;
895: }
896: }
897: }
899: /* get Q - row oriented */
900: LAPACKungqr_( &Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO );
901: if( INFO != 0 ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR error arg %d",-INFO);
903: for( ii = 0 ; ii < M ; ii++ ){
904: for( jj = 0 ; jj < N ; jj++ ) {
905: qqr[N*ii + jj] = qqc[jj*Mdata + ii];
906: }
907: }
909: /* add diagonal block of P0 */
910: for(kk=0;kk<N;kk++) {
911: cids[kk] = N*cgid + kk; /* global col IDs in P0 */
912: }
913: MatSetValues(a_Prol,M,fids,N,cids,qqr,INSERT_VALUES);
915: PetscFree( qqc );
916: PetscFree( qqr );
917: PetscFree( TAU );
918: PetscFree( WORK );
919: PetscFree( fids );
920: clid++;
921: } /* coarse agg */
922: } /* for all fine nodes */
923: MatAssemblyBegin(a_Prol,MAT_FINAL_ASSEMBLY);
924: MatAssemblyEnd(a_Prol,MAT_FINAL_ASSEMBLY);
926: /* MPI_Allreduce( &ndone, &ii, 1, MPIU_INT, MPIU_SUM, wcomm ); */
927: /* MatGetSize( a_Prol, &kk, &jj ); */
928: /* MPI_Allreduce( &minsz, &jj, 1, MPIU_INT, MPIU_MIN, wcomm ); */
929: /* PetscPrintf(wcomm," **** [%d]%s %d total done, %d nodes (%d local done), min agg. size = %d\n",mype,__FUNCT__,ii,kk/bs,ndone,jj); */
931: #ifdef OUT_AGGS
932: if(llev==1) fclose(file);
933: #endif
934: GAMGTableDestroy( &fgid_flid );
936: return(0);
937: }
939: /* -------------------------------------------------------------------------- */
940: /*
941: PCGAMGgraph_AGG
943: Input Parameter:
944: . pc - this
945: . Amat - matrix on this fine level
946: Output Parameter:
947: . a_Gmat -
948: */
951: PetscErrorCode PCGAMGgraph_AGG( PC pc,
952: const Mat Amat,
953: Mat *a_Gmat
954: )
955: {
957: PC_MG *mg = (PC_MG*)pc->data;
958: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
959: const PetscInt verbose = pc_gamg->verbose;
960: const PetscReal vfilter = pc_gamg->threshold;
961: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
962: PetscMPIInt mype,npe;
963: Mat Gmat;
964: MPI_Comm wcomm = ((PetscObject)Amat)->comm;
965: PetscBool set,flg,symm;
968: #if defined PETSC_USE_LOG
969: PetscLogEventBegin(PC_GAMGGgraph_AGG,0,0,0,0);
970: #endif
971: MPI_Comm_rank( wcomm, &mype);
972: MPI_Comm_size( wcomm, &npe);
974: MatIsSymmetricKnown(Amat, &set, &flg);
975: symm = (PetscBool)(pc_gamg_agg->sym_graph || !(set && flg));
977: PCGAMGCreateGraph( Amat, &Gmat ); CHKERRQ( ierr );
978: PCGAMGFilterGraph( &Gmat, vfilter, symm, verbose ); CHKERRQ( ierr );
979:
980: *a_Gmat = Gmat;
982: #if defined PETSC_USE_LOG
983: PetscLogEventEnd(PC_GAMGGgraph_AGG,0,0,0,0);
984: #endif
985: return(0);
986: }
988: /* -------------------------------------------------------------------------- */
989: /*
990: PCGAMGCoarsen_AGG
992: Input Parameter:
993: . a_pc - this
994: Input/Output Parameter:
995: . a_Gmat1 - graph on this fine level - coarsening can change this (squares it)
996: Output Parameter:
997: . agg_lists - list of aggregates
998: */
1001: PetscErrorCode PCGAMGCoarsen_AGG( PC a_pc,
1002: Mat *a_Gmat1,
1003: PetscCoarsenData **agg_lists
1004: )
1005: {
1007: PC_MG *mg = (PC_MG*)a_pc->data;
1008: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1009: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
1010: Mat mat,Gmat2, Gmat1 = *a_Gmat1; /* squared graph */
1011: IS perm;
1012: PetscInt Ii,nloc,bs,n,m;
1013: PetscInt *permute;
1014: PetscBool *bIndexSet;
1015: MatCoarsen crs;
1016: MPI_Comm wcomm = ((PetscObject)Gmat1)->comm;
1017: PetscMPIInt mype,npe;
1020: #if defined PETSC_USE_LOG
1021: PetscLogEventBegin(PC_GAMGCoarsen_AGG,0,0,0,0);
1022: #endif
1023: MPI_Comm_rank( wcomm, &mype);
1024: MPI_Comm_size( wcomm, &npe);
1025: MatGetLocalSize( Gmat1, &n, &m );
1026: MatGetBlockSize( Gmat1, &bs ); assert(bs==1);
1027: nloc = n/bs;
1029: if( pc_gamg_agg->square_graph ) {
1030: MatTransposeMatMult( Gmat1, Gmat1, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat2 );
1031:
1032: }
1033: else Gmat2 = Gmat1;
1035: /* get MIS aggs */
1036: /* randomize */
1037: PetscMalloc( nloc*sizeof(PetscInt), &permute );
1038: PetscMalloc( nloc*sizeof(PetscBool), &bIndexSet );
1039: for ( Ii = 0; Ii < nloc ; Ii++ ){
1040: bIndexSet[Ii] = PETSC_FALSE;
1041: permute[Ii] = Ii;
1042: }
1043: srand(1); /* make deterministic */
1044: for ( Ii = 0; Ii < nloc ; Ii++ ) {
1045: PetscInt iSwapIndex = rand()%nloc;
1046: if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) {
1047: PetscInt iTemp = permute[iSwapIndex];
1048: permute[iSwapIndex] = permute[Ii];
1049: permute[Ii] = iTemp;
1050: bIndexSet[iSwapIndex] = PETSC_TRUE;
1051: }
1052: }
1053: PetscFree( bIndexSet );
1054:
1055: ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm);
1056:
1057: #if defined PETSC_GAMG_USE_LOG
1058: PetscLogEventBegin(petsc_gamg_setup_events[SET4],0,0,0,0);
1059: #endif
1060: MatCoarsenCreate( wcomm, &crs );
1061: /* MatCoarsenSetType( crs, MATCOARSENMIS ); */
1062: MatCoarsenSetFromOptions( crs );
1063: MatCoarsenSetGreedyOrdering( crs, perm );
1064: MatCoarsenSetAdjacency( crs, Gmat2 );
1065: MatCoarsenSetVerbose( crs, pc_gamg->verbose );
1066: MatCoarsenSetStrictAggs( crs, PETSC_TRUE );
1067: MatCoarsenApply( crs );
1068: MatCoarsenGetData( crs, agg_lists ); /* output */
1069: MatCoarsenDestroy( &crs );
1071: ISDestroy( &perm );
1072: PetscFree( permute );
1073: #if defined PETSC_GAMG_USE_LOG
1074: PetscLogEventEnd(petsc_gamg_setup_events[SET4],0,0,0,0);
1075: #endif
1076: /* smooth aggs */
1077: if( Gmat2 != Gmat1 ) {
1078: const PetscCoarsenData *llist = *agg_lists;
1079: smoothAggs( Gmat2, Gmat1, *agg_lists );
1080: MatDestroy( &Gmat1 );
1081: *a_Gmat1 = Gmat2; /* output */
1082: PetscCDGetMat( llist, &mat );
1083: if(mat) SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "Auxilary matrix with squared graph????");
1084: }
1085: else {
1086: const PetscCoarsenData *llist = *agg_lists;
1087: /* see if we have a matrix that takes pecedence (returned from MatCoarsenAppply) */
1088: PetscCDGetMat( llist, &mat );
1089: if( mat ) {
1090: MatDestroy( &Gmat1 );
1091: *a_Gmat1 = mat; /* output */
1092: }
1093: }
1094: #if defined PETSC_USE_LOG
1095: PetscLogEventEnd(PC_GAMGCoarsen_AGG,0,0,0,0);
1096: #endif
1097: return(0);
1098: }
1100: /* -------------------------------------------------------------------------- */
1101: /*
1102: PCGAMGProlongator_AGG
1103:
1104: Input Parameter:
1105: . pc - this
1106: . Amat - matrix on this fine level
1107: . Graph - used to get ghost data for nodes in
1108: . agg_lists - list of aggregates
1109: Output Parameter:
1110: . a_P_out - prolongation operator to the next level
1111: */
1114: PetscErrorCode PCGAMGProlongator_AGG( PC pc,
1115: const Mat Amat,
1116: const Mat Gmat,
1117: PetscCoarsenData *agg_lists,
1118: Mat *a_P_out
1119: )
1120: {
1121: PC_MG *mg = (PC_MG*)pc->data;
1122: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1123: const PetscInt verbose = pc_gamg->verbose;
1124: const PetscInt data_cols = pc_gamg->data_cell_cols;
1126: PetscInt Istart,Iend,nloc,ii,jj,kk,my0,nLocalSelected,bs;
1127: Mat Prol;
1128: PetscMPIInt mype, npe;
1129: MPI_Comm wcomm = ((PetscObject)Amat)->comm;
1130: const PetscInt col_bs = data_cols;
1131: PetscReal *data_w_ghost;
1132: PetscInt myCrs0, nbnodes=0, *flid_fgid;
1135: #if defined PETSC_USE_LOG
1136: PetscLogEventBegin(PC_GAMGProlongator_AGG,0,0,0,0);
1137: #endif
1138: MPI_Comm_rank( wcomm, &mype);
1139: MPI_Comm_size( wcomm, &npe);
1140: MatGetOwnershipRange( Amat, &Istart, &Iend );
1141: MatGetBlockSize( Amat, &bs ); CHKERRQ( ierr );
1142: nloc = (Iend-Istart)/bs; my0 = Istart/bs; assert((Iend-Istart)%bs==0);
1144: /* get 'nLocalSelected' */
1145: for( ii=0, nLocalSelected = 0 ; ii < nloc ; ii++ ){
1146: PetscBool ise;
1147: /* filter out singletons 0 or 1? */
1148: PetscCDEmptyAt( agg_lists, ii, &ise );
1149: if( !ise ) nLocalSelected++;
1150: }
1151:
1152: /* create prolongator, create P matrix */
1153: MatCreate( wcomm, &Prol );
1154: MatSetSizes(Prol,nloc*bs,nLocalSelected*col_bs,PETSC_DETERMINE,PETSC_DETERMINE);
1155:
1156: MatSetBlockSizes( Prol, bs, col_bs );
1157: MatSetType( Prol, MATAIJ );
1158: MatSeqAIJSetPreallocation( Prol, data_cols, PETSC_NULL);
1159: MatMPIAIJSetPreallocation(Prol,data_cols, PETSC_NULL,data_cols, PETSC_NULL);
1160: /* nloc*bs, nLocalSelected*col_bs, */
1161: /* PETSC_DETERMINE, PETSC_DETERMINE, */
1162: /* data_cols, PETSC_NULL, data_cols, PETSC_NULL, */
1163: /* &Prol ); */
1165: /* can get all points "removed" */
1166: MatGetSize( Prol, &kk, &ii );
1167: if( ii==0 ) {
1168: if( verbose ) {
1169: PetscPrintf(wcomm,"[%d]%s no selected points on coarse grid\n",mype,__FUNCT__);
1170: }
1171: MatDestroy( &Prol );
1172: *a_P_out = PETSC_NULL; /* out */
1173: return(0);
1174: }
1175: if( verbose ) {
1176: PetscPrintf(wcomm,"\t\t[%d]%s New grid %d nodes\n",mype,__FUNCT__,ii/col_bs);
1177: }
1178: MatGetOwnershipRangeColumn( Prol, &myCrs0, &kk );
1180: assert((kk-myCrs0)%col_bs==0);
1181: myCrs0 = myCrs0/col_bs;
1182: assert((kk/col_bs-myCrs0)==nLocalSelected);
1184: /* create global vector of data in 'data_w_ghost' */
1185: #if defined PETSC_GAMG_USE_LOG
1186: PetscLogEventBegin(petsc_gamg_setup_events[SET7],0,0,0,0);
1187: #endif
1188: if (npe > 1) { /* */
1189: PetscReal *tmp_gdata,*tmp_ldata,*tp2;
1190: PetscMalloc( nloc*sizeof(PetscReal), &tmp_ldata );
1191: for( jj = 0 ; jj < data_cols ; jj++ ){
1192: for( kk = 0 ; kk < bs ; kk++) {
1193: PetscInt ii,stride;
1194: const PetscReal *tp = pc_gamg->data + jj*bs*nloc + kk;
1195: for( ii = 0 ; ii < nloc ; ii++, tp += bs ){
1196: tmp_ldata[ii] = *tp;
1197: }
1198: PCGAMGGetDataWithGhosts( Gmat, 1, tmp_ldata, &stride, &tmp_gdata );
1199:
1201: if(jj==0 && kk==0) { /* now I know how many todal nodes - allocate */
1202: PetscMalloc( stride*bs*data_cols*sizeof(PetscReal), &data_w_ghost );
1203: nbnodes = bs*stride;
1204: }
1205: tp2 = data_w_ghost + jj*bs*stride + kk;
1206: for( ii = 0 ; ii < stride ; ii++, tp2 += bs ){
1207: *tp2 = tmp_gdata[ii];
1208: }
1209: PetscFree( tmp_gdata );
1210: }
1211: }
1212: PetscFree( tmp_ldata );
1213: }
1214: else {
1215: nbnodes = bs*nloc;
1216: data_w_ghost = (PetscReal*)pc_gamg->data;
1217: }
1218:
1219: /* get P0 */
1220: if( npe > 1 ){
1221: PetscReal *fid_glid_loc,*fiddata;
1222: PetscInt stride;
1223:
1224: PetscMalloc( nloc*sizeof(PetscReal), &fid_glid_loc );
1225: for(kk=0;kk<nloc;kk++) fid_glid_loc[kk] = (PetscReal)(my0+kk);
1226: PCGAMGGetDataWithGhosts( Gmat, 1, fid_glid_loc, &stride, &fiddata );
1227:
1228: PetscMalloc( stride*sizeof(PetscInt), &flid_fgid );
1229: for(kk=0;kk<stride;kk++) flid_fgid[kk] = (PetscInt)fiddata[kk];
1230: PetscFree( fiddata );
1232: assert(stride==nbnodes/bs);
1233: PetscFree( fid_glid_loc );
1234: }
1235: else {
1236: PetscMalloc( nloc*sizeof(PetscInt), &flid_fgid );
1237: for(kk=0;kk<nloc;kk++) flid_fgid[kk] = my0 + kk;
1238: }
1239: #if defined PETSC_GAMG_USE_LOG
1240: PetscLogEventEnd(petsc_gamg_setup_events[SET7],0,0,0,0);
1241: PetscLogEventBegin(petsc_gamg_setup_events[SET8],0,0,0,0);
1242: #endif
1243: {
1244: PetscReal *data_out = PETSC_NULL;
1245: formProl0( agg_lists, bs, data_cols, myCrs0, nbnodes,
1246: data_w_ghost, flid_fgid, &data_out, Prol );
1247:
1248: PetscFree( pc_gamg->data ); CHKERRQ( ierr );
1249: pc_gamg->data = data_out;
1250: pc_gamg->data_cell_rows = data_cols;
1251: pc_gamg->data_sz = data_cols*data_cols*nLocalSelected;
1252: }
1253: #if defined PETSC_GAMG_USE_LOG
1254: PetscLogEventEnd(petsc_gamg_setup_events[SET8],0,0,0,0);
1255: #endif
1256: if (npe > 1) PetscFree( data_w_ghost );
1257: PetscFree( flid_fgid );
1258:
1259: *a_P_out = Prol; /* out */
1260: #if defined PETSC_USE_LOG
1261: PetscLogEventEnd(PC_GAMGProlongator_AGG,0,0,0,0);
1262: #endif
1263: return(0);
1264: }
1266: /* -------------------------------------------------------------------------- */
1267: /*
1268: PCGAMGOptprol_AGG
1270: Input Parameter:
1271: . pc - this
1272: . Amat - matrix on this fine level
1273: In/Output Parameter:
1274: . a_P_out - prolongation operator to the next level
1275: */
1278: PetscErrorCode PCGAMGOptprol_AGG( PC pc,
1279: const Mat Amat,
1280: Mat *a_P
1281: )
1282: {
1284: PC_MG *mg = (PC_MG*)pc->data;
1285: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1286: const PetscInt verbose = pc_gamg->verbose;
1287: PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx;
1288: PetscInt jj;
1289: PetscMPIInt mype,npe;
1290: Mat Prol = *a_P;
1291: MPI_Comm wcomm = ((PetscObject)Amat)->comm;
1294: #if defined PETSC_USE_LOG
1295: PetscLogEventBegin(PC_GAMGOptprol_AGG,0,0,0,0);
1296: #endif
1297: MPI_Comm_rank( wcomm, &mype);
1298: MPI_Comm_size( wcomm, &npe);
1300: /* smooth P0 */
1301: for( jj = 0 ; jj < pc_gamg_agg->nsmooths ; jj++ ){
1302: Mat tMat;
1303: Vec diag;
1304: PetscReal alpha, emax, emin;
1305: #if defined PETSC_GAMG_USE_LOG
1306: PetscLogEventBegin(petsc_gamg_setup_events[SET9],0,0,0,0);
1307: #endif
1308: if( jj == 0 ) {
1309: KSP eksp;
1310: Vec bb, xx;
1311: PC pc;
1312: MatGetVecs( Amat, &bb, 0 );
1313: MatGetVecs( Amat, &xx, 0 );
1314: {
1315: PetscRandom rctx;
1316: PetscRandomCreate(wcomm,&rctx);
1317: PetscRandomSetFromOptions(rctx);
1318: VecSetRandom(bb,rctx);
1319: PetscRandomDestroy( &rctx );
1320: }
1321: KSPCreate(wcomm,&eksp);
1322: KSPAppendOptionsPrefix( eksp, "gamg_est_");
1323: KSPSetFromOptions( eksp );
1324: KSPSetInitialGuessNonzero( eksp, PETSC_FALSE );
1325: KSPSetOperators( eksp, Amat, Amat, SAME_NONZERO_PATTERN );
1326: CHKERRQ( ierr );
1327: KSPGetPC( eksp, &pc ); CHKERRQ( ierr );
1328: PCSetType( pc, PCJACOBI ); /* smoother */
1329: KSPSetTolerances(eksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
1330:
1331: KSPSetNormType( eksp, KSP_NORM_NONE );
1332: KSPSetComputeSingularValues( eksp,PETSC_TRUE );
1333:
1334: /* solve - keep stuff out of logging */
1335: PetscLogEventDeactivate(KSP_Solve);
1336: PetscLogEventDeactivate(PC_Apply);
1337: KSPSolve( eksp, bb, xx );
1338: PetscLogEventActivate(KSP_Solve);
1339: PetscLogEventActivate(PC_Apply);
1340:
1341: KSPComputeExtremeSingularValues( eksp, &emax, &emin );
1342: if( verbose ) {
1343: PetscPrintf(wcomm,"\t\t\t%s smooth P0: max eigen=%e min=%e PC=%s\n",
1344: __FUNCT__,emax,emin,PCJACOBI);
1345: }
1346: VecDestroy( &xx );
1347: VecDestroy( &bb );
1348: KSPDestroy( &eksp );
1350: if( pc_gamg->emax_id == -1 ) {
1351: PetscObjectComposedDataRegister( &pc_gamg->emax_id );
1352: assert(pc_gamg->emax_id != -1 );
1353: }
1354: PetscObjectComposedDataSetScalar( (PetscObject)Amat, pc_gamg->emax_id, emax );
1355: }
1357: /* smooth P1 := (I - omega/lam D^{-1}A)P0 */
1358: MatMatMult( Amat, Prol, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &tMat );
1359: MatGetVecs( Amat, &diag, 0 );
1360: MatGetDiagonal( Amat, diag ); /* effectively PCJACOBI */
1361: VecReciprocal( diag );
1362: MatDiagonalScale( tMat, diag, 0 );
1363: VecDestroy( &diag );
1364: alpha = -1.5/emax;
1365: MatAYPX( tMat, alpha, Prol, SUBSET_NONZERO_PATTERN );
1366: MatDestroy( &Prol );
1367: Prol = tMat;
1368: #if defined PETSC_GAMG_USE_LOG
1369: PetscLogEventEnd(petsc_gamg_setup_events[SET9],0,0,0,0);
1370: #endif
1371: }
1372: #if defined PETSC_USE_LOG
1373: PetscLogEventEnd(PC_GAMGOptprol_AGG,0,0,0,0);
1374: #endif
1375: *a_P = Prol;
1377: return(0);
1378: }
1380: /* -------------------------------------------------------------------------- */
1381: /*
1382: PCGAMGKKTProl_AGG
1384: Input Parameter:
1385: . pc - this
1386: . Prol11 - matrix on this fine level
1387: . A21 - matrix on this fine level
1388: In/Output Parameter:
1389: . a_P22 - prolongation operator to the next level
1390: */
1393: PetscErrorCode PCGAMGKKTProl_AGG( PC pc,
1394: const Mat Prol11,
1395: const Mat A21,
1396: Mat *a_P22
1397: )
1398: {
1400: PC_MG *mg = (PC_MG*)pc->data;
1401: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1402: const PetscInt verbose = pc_gamg->verbose;
1403: /* PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG*)pc_gamg->subctx; */
1404: PetscMPIInt mype,npe;
1405: Mat Prol22,Tmat,Gmat;
1406: MPI_Comm wcomm = ((PetscObject)pc)->comm;
1407: PetscCoarsenData *agg_lists;
1410: #if defined PETSC_USE_LOG
1411: PetscLogEventBegin(PC_GAMGKKTProl_AGG,0,0,0,0);
1412: #endif
1413: MPI_Comm_rank( wcomm, &mype);
1414: MPI_Comm_size( wcomm, &npe);
1416: /* form C graph */
1417: MatMatMult( A21, Prol11, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Tmat);
1418: MatMatTransposeMult(Tmat, Tmat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Gmat );
1419: MatDestroy(&Tmat);
1420: PCGAMGFilterGraph(&Gmat, 0.0, PETSC_FALSE, verbose);
1422: /* coarsen constraints */
1423: {
1424: MatCoarsen crs;
1425: MatCoarsenCreate( wcomm, &crs );
1426: MatCoarsenSetType( crs, MATCOARSENMIS );
1427: MatCoarsenSetAdjacency( crs, Gmat );
1428: MatCoarsenSetVerbose( crs, verbose );
1429: MatCoarsenSetStrictAggs( crs, PETSC_TRUE );
1430: MatCoarsenApply( crs );
1431: MatCoarsenGetData( crs, &agg_lists );
1432: MatCoarsenDestroy( &crs );
1433: }
1435: /* form simple prolongation 'Prol22' */
1436: {
1437: PetscInt ii,mm,clid,my0,nloc,nLocalSelected;
1438: PetscScalar val = 1.0;
1439: /* get 'nLocalSelected' */
1440: MatGetLocalSize( Gmat, &nloc, &ii );
1441: for( ii=0, nLocalSelected = 0 ; ii < nloc ; ii++ ){
1442: PetscBool ise;
1443: /* filter out singletons 0 or 1? */
1444: PetscCDEmptyAt( agg_lists, ii, &ise );
1445: if( !ise ) nLocalSelected++;
1446: }
1448: MatCreate(wcomm,&Prol22);
1449: MatSetSizes( Prol22,nloc, nLocalSelected,
1450: PETSC_DETERMINE, PETSC_DETERMINE);
1451:
1452: MatSetType( Prol22, MATAIJ );
1453: MatSeqAIJSetPreallocation(Prol22,1,PETSC_NULL);
1454: MatMPIAIJSetPreallocation(Prol22,1,PETSC_NULL,1,PETSC_NULL);
1455:
1456: /* MatCreateAIJ( wcomm, */
1457: /* nloc, nLocalSelected, */
1458: /* PETSC_DETERMINE, PETSC_DETERMINE, */
1459: /* 1, PETSC_NULL, 1, PETSC_NULL, */
1460: /* &Prol22 ); */
1461:
1462: MatGetOwnershipRange( Prol22, &my0, &ii );
1463: nloc = ii - my0;
1464:
1465: /* make aggregates */
1466: for( mm = clid = 0 ; mm < nloc ; mm++ ){
1467: PetscCDSizeAt( agg_lists, mm, &ii );
1468: if( ii > 0 ) {
1469: PetscInt asz=ii,cgid=my0+clid,rids[1000];
1470: PetscCDPos pos;
1471: if(asz>1000)SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Very large aggregate: %d",asz);
1472: ii = 0;
1473: PetscCDGetHeadPos(agg_lists,mm,&pos);
1474: while(pos){
1475: PetscInt gid1;
1476: PetscLLNGetID( pos, &gid1 );
1477: PetscCDGetNextPos(agg_lists,mm,&pos);
1478:
1479: rids[ii++] = gid1;
1480: }
1481: assert(ii==asz);
1482: /* add diagonal block of P0 */
1483: MatSetValues(Prol22,asz,rids,1,&cgid,&val,INSERT_VALUES);
1484:
1485: clid++;
1486: } /* coarse agg */
1487: } /* for all fine nodes */
1488: MatAssemblyBegin(Prol22,MAT_FINAL_ASSEMBLY);
1489: MatAssemblyEnd(Prol22,MAT_FINAL_ASSEMBLY);
1490: }
1492: /* clean up */
1493: MatDestroy( &Gmat );
1494: PetscCDDestroy( agg_lists );
1495: #if defined PETSC_USE_LOG
1496: PetscLogEventEnd(PC_GAMGKKTProl_AGG,0,0,0,0);
1497: #endif
1498: *a_P22 = Prol22;
1500: return(0);
1501: }
1503: /* -------------------------------------------------------------------------- */
1504: /*
1505: PCCreateGAMG_AGG
1507: Input Parameter:
1508: . pc -
1509: */
1512: PetscErrorCode PCCreateGAMG_AGG( PC pc )
1513: {
1514: PetscErrorCode ierr;
1515: PC_MG *mg = (PC_MG*)pc->data;
1516: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1517: PC_GAMG_AGG *pc_gamg_agg;
1520: /* create sub context for SA */
1521: PetscNewLog( pc, PC_GAMG_AGG, &pc_gamg_agg );
1522: assert(!pc_gamg->subctx);
1523: pc_gamg->subctx = pc_gamg_agg;
1524:
1525: pc->ops->setfromoptions = PCSetFromOptions_GAMG_AGG;
1526: pc->ops->destroy = PCDestroy_AGG;
1527: /* reset does not do anything; setup not virtual */
1529: /* set internal function pointers */
1530: pc_gamg->graph = PCGAMGgraph_AGG;
1531: pc_gamg->coarsen = PCGAMGCoarsen_AGG;
1532: pc_gamg->prolongator = PCGAMGProlongator_AGG;
1533: pc_gamg->optprol = PCGAMGOptprol_AGG;
1534: pc_gamg->formkktprol = PCGAMGKKTProl_AGG;
1535:
1536: pc_gamg->createdefaultdata = PCSetData_AGG;
1538: PetscObjectComposeFunctionDynamic( (PetscObject)pc,
1539: "PCSetCoordinates_C",
1540: "PCSetCoordinates_AGG",
1541: PCSetCoordinates_AGG);
1542: return(0);
1543: }