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