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