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