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