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