Actual source code: classical.c
petsc-3.8.4 2018-03-24
1: #include <../src/ksp/pc/impls/gamg/gamg.h>
2: #include <petscsf.h>
4: PetscFunctionList PCGAMGClassicalProlongatorList = NULL;
5: PetscBool PCGAMGClassicalPackageInitialized = PETSC_FALSE;
7: typedef struct {
8: PetscReal interp_threshold; /* interpolation threshold */
9: char prolongtype[256];
10: PetscInt nsmooths; /* number of jacobi smoothings on the prolongator */
11: } PC_GAMG_Classical;
13: /*@C
14: PCGAMGClassicalSetType - Sets the type of classical interpolation to use
16: Collective on PC
18: Input Parameters:
19: . pc - the preconditioner context
21: Options Database Key:
22: . -pc_gamg_classical_type
24: Level: intermediate
26: .seealso: ()
27: @*/
28: PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type)
29: {
34: PetscTryMethod(pc,"PCGAMGClassicalSetType_C",(PC,PCGAMGType),(pc,type));
35: return(0);
36: }
38: /*@C
39: PCGAMGClassicalGetType - Gets the type of classical interpolation to use
41: Collective on PC
43: Input Parameter:
44: . pc - the preconditioner context
46: Output Parameter:
47: . type - the type used
49: Level: intermediate
51: .seealso: ()
52: @*/
53: PetscErrorCode PCGAMGClassicalGetType(PC pc, PCGAMGClassicalType *type)
54: {
59: PetscUseMethod(pc,"PCGAMGClassicalGetType_C",(PC,PCGAMGType*),(pc,type));
60: return(0);
61: }
63: static PetscErrorCode PCGAMGClassicalSetType_GAMG(PC pc, PCGAMGClassicalType type)
64: {
65: PetscErrorCode ierr;
66: PC_MG *mg = (PC_MG*)pc->data;
67: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
68: PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
71: PetscStrcpy(cls->prolongtype,type);
72: return(0);
73: }
75: static PetscErrorCode PCGAMGClassicalGetType_GAMG(PC pc, PCGAMGClassicalType *type)
76: {
77: PC_MG *mg = (PC_MG*)pc->data;
78: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
79: PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
82: *type = cls->prolongtype;
83: return(0);
84: }
86: PetscErrorCode PCGAMGGraph_Classical(PC pc,Mat A,Mat *G)
87: {
88: PetscInt s,f,n,idx,lidx,gidx;
89: PetscInt r,c,ncols;
90: const PetscInt *rcol;
91: const PetscScalar *rval;
92: PetscInt *gcol;
93: PetscScalar *gval;
94: PetscReal rmax;
95: PetscInt cmax = 0;
96: PC_MG *mg;
97: PC_GAMG *gamg;
98: PetscErrorCode ierr;
99: PetscInt *gsparse,*lsparse;
100: PetscScalar *Amax;
101: MatType mtype;
104: mg = (PC_MG *)pc->data;
105: gamg = (PC_GAMG *)mg->innerctx;
107: MatGetOwnershipRange(A,&s,&f);
108: n=f-s;
109: PetscMalloc3(n,&lsparse,n,&gsparse,n,&Amax);
111: for (r = 0;r < n;r++) {
112: lsparse[r] = 0;
113: gsparse[r] = 0;
114: }
116: for (r = s;r < f;r++) {
117: /* determine the maximum off-diagonal in each row */
118: rmax = 0.;
119: MatGetRow(A,r,&ncols,&rcol,&rval);
120: for (c = 0; c < ncols; c++) {
121: if (PetscRealPart(-rval[c]) > rmax && rcol[c] != r) {
122: rmax = PetscRealPart(-rval[c]);
123: }
124: }
125: Amax[r-s] = rmax;
126: if (ncols > cmax) cmax = ncols;
127: lidx = 0;
128: gidx = 0;
129: /* create the local and global sparsity patterns */
130: for (c = 0; c < ncols; c++) {
131: if (PetscRealPart(-rval[c]) > gamg->threshold[0]*PetscRealPart(Amax[r-s]) || rcol[c] == r) {
132: if (rcol[c] < f && rcol[c] >= s) {
133: lidx++;
134: } else {
135: gidx++;
136: }
137: }
138: }
139: MatRestoreRow(A,r,&ncols,&rcol,&rval);
140: lsparse[r-s] = lidx;
141: gsparse[r-s] = gidx;
142: }
143: PetscMalloc2(cmax,&gval,cmax,&gcol);
145: MatCreate(PetscObjectComm((PetscObject)A),G);
146: MatGetType(A,&mtype);
147: MatSetType(*G,mtype);
148: MatSetSizes(*G,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
149: MatMPIAIJSetPreallocation(*G,0,lsparse,0,gsparse);
150: MatSeqAIJSetPreallocation(*G,0,lsparse);
151: for (r = s;r < f;r++) {
152: MatGetRow(A,r,&ncols,&rcol,&rval);
153: idx = 0;
154: for (c = 0; c < ncols; c++) {
155: /* classical strength of connection */
156: if (PetscRealPart(-rval[c]) > gamg->threshold[0]*PetscRealPart(Amax[r-s]) || rcol[c] == r) {
157: gcol[idx] = rcol[c];
158: gval[idx] = rval[c];
159: idx++;
160: }
161: }
162: MatSetValues(*G,1,&r,idx,gcol,gval,INSERT_VALUES);
163: MatRestoreRow(A,r,&ncols,&rcol,&rval);
164: }
165: MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY);
166: MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY);
168: PetscFree2(gval,gcol);
169: PetscFree3(lsparse,gsparse,Amax);
170: return(0);
171: }
174: PetscErrorCode PCGAMGCoarsen_Classical(PC pc,Mat *G,PetscCoarsenData **agg_lists)
175: {
176: PetscErrorCode ierr;
177: MatCoarsen crs;
178: MPI_Comm fcomm = ((PetscObject)pc)->comm;
181: if (!G) SETERRQ(fcomm,PETSC_ERR_ARG_WRONGSTATE,"Must set Graph in PC in PCGAMG before coarsening");
183: MatCoarsenCreate(fcomm,&crs);
184: MatCoarsenSetFromOptions(crs);
185: MatCoarsenSetAdjacency(crs,*G);
186: MatCoarsenSetStrictAggs(crs,PETSC_TRUE);
187: MatCoarsenApply(crs);
188: MatCoarsenGetData(crs,agg_lists);
189: MatCoarsenDestroy(&crs);
191: return(0);
192: }
194: PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
195: {
196: PetscErrorCode ierr;
197: PC_MG *mg = (PC_MG*)pc->data;
198: PC_GAMG *gamg = (PC_GAMG*)mg->innerctx;
199: PetscBool iscoarse,isMPIAIJ,isSEQAIJ;
200: PetscInt fn,cn,fs,fe,cs,ce,i,j,ncols,col,row_f,row_c,cmax=0,idx,noff;
201: PetscInt *lcid,*gcid,*lsparse,*gsparse,*colmap,*pcols;
202: const PetscInt *rcol;
203: PetscReal *Amax_pos,*Amax_neg;
204: PetscScalar g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta,pij;
205: PetscScalar *pvals;
206: const PetscScalar *rval;
207: Mat lA,gA=NULL;
208: MatType mtype;
209: Vec C,lvec;
210: PetscLayout clayout;
211: PetscSF sf;
212: Mat_MPIAIJ *mpiaij;
215: MatGetOwnershipRange(A,&fs,&fe);
216: fn = fe-fs;
217: PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);
218: PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSEQAIJ);
219: if (!isMPIAIJ && !isSEQAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Classical AMG requires MPIAIJ matrix");
220: if (isMPIAIJ) {
221: mpiaij = (Mat_MPIAIJ*)A->data;
222: lA = mpiaij->A;
223: gA = mpiaij->B;
224: lvec = mpiaij->lvec;
225: VecGetSize(lvec,&noff);
226: colmap = mpiaij->garray;
227: MatGetLayouts(A,NULL,&clayout);
228: PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
229: PetscSFSetGraphLayout(sf,clayout,noff,NULL,PETSC_COPY_VALUES,colmap);
230: PetscMalloc1(noff,&gcid);
231: } else {
232: lA = A;
233: }
234: PetscMalloc5(fn,&lsparse,fn,&gsparse,fn,&lcid,fn,&Amax_pos,fn,&Amax_neg);
236: /* count the number of coarse unknowns */
237: cn = 0;
238: for (i=0;i<fn;i++) {
239: /* filter out singletons */
240: PetscCDEmptyAt(agg_lists,i,&iscoarse);
241: lcid[i] = -1;
242: if (!iscoarse) {
243: cn++;
244: }
245: }
247: /* create the coarse vector */
248: VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&C);
249: VecGetOwnershipRange(C,&cs,&ce);
251: cn = 0;
252: for (i=0;i<fn;i++) {
253: PetscCDEmptyAt(agg_lists,i,&iscoarse);
254: if (!iscoarse) {
255: lcid[i] = cs+cn;
256: cn++;
257: } else {
258: lcid[i] = -1;
259: }
260: }
262: if (gA) {
263: PetscSFBcastBegin(sf,MPIU_INT,lcid,gcid);
264: PetscSFBcastEnd(sf,MPIU_INT,lcid,gcid);
265: }
267: /* determine the biggest off-diagonal entries in each row */
268: for (i=fs;i<fe;i++) {
269: Amax_pos[i-fs] = 0.;
270: Amax_neg[i-fs] = 0.;
271: MatGetRow(A,i,&ncols,&rcol,&rval);
272: for(j=0;j<ncols;j++){
273: if ((PetscRealPart(-rval[j]) > Amax_neg[i-fs]) && i != rcol[j]) Amax_neg[i-fs] = PetscAbsScalar(rval[j]);
274: if ((PetscRealPart(rval[j]) > Amax_pos[i-fs]) && i != rcol[j]) Amax_pos[i-fs] = PetscAbsScalar(rval[j]);
275: }
276: if (ncols > cmax) cmax = ncols;
277: MatRestoreRow(A,i,&ncols,&rcol,&rval);
278: }
279: PetscMalloc2(cmax,&pcols,cmax,&pvals);
280: VecDestroy(&C);
282: /* count the on and off processor sparsity patterns for the prolongator */
283: for (i=0;i<fn;i++) {
284: /* on */
285: lsparse[i] = 0;
286: gsparse[i] = 0;
287: if (lcid[i] >= 0) {
288: lsparse[i] = 1;
289: gsparse[i] = 0;
290: } else {
291: MatGetRow(lA,i,&ncols,&rcol,&rval);
292: for (j = 0;j < ncols;j++) {
293: col = rcol[j];
294: if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
295: lsparse[i] += 1;
296: }
297: }
298: MatRestoreRow(lA,i,&ncols,&rcol,&rval);
299: /* off */
300: if (gA) {
301: MatGetRow(gA,i,&ncols,&rcol,&rval);
302: for (j = 0; j < ncols; j++) {
303: col = rcol[j];
304: if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
305: gsparse[i] += 1;
306: }
307: }
308: MatRestoreRow(gA,i,&ncols,&rcol,&rval);
309: }
310: }
311: }
313: /* preallocate and create the prolongator */
314: MatCreate(PetscObjectComm((PetscObject)A),P);
315: MatGetType(G,&mtype);
316: MatSetType(*P,mtype);
317: MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);
318: MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);
319: MatSeqAIJSetPreallocation(*P,0,lsparse);
321: /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
322: for (i = 0;i < fn;i++) {
323: /* determine on or off */
324: row_f = i + fs;
325: row_c = lcid[i];
326: if (row_c >= 0) {
327: pij = 1.;
328: MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);
329: } else {
330: g_pos = 0.;
331: g_neg = 0.;
332: a_pos = 0.;
333: a_neg = 0.;
334: diag = 0.;
336: /* local connections */
337: MatGetRow(lA,i,&ncols,&rcol,&rval);
338: for (j = 0; j < ncols; j++) {
339: col = rcol[j];
340: if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
341: if (PetscRealPart(rval[j]) > 0.) {
342: g_pos += rval[j];
343: } else {
344: g_neg += rval[j];
345: }
346: }
347: if (col != i) {
348: if (PetscRealPart(rval[j]) > 0.) {
349: a_pos += rval[j];
350: } else {
351: a_neg += rval[j];
352: }
353: } else {
354: diag = rval[j];
355: }
356: }
357: MatRestoreRow(lA,i,&ncols,&rcol,&rval);
359: /* ghosted connections */
360: if (gA) {
361: MatGetRow(gA,i,&ncols,&rcol,&rval);
362: for (j = 0; j < ncols; j++) {
363: col = rcol[j];
364: if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
365: if (PetscRealPart(rval[j]) > 0.) {
366: g_pos += rval[j];
367: } else {
368: g_neg += rval[j];
369: }
370: }
371: if (PetscRealPart(rval[j]) > 0.) {
372: a_pos += rval[j];
373: } else {
374: a_neg += rval[j];
375: }
376: }
377: MatRestoreRow(gA,i,&ncols,&rcol,&rval);
378: }
380: if (g_neg == 0.) {
381: alpha = 0.;
382: } else {
383: alpha = -a_neg/g_neg;
384: }
386: if (g_pos == 0.) {
387: diag += a_pos;
388: beta = 0.;
389: } else {
390: beta = -a_pos/g_pos;
391: }
392: if (diag == 0.) {
393: invdiag = 0.;
394: } else invdiag = 1. / diag;
395: /* on */
396: MatGetRow(lA,i,&ncols,&rcol,&rval);
397: idx = 0;
398: for (j = 0;j < ncols;j++) {
399: col = rcol[j];
400: if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
401: row_f = i + fs;
402: row_c = lcid[col];
403: /* set the values for on-processor ones */
404: if (PetscRealPart(rval[j]) < 0.) {
405: pij = rval[j]*alpha*invdiag;
406: } else {
407: pij = rval[j]*beta*invdiag;
408: }
409: if (PetscAbsScalar(pij) != 0.) {
410: pvals[idx] = pij;
411: pcols[idx] = row_c;
412: idx++;
413: }
414: }
415: }
416: MatRestoreRow(lA,i,&ncols,&rcol,&rval);
417: /* off */
418: if (gA) {
419: MatGetRow(gA,i,&ncols,&rcol,&rval);
420: for (j = 0; j < ncols; j++) {
421: col = rcol[j];
422: if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold[0]*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold[0]*Amax_neg[i])) {
423: row_f = i + fs;
424: row_c = gcid[col];
425: /* set the values for on-processor ones */
426: if (PetscRealPart(rval[j]) < 0.) {
427: pij = rval[j]*alpha*invdiag;
428: } else {
429: pij = rval[j]*beta*invdiag;
430: }
431: if (PetscAbsScalar(pij) != 0.) {
432: pvals[idx] = pij;
433: pcols[idx] = row_c;
434: idx++;
435: }
436: }
437: }
438: MatRestoreRow(gA,i,&ncols,&rcol,&rval);
439: }
440: MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);
441: }
442: }
444: MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
445: MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);
447: PetscFree5(lsparse,gsparse,lcid,Amax_pos,Amax_neg);
449: PetscFree2(pcols,pvals);
450: if (gA) {
451: PetscSFDestroy(&sf);
452: PetscFree(gcid);
453: }
455: return(0);
456: }
458: PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc,Mat *P)
459: {
460: PetscInt j,i,ps,pf,pn,pcs,pcf,pcn,idx,cmax;
461: PetscErrorCode ierr;
462: const PetscScalar *pval;
463: const PetscInt *pcol;
464: PetscScalar *pnval;
465: PetscInt *pncol;
466: PetscInt ncols;
467: Mat Pnew;
468: PetscInt *lsparse,*gsparse;
469: PetscReal pmax_pos,pmax_neg,ptot_pos,ptot_neg,pthresh_pos,pthresh_neg;
470: PC_MG *mg = (PC_MG*)pc->data;
471: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
472: PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
473: MatType mtype;
476: /* trim and rescale with reallocation */
477: MatGetOwnershipRange(*P,&ps,&pf);
478: MatGetOwnershipRangeColumn(*P,&pcs,&pcf);
479: pn = pf-ps;
480: pcn = pcf-pcs;
481: PetscMalloc2(pn,&lsparse,pn,&gsparse);
482: /* allocate */
483: cmax = 0;
484: for (i=ps;i<pf;i++) {
485: lsparse[i-ps] = 0;
486: gsparse[i-ps] = 0;
487: MatGetRow(*P,i,&ncols,&pcol,&pval);
488: if (ncols > cmax) {
489: cmax = ncols;
490: }
491: pmax_pos = 0.;
492: pmax_neg = 0.;
493: for (j=0;j<ncols;j++) {
494: if (PetscRealPart(pval[j]) > pmax_pos) {
495: pmax_pos = PetscRealPart(pval[j]);
496: } else if (PetscRealPart(pval[j]) < pmax_neg) {
497: pmax_neg = PetscRealPart(pval[j]);
498: }
499: }
500: for (j=0;j<ncols;j++) {
501: if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
502: if (pcol[j] >= pcs && pcol[j] < pcf) {
503: lsparse[i-ps]++;
504: } else {
505: gsparse[i-ps]++;
506: }
507: }
508: }
509: MatRestoreRow(*P,i,&ncols,&pcol,&pval);
510: }
512: PetscMalloc2(cmax,&pnval,cmax,&pncol);
514: MatGetType(*P,&mtype);
515: MatCreate(PetscObjectComm((PetscObject)*P),&Pnew);
516: MatSetType(Pnew, mtype);
517: MatSetSizes(Pnew,pn,pcn,PETSC_DETERMINE,PETSC_DETERMINE);
518: MatSeqAIJSetPreallocation(Pnew,0,lsparse);
519: MatMPIAIJSetPreallocation(Pnew,0,lsparse,0,gsparse);
521: for (i=ps;i<pf;i++) {
522: MatGetRow(*P,i,&ncols,&pcol,&pval);
523: pmax_pos = 0.;
524: pmax_neg = 0.;
525: for (j=0;j<ncols;j++) {
526: if (PetscRealPart(pval[j]) > pmax_pos) {
527: pmax_pos = PetscRealPart(pval[j]);
528: } else if (PetscRealPart(pval[j]) < pmax_neg) {
529: pmax_neg = PetscRealPart(pval[j]);
530: }
531: }
532: pthresh_pos = 0.;
533: pthresh_neg = 0.;
534: ptot_pos = 0.;
535: ptot_neg = 0.;
536: for (j=0;j<ncols;j++) {
537: if (PetscRealPart(pval[j]) >= cls->interp_threshold*pmax_pos) {
538: pthresh_pos += PetscRealPart(pval[j]);
539: } else if (PetscRealPart(pval[j]) <= cls->interp_threshold*pmax_neg) {
540: pthresh_neg += PetscRealPart(pval[j]);
541: }
542: if (PetscRealPart(pval[j]) > 0.) {
543: ptot_pos += PetscRealPart(pval[j]);
544: } else {
545: ptot_neg += PetscRealPart(pval[j]);
546: }
547: }
548: if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
549: if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
550: idx=0;
551: for (j=0;j<ncols;j++) {
552: if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold) {
553: pnval[idx] = ptot_pos*pval[j];
554: pncol[idx] = pcol[j];
555: idx++;
556: } else if (PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
557: pnval[idx] = ptot_neg*pval[j];
558: pncol[idx] = pcol[j];
559: idx++;
560: }
561: }
562: MatRestoreRow(*P,i,&ncols,&pcol,&pval);
563: MatSetValues(Pnew,1,&i,idx,pncol,pnval,INSERT_VALUES);
564: }
566: MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY);
567: MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY);
568: MatDestroy(P);
570: *P = Pnew;
571: PetscFree2(lsparse,gsparse);
572: PetscFree2(pnval,pncol);
573: return(0);
574: }
576: PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
577: {
578: PetscErrorCode ierr;
579: Mat lA,*lAs;
580: MatType mtype;
581: Vec cv;
582: PetscInt *gcid,*lcid,*lsparse,*gsparse,*picol;
583: PetscInt fs,fe,cs,ce,nl,i,j,k,li,lni,ci,ncols,maxcols,fn,cn,cid;
584: PetscMPIInt size;
585: const PetscInt *lidx,*icol,*gidx;
586: PetscBool iscoarse;
587: PetscScalar vi,pentry,pjentry;
588: PetscScalar *pcontrib,*pvcol;
589: const PetscScalar *vcol;
590: PetscReal diag,jdiag,jwttotal;
591: PetscInt pncols;
592: PetscSF sf;
593: PetscLayout clayout;
594: IS lis;
597: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
598: MatGetOwnershipRange(A,&fs,&fe);
599: fn = fe-fs;
600: ISCreateStride(PETSC_COMM_SELF,fe-fs,fs,1,&lis);
601: if (size > 1) {
602: MatGetLayouts(A,NULL,&clayout);
603: /* increase the overlap by two to get neighbors of neighbors */
604: MatIncreaseOverlap(A,1,&lis,2);
605: ISSort(lis);
606: /* get the local part of A */
607: MatCreateSubMatrices(A,1,&lis,&lis,MAT_INITIAL_MATRIX,&lAs);
608: lA = lAs[0];
609: /* build an SF out of it */
610: ISGetLocalSize(lis,&nl);
611: PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
612: ISGetIndices(lis,&lidx);
613: PetscSFSetGraphLayout(sf,clayout,nl,NULL,PETSC_COPY_VALUES,lidx);
614: ISRestoreIndices(lis,&lidx);
615: } else {
616: lA = A;
617: nl = fn;
618: }
619: /* create a communication structure for the overlapped portion and transmit coarse indices */
620: PetscMalloc3(fn,&lsparse,fn,&gsparse,nl,&pcontrib);
621: /* create coarse vector */
622: cn = 0;
623: for (i=0;i<fn;i++) {
624: PetscCDEmptyAt(agg_lists,i,&iscoarse);
625: if (!iscoarse) {
626: cn++;
627: }
628: }
629: PetscMalloc1(fn,&gcid);
630: VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&cv);
631: VecGetOwnershipRange(cv,&cs,&ce);
632: cn = 0;
633: for (i=0;i<fn;i++) {
634: PetscCDEmptyAt(agg_lists,i,&iscoarse);
635: if (!iscoarse) {
636: gcid[i] = cs+cn;
637: cn++;
638: } else {
639: gcid[i] = -1;
640: }
641: }
642: if (size > 1) {
643: PetscMalloc1(nl,&lcid);
644: PetscSFBcastBegin(sf,MPIU_INT,gcid,lcid);
645: PetscSFBcastEnd(sf,MPIU_INT,gcid,lcid);
646: } else {
647: lcid = gcid;
648: }
649: /* count to preallocate the prolongator */
650: ISGetIndices(lis,&gidx);
651: maxcols = 0;
652: /* count the number of unique contributing coarse cells for each fine */
653: for (i=0;i<nl;i++) {
654: pcontrib[i] = 0.;
655: MatGetRow(lA,i,&ncols,&icol,NULL);
656: if (gidx[i] >= fs && gidx[i] < fe) {
657: li = gidx[i] - fs;
658: lsparse[li] = 0;
659: gsparse[li] = 0;
660: cid = lcid[i];
661: if (cid >= 0) {
662: lsparse[li] = 1;
663: } else {
664: for (j=0;j<ncols;j++) {
665: if (lcid[icol[j]] >= 0) {
666: pcontrib[icol[j]] = 1.;
667: } else {
668: ci = icol[j];
669: MatRestoreRow(lA,i,&ncols,&icol,NULL);
670: MatGetRow(lA,ci,&ncols,&icol,NULL);
671: for (k=0;k<ncols;k++) {
672: if (lcid[icol[k]] >= 0) {
673: pcontrib[icol[k]] = 1.;
674: }
675: }
676: MatRestoreRow(lA,ci,&ncols,&icol,NULL);
677: MatGetRow(lA,i,&ncols,&icol,NULL);
678: }
679: }
680: for (j=0;j<ncols;j++) {
681: if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
682: lni = lcid[icol[j]];
683: if (lni >= cs && lni < ce) {
684: lsparse[li]++;
685: } else {
686: gsparse[li]++;
687: }
688: pcontrib[icol[j]] = 0.;
689: } else {
690: ci = icol[j];
691: MatRestoreRow(lA,i,&ncols,&icol,NULL);
692: MatGetRow(lA,ci,&ncols,&icol,NULL);
693: for (k=0;k<ncols;k++) {
694: if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
695: lni = lcid[icol[k]];
696: if (lni >= cs && lni < ce) {
697: lsparse[li]++;
698: } else {
699: gsparse[li]++;
700: }
701: pcontrib[icol[k]] = 0.;
702: }
703: }
704: MatRestoreRow(lA,ci,&ncols,&icol,NULL);
705: MatGetRow(lA,i,&ncols,&icol,NULL);
706: }
707: }
708: }
709: if (lsparse[li] + gsparse[li] > maxcols) maxcols = lsparse[li]+gsparse[li];
710: }
711: MatRestoreRow(lA,i,&ncols,&icol,&vcol);
712: }
713: PetscMalloc2(maxcols,&picol,maxcols,&pvcol);
714: MatCreate(PetscObjectComm((PetscObject)A),P);
715: MatGetType(A,&mtype);
716: MatSetType(*P,mtype);
717: MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);
718: MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);
719: MatSeqAIJSetPreallocation(*P,0,lsparse);
720: for (i=0;i<nl;i++) {
721: diag = 0.;
722: if (gidx[i] >= fs && gidx[i] < fe) {
723: pncols=0;
724: cid = lcid[i];
725: if (cid >= 0) {
726: pncols = 1;
727: picol[0] = cid;
728: pvcol[0] = 1.;
729: } else {
730: MatGetRow(lA,i,&ncols,&icol,&vcol);
731: for (j=0;j<ncols;j++) {
732: pentry = vcol[j];
733: if (lcid[icol[j]] >= 0) {
734: /* coarse neighbor */
735: pcontrib[icol[j]] += pentry;
736: } else if (icol[j] != i) {
737: /* the neighbor is a strongly connected fine node */
738: ci = icol[j];
739: vi = vcol[j];
740: MatRestoreRow(lA,i,&ncols,&icol,&vcol);
741: MatGetRow(lA,ci,&ncols,&icol,&vcol);
742: jwttotal=0.;
743: jdiag = 0.;
744: for (k=0;k<ncols;k++) {
745: if (ci == icol[k]) {
746: jdiag = PetscRealPart(vcol[k]);
747: }
748: }
749: for (k=0;k<ncols;k++) {
750: if (lcid[icol[k]] >= 0 && jdiag*PetscRealPart(vcol[k]) < 0.) {
751: pjentry = vcol[k];
752: jwttotal += PetscRealPart(pjentry);
753: }
754: }
755: if (jwttotal != 0.) {
756: jwttotal = PetscRealPart(vi)/jwttotal;
757: for (k=0;k<ncols;k++) {
758: if (lcid[icol[k]] >= 0 && jdiag*PetscRealPart(vcol[k]) < 0.) {
759: pjentry = vcol[k]*jwttotal;
760: pcontrib[icol[k]] += pjentry;
761: }
762: }
763: } else {
764: diag += PetscRealPart(vi);
765: }
766: MatRestoreRow(lA,ci,&ncols,&icol,&vcol);
767: MatGetRow(lA,i,&ncols,&icol,&vcol);
768: } else {
769: diag += PetscRealPart(vcol[j]);
770: }
771: }
772: if (diag != 0.) {
773: diag = 1./diag;
774: for (j=0;j<ncols;j++) {
775: if (lcid[icol[j]] >= 0 && pcontrib[icol[j]] != 0.) {
776: /* the neighbor is a coarse node */
777: if (PetscAbsScalar(pcontrib[icol[j]]) > 0.0) {
778: lni = lcid[icol[j]];
779: pvcol[pncols] = -pcontrib[icol[j]]*diag;
780: picol[pncols] = lni;
781: pncols++;
782: }
783: pcontrib[icol[j]] = 0.;
784: } else {
785: /* the neighbor is a strongly connected fine node */
786: ci = icol[j];
787: MatRestoreRow(lA,i,&ncols,&icol,&vcol);
788: MatGetRow(lA,ci,&ncols,&icol,&vcol);
789: for (k=0;k<ncols;k++) {
790: if (lcid[icol[k]] >= 0 && pcontrib[icol[k]] != 0.) {
791: if (PetscAbsScalar(pcontrib[icol[k]]) > 0.0) {
792: lni = lcid[icol[k]];
793: pvcol[pncols] = -pcontrib[icol[k]]*diag;
794: picol[pncols] = lni;
795: pncols++;
796: }
797: pcontrib[icol[k]] = 0.;
798: }
799: }
800: MatRestoreRow(lA,ci,&ncols,&icol,&vcol);
801: MatGetRow(lA,i,&ncols,&icol,&vcol);
802: }
803: pcontrib[icol[j]] = 0.;
804: }
805: MatRestoreRow(lA,i,&ncols,&icol,&vcol);
806: }
807: }
808: ci = gidx[i];
809: if (pncols > 0) {
810: MatSetValues(*P,1,&ci,pncols,picol,pvcol,INSERT_VALUES);
811: }
812: }
813: }
814: ISRestoreIndices(lis,&gidx);
815: PetscFree2(picol,pvcol);
816: PetscFree3(lsparse,gsparse,pcontrib);
817: ISDestroy(&lis);
818: PetscFree(gcid);
819: if (size > 1) {
820: PetscFree(lcid);
821: MatDestroyMatrices(1,&lAs);
822: PetscSFDestroy(&sf);
823: }
824: VecDestroy(&cv);
825: MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);
826: MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);
827: return(0);
828: }
830: PetscErrorCode PCGAMGOptProlongator_Classical_Jacobi(PC pc,Mat A,Mat *P)
831: {
833: PetscErrorCode ierr;
834: PetscInt f,s,n,cf,cs,i,idx;
835: PetscInt *coarserows;
836: PetscInt ncols;
837: const PetscInt *pcols;
838: const PetscScalar *pvals;
839: Mat Pnew;
840: Vec diag;
841: PC_MG *mg = (PC_MG*)pc->data;
842: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
843: PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
846: if (cls->nsmooths == 0) {
847: PCGAMGTruncateProlongator_Private(pc,P);
848: return(0);
849: }
850: MatGetOwnershipRange(*P,&s,&f);
851: n = f-s;
852: MatGetOwnershipRangeColumn(*P,&cs,&cf);
853: PetscMalloc1(n,&coarserows);
854: /* identify the rows corresponding to coarse unknowns */
855: idx = 0;
856: for (i=s;i<f;i++) {
857: MatGetRow(*P,i,&ncols,&pcols,&pvals);
858: /* assume, for now, that it's a coarse unknown if it has a single unit entry */
859: if (ncols == 1) {
860: if (pvals[0] == 1.) {
861: coarserows[idx] = i;
862: idx++;
863: }
864: }
865: MatRestoreRow(*P,i,&ncols,&pcols,&pvals);
866: }
867: MatCreateVecs(A,&diag,0);
868: MatGetDiagonal(A,diag);
869: VecReciprocal(diag);
870: for (i=0;i<cls->nsmooths;i++) {
871: MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew);
872: MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL);
873: MatDiagonalScale(Pnew,diag,0);
874: MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN);
875: MatDestroy(P);
876: *P = Pnew;
877: Pnew = NULL;
878: }
879: VecDestroy(&diag);
880: PetscFree(coarserows);
881: PCGAMGTruncateProlongator_Private(pc,P);
882: return(0);
883: }
885: PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
886: {
887: PetscErrorCode ierr;
888: PetscErrorCode (*f)(PC,Mat,Mat,PetscCoarsenData*,Mat*);
889: PC_MG *mg = (PC_MG*)pc->data;
890: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
891: PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
894: PetscFunctionListFind(PCGAMGClassicalProlongatorList,cls->prolongtype,&f);
895: if (!f)SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot find PCGAMG Classical prolongator type");
896: (*f)(pc,A,G,agg_lists,P);
897: return(0);
898: }
900: PetscErrorCode PCGAMGDestroy_Classical(PC pc)
901: {
903: PC_MG *mg = (PC_MG*)pc->data;
904: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
907: PetscFree(pc_gamg->subctx);
908: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",NULL);
909: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",NULL);
910: return(0);
911: }
913: PetscErrorCode PCGAMGSetFromOptions_Classical(PetscOptionItems *PetscOptionsObject,PC pc)
914: {
915: PC_MG *mg = (PC_MG*)pc->data;
916: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
917: PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
918: char tname[256];
919: PetscErrorCode ierr;
920: PetscBool flg;
923: PetscOptionsHead(PetscOptionsObject,"GAMG-Classical options");
924: PetscOptionsFList("-pc_gamg_classical_type","Type of Classical AMG prolongation","PCGAMGClassicalSetType",PCGAMGClassicalProlongatorList,cls->prolongtype, tname, sizeof(tname), &flg);
925: if (flg) {
926: PCGAMGClassicalSetType(pc,tname);
927: }
928: PetscOptionsReal("-pc_gamg_classical_interp_threshold","Threshold for classical interpolator entries","",cls->interp_threshold,&cls->interp_threshold,NULL);
929: PetscOptionsInt("-pc_gamg_classical_nsmooths","Threshold for classical interpolator entries","",cls->nsmooths,&cls->nsmooths,NULL);
930: PetscOptionsTail();
931: return(0);
932: }
934: PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
935: {
936: PC_MG *mg = (PC_MG*)pc->data;
937: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
940: /* no data for classical AMG */
941: pc_gamg->data = NULL;
942: pc_gamg->data_cell_cols = 0;
943: pc_gamg->data_cell_rows = 0;
944: pc_gamg->data_sz = 0;
945: return(0);
946: }
949: PetscErrorCode PCGAMGClassicalFinalizePackage(void)
950: {
954: PCGAMGClassicalPackageInitialized = PETSC_FALSE;
955: PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList);
956: return(0);
957: }
959: PetscErrorCode PCGAMGClassicalInitializePackage(void)
960: {
964: if (PCGAMGClassicalPackageInitialized) return(0);
965: PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALDIRECT,PCGAMGProlongator_Classical_Direct);
966: PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALSTANDARD,PCGAMGProlongator_Classical_Standard);
967: PetscRegisterFinalize(PCGAMGClassicalFinalizePackage);
968: return(0);
969: }
971: /* -------------------------------------------------------------------------- */
972: /*
973: PCCreateGAMG_Classical
975: */
976: PetscErrorCode PCCreateGAMG_Classical(PC pc)
977: {
979: PC_MG *mg = (PC_MG*)pc->data;
980: PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
981: PC_GAMG_Classical *pc_gamg_classical;
984: PCGAMGClassicalInitializePackage();
985: if (pc_gamg->subctx) {
986: /* call base class */
987: PCDestroy_GAMG(pc);
988: }
990: /* create sub context for SA */
991: PetscNewLog(pc,&pc_gamg_classical);
992: pc_gamg->subctx = pc_gamg_classical;
993: pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
994: /* reset does not do anything; setup not virtual */
996: /* set internal function pointers */
997: pc_gamg->ops->destroy = PCGAMGDestroy_Classical;
998: pc_gamg->ops->graph = PCGAMGGraph_Classical;
999: pc_gamg->ops->coarsen = PCGAMGCoarsen_Classical;
1000: pc_gamg->ops->prolongator = PCGAMGProlongator_Classical;
1001: pc_gamg->ops->optprolongator = PCGAMGOptProlongator_Classical_Jacobi;
1002: pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
1004: pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
1005: pc_gamg_classical->interp_threshold = 0.2;
1006: pc_gamg_classical->nsmooths = 0;
1007: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",PCGAMGClassicalSetType_GAMG);
1008: PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalGetType_C",PCGAMGClassicalGetType_GAMG);
1009: PCGAMGClassicalSetType(pc,PCGAMGCLASSICALSTANDARD);
1010: return(0);
1011: }