xref: /petsc/src/ksp/pc/impls/gamg/classical.c (revision 87b9b13c6708a27ea1190611c9e657f02c5a6d93) 
1 #include <../src/ksp/pc/impls/gamg/gamg.h>        /*I "petscpc.h" I*/
2 #include <petsc-private/kspimpl.h>
3 #include <petscsf.h>
4 
5 PetscFunctionList PCGAMGClassicalProlongatorList    = NULL;
6 PetscBool         PCGAMGClassicalPackageInitialized = PETSC_FALSE;
7 
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;
13 
14 #undef __FUNCT__
15 #define __FUNCT__ "PCGAMGClassicalSetType"
16 /*@C
17    PCGAMGClassicalSetType - Sets the type of classical interpolation to use
18 
19    Collective on PC
20 
21    Input Parameters:
22 .  pc - the preconditioner context
23 
24    Options Database Key:
25 .  -pc_gamg_classical_type
26 
27    Level: intermediate
28 
29 .seealso: ()
30 @*/
31 PetscErrorCode PCGAMGClassicalSetType(PC pc, PCGAMGClassicalType type)
32 {
33   PetscErrorCode ierr;
34 
35   PetscFunctionBegin;
36   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
37   ierr = PetscTryMethod(pc,"PCGAMGClassicalSetType_C",(PC,PCGAMGType),(pc,type));CHKERRQ(ierr);
38   PetscFunctionReturn(0);
39 }
40 
41 #undef __FUNCT__
42 #define __FUNCT__ "PCGAMGClassicalSetType_GAMG"
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;
49 
50   PetscFunctionBegin;
51   ierr = PetscStrcpy(cls->prolongtype,type);CHKERRQ(ierr);
52   PetscFunctionReturn(0);
53 }
54 
55 #undef __FUNCT__
56 #define __FUNCT__ "PCGAMGGraph_Classical"
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;
73 
74   PetscFunctionBegin;
75   mg   = (PC_MG *)pc->data;
76   gamg = (PC_GAMG *)mg->innerctx;
77 
78   ierr = MatGetOwnershipRange(A,&s,&f);CHKERRQ(ierr);
79   n=f-s;
80   ierr = PetscMalloc1(n,&lsparse);CHKERRQ(ierr);
81   ierr = PetscMalloc1(n,&gsparse);CHKERRQ(ierr);
82   ierr = PetscMalloc1(n,&Amax);CHKERRQ(ierr);
83 
84   for (r = 0;r < n;r++) {
85     lsparse[r] = 0;
86     gsparse[r] = 0;
87   }
88 
89   for (r = s;r < f;r++) {
90     /* determine the maximum off-diagonal in each row */
91     rmax = 0.;
92     ierr = MatGetRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
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     ierr = MatRestoreRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
113     lsparse[r-s] = lidx;
114     gsparse[r-s] = gidx;
115   }
116   ierr = PetscMalloc1(cmax,&gval);CHKERRQ(ierr);
117   ierr = PetscMalloc1(cmax,&gcol);CHKERRQ(ierr);
118 
119   ierr = MatCreate(PetscObjectComm((PetscObject)A),G); CHKERRQ(ierr);
120   ierr = MatGetType(A,&mtype);CHKERRQ(ierr);
121   ierr = MatSetType(*G,mtype);CHKERRQ(ierr);
122   ierr = MatSetSizes(*G,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
123   ierr = MatMPIAIJSetPreallocation(*G,0,lsparse,0,gsparse);CHKERRQ(ierr);
124   ierr = MatSeqAIJSetPreallocation(*G,0,lsparse);CHKERRQ(ierr);
125   for (r = s;r < f;r++) {
126     ierr = MatGetRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
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     ierr = MatSetValues(*G,1,&r,idx,gcol,gval,INSERT_VALUES);CHKERRQ(ierr);
137     ierr = MatRestoreRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
138   }
139   ierr = MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
140   ierr = MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
141 
142   ierr = PetscFree(gval);CHKERRQ(ierr);
143   ierr = PetscFree(gcol);CHKERRQ(ierr);
144   ierr = PetscFree(lsparse);CHKERRQ(ierr);
145   ierr = PetscFree(gsparse);CHKERRQ(ierr);
146   ierr = PetscFree(Amax);CHKERRQ(ierr);
147   PetscFunctionReturn(0);
148 }
149 
150 
151 #undef __FUNCT__
152 #define __FUNCT__ "PCGAMGCoarsen_Classical"
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;
158 
159   PetscFunctionBegin;
160 
161 
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   }
166 
167   ierr = MatCoarsenCreate(fcomm,&crs);CHKERRQ(ierr);
168   ierr = MatCoarsenSetFromOptions(crs);CHKERRQ(ierr);
169   ierr = MatCoarsenSetAdjacency(crs,*G);CHKERRQ(ierr);
170   ierr = MatCoarsenSetStrictAggs(crs,PETSC_TRUE);CHKERRQ(ierr);
171   ierr = MatCoarsenApply(crs);CHKERRQ(ierr);
172   ierr = MatCoarsenGetData(crs,agg_lists);CHKERRQ(ierr);
173   ierr = MatCoarsenDestroy(&crs);CHKERRQ(ierr);
174 
175   PetscFunctionReturn(0);
176 }
177 
178 #undef __FUNCT__
179 #define __FUNCT__ "PCGAMGProlongator_Classical_Direct"
180 PetscErrorCode PCGAMGProlongator_Classical_Direct(PC pc, const Mat A, const Mat G, PetscCoarsenData *agg_lists,Mat *P)
181 {
182   PetscErrorCode    ierr;
183   PetscReal         *Amax_pos,*Amax_neg;
184   Mat               lA,gA=NULL;                /* on and off diagonal matrices */
185   PetscInt          fn;                        /* fine local blocked sizes */
186   PetscInt          cn;                        /* coarse local blocked sizes */
187   PetscInt          fs,fe;                     /* fine (row) ownership range*/
188   PetscInt          cs,ce;                     /* coarse (column) ownership range */
189   PetscInt          i,j;                       /* indices! */
190   PetscBool         iscoarse;                  /* flag for determining if a node is coarse */
191   PetscInt          *lcid,*gcid;               /* on and off-processor coarse unknown IDs */
192   PetscInt          *lsparse,*gsparse;         /* on and off-processor sparsity patterns for prolongator */
193   PetscScalar       pij;
194   const PetscScalar *rval;
195   const PetscInt    *rcol;
196   PetscScalar       g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta;
197   Vec               C;   /* vec of fine size */
198   MatType           mtype;
199   PetscInt          ncols,col;
200   PetscInt          row_f,row_c;
201   PetscInt          cmax=0,idx;
202   PetscScalar       *pvals;
203   PetscInt          *pcols;
204   PC_MG             *mg          = (PC_MG*)pc->data;
205   PC_GAMG           *gamg        = (PC_GAMG*)mg->innerctx;
206   PetscLayout       clayout;
207   PetscSF           sf;
208   Vec               lvec;
209   PetscInt          *colmap,noff;
210   PetscBool         isMPIAIJ,isSEQAIJ;
211   Mat_MPIAIJ        *mpiaij;
212 
213   PetscFunctionBegin;
214   ierr = MatGetOwnershipRange(A,&fs,&fe);CHKERRQ(ierr);
215   fn = fe-fs;
216   ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr);
217   ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSEQAIJ);CHKERRQ(ierr);
218   if (!isMPIAIJ && !isSEQAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Classical AMG requires MPIAIJ matrix");
219   if (isMPIAIJ) {
220     mpiaij = (Mat_MPIAIJ*)A->data;
221     lA = mpiaij->A;
222     gA = mpiaij->B;
223     lvec = mpiaij->lvec;
224     ierr = VecGetSize(lvec,&noff);CHKERRQ(ierr);
225     colmap = mpiaij->garray;
226     ierr = MatGetLayouts(A,NULL,&clayout);CHKERRQ(ierr);
227     ierr = PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);CHKERRQ(ierr);
228     ierr = PetscSFSetGraphLayout(sf,clayout,noff,NULL,PETSC_COPY_VALUES,colmap);CHKERRQ(ierr);
229     ierr = PetscMalloc1(noff,&gcid);CHKERRQ(ierr);
230   } else {
231     lA = A;
232   }
233   ierr = PetscMalloc1(fn,&lsparse);CHKERRQ(ierr);
234   ierr = PetscMalloc1(fn,&gsparse);CHKERRQ(ierr);
235   ierr = PetscMalloc1(fn,&lcid);CHKERRQ(ierr);
236   ierr = PetscMalloc1(fn,&Amax_pos);CHKERRQ(ierr);
237   ierr = PetscMalloc1(fn,&Amax_neg);CHKERRQ(ierr);
238 
239   /* count the number of coarse unknowns */
240   cn = 0;
241   for (i=0;i<fn;i++) {
242     /* filter out singletons */
243     ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr);
244     lcid[i] = -1;
245     if (!iscoarse) {
246       cn++;
247     }
248   }
249 
250    /* create the coarse vector */
251   ierr = VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&C);CHKERRQ(ierr);
252   ierr = VecGetOwnershipRange(C,&cs,&ce);CHKERRQ(ierr);
253 
254   cn = 0;
255   for (i=0;i<fn;i++) {
256     ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr);
257     if (!iscoarse) {
258       lcid[i] = cs+cn;
259       cn++;
260     } else {
261       lcid[i] = -1;
262     }
263   }
264 
265   if (gA) {
266     ierr = PetscSFBcastBegin(sf,MPIU_INT,lcid,gcid);CHKERRQ(ierr);
267     ierr = PetscSFBcastEnd(sf,MPIU_INT,lcid,gcid);CHKERRQ(ierr);
268   }
269 
270   /* determine the biggest off-diagonal entries in each row */
271   for (i=fs;i<fe;i++) {
272     Amax_pos[i-fs] = 0.;
273     Amax_neg[i-fs] = 0.;
274     ierr = MatGetRow(A,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
275     for(j=0;j<ncols;j++){
276       if ((PetscRealPart(-rval[j]) > Amax_neg[i-fs]) && i != rcol[j]) Amax_neg[i-fs] = PetscAbsScalar(rval[j]);
277       if ((PetscRealPart(rval[j])  > Amax_pos[i-fs]) && i != rcol[j]) Amax_pos[i-fs] = PetscAbsScalar(rval[j]);
278     }
279     if (ncols > cmax) cmax = ncols;
280     ierr = MatRestoreRow(A,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
281   }
282   ierr = PetscMalloc1(cmax,&pcols);CHKERRQ(ierr);
283   ierr = PetscMalloc1(cmax,&pvals);CHKERRQ(ierr);
284   ierr = VecDestroy(&C);CHKERRQ(ierr);
285 
286   /* count the on and off processor sparsity patterns for the prolongator */
287   for (i=0;i<fn;i++) {
288     /* on */
289     lsparse[i] = 0;
290     gsparse[i] = 0;
291     if (lcid[i] >= 0) {
292       lsparse[i] = 1;
293       gsparse[i] = 0;
294     } else {
295       ierr = MatGetRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
296       for (j = 0;j < ncols;j++) {
297         col = rcol[j];
298         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
299           lsparse[i] += 1;
300         }
301       }
302       ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
303       /* off */
304       if (gA) {
305         ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
306         for (j = 0; j < ncols; j++) {
307           col = rcol[j];
308           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
309             gsparse[i] += 1;
310           }
311         }
312         ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
313       }
314     }
315   }
316 
317   /* preallocate and create the prolongator */
318   ierr = MatCreate(PetscObjectComm((PetscObject)A),P); CHKERRQ(ierr);
319   ierr = MatGetType(G,&mtype);CHKERRQ(ierr);
320   ierr = MatSetType(*P,mtype);CHKERRQ(ierr);
321 
322   ierr = MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
323   ierr = MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);CHKERRQ(ierr);
324   ierr = MatSeqAIJSetPreallocation(*P,0,lsparse);CHKERRQ(ierr);
325 
326   /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
327   for (i = 0;i < fn;i++) {
328     /* determine on or off */
329     row_f = i + fs;
330     row_c = lcid[i];
331     if (row_c >= 0) {
332       pij = 1.;
333       ierr = MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);CHKERRQ(ierr);
334     } else {
335       g_pos = 0.;
336       g_neg = 0.;
337       a_pos = 0.;
338       a_neg = 0.;
339       diag = 0.;
340 
341       /* local connections */
342       ierr = MatGetRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
343       for (j = 0; j < ncols; j++) {
344         col = rcol[j];
345         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
346           if (PetscRealPart(rval[j]) > 0.) {
347             g_pos += rval[j];
348           } else {
349             g_neg += rval[j];
350           }
351         }
352         if (col != i) {
353           if (PetscRealPart(rval[j]) > 0.) {
354             a_pos += rval[j];
355           } else {
356             a_neg += rval[j];
357           }
358         } else {
359           diag = rval[j];
360         }
361       }
362       ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
363 
364       /* ghosted connections */
365       if (gA) {
366         ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
367         for (j = 0; j < ncols; j++) {
368           col = rcol[j];
369           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
370             if (PetscRealPart(rval[j]) > 0.) {
371               g_pos += rval[j];
372             } else {
373               g_neg += rval[j];
374             }
375           }
376           if (PetscRealPart(rval[j]) > 0.) {
377             a_pos += rval[j];
378           } else {
379             a_neg += rval[j];
380           }
381         }
382         ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
383       }
384 
385       if (g_neg == 0.) {
386         alpha = 0.;
387       } else {
388         alpha = -a_neg/g_neg;
389       }
390 
391       if (g_pos == 0.) {
392         diag += a_pos;
393         beta = 0.;
394       } else {
395         beta = -a_pos/g_pos;
396       }
397       if (diag == 0.) {
398         invdiag = 0.;
399       } else invdiag = 1. / diag;
400       /* on */
401       ierr = MatGetRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
402       idx = 0;
403       for (j = 0;j < ncols;j++) {
404         col = rcol[j];
405         if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
406           row_f = i + fs;
407           row_c = lcid[col];
408           /* set the values for on-processor ones */
409           if (PetscRealPart(rval[j]) < 0.) {
410             pij = rval[j]*alpha*invdiag;
411           } else {
412             pij = rval[j]*beta*invdiag;
413           }
414           if (PetscAbsScalar(pij) != 0.) {
415             pvals[idx] = pij;
416             pcols[idx] = row_c;
417             idx++;
418           }
419         }
420       }
421       ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
422       /* off */
423       if (gA) {
424         ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
425         for (j = 0; j < ncols; j++) {
426           col = rcol[j];
427           if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
428             row_f = i + fs;
429             row_c = gcid[col];
430             /* set the values for on-processor ones */
431             if (PetscRealPart(rval[j]) < 0.) {
432               pij = rval[j]*alpha*invdiag;
433             } else {
434               pij = rval[j]*beta*invdiag;
435             }
436             if (PetscAbsScalar(pij) != 0.) {
437               pvals[idx] = pij;
438               pcols[idx] = row_c;
439               idx++;
440             }
441           }
442         }
443         ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
444       }
445       ierr = MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);CHKERRQ(ierr);
446     }
447   }
448 
449   ierr = MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
450   ierr = MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
451 
452   ierr = PetscFree(lsparse);CHKERRQ(ierr);
453   ierr = PetscFree(gsparse);CHKERRQ(ierr);
454   ierr = PetscFree(pcols);CHKERRQ(ierr);
455   ierr = PetscFree(pvals);CHKERRQ(ierr);
456   ierr = PetscFree(Amax_pos);CHKERRQ(ierr);
457   ierr = PetscFree(Amax_neg);CHKERRQ(ierr);
458   ierr = PetscFree(lcid);CHKERRQ(ierr);
459   if (gA) {
460     ierr = PetscSFDestroy(&sf);CHKERRQ(ierr);
461     ierr = PetscFree(gcid);CHKERRQ(ierr);
462   }
463 
464   PetscFunctionReturn(0);
465 }
466 
467 #undef __FUNCT__
468 #define __FUNCT__ "PCGAMGTruncateProlongator_Private"
469 PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc,Mat *P)
470 {
471   PetscInt          j,i,ps,pf,pn,pcs,pcf,pcn,idx,cmax;
472   PetscErrorCode    ierr;
473   const PetscScalar *pval;
474   const PetscInt    *pcol;
475   PetscScalar       *pnval;
476   PetscInt          *pncol;
477   PetscInt          ncols;
478   Mat               Pnew;
479   PetscInt          *lsparse,*gsparse;
480   PetscReal         pmax_pos,pmax_neg,ptot_pos,ptot_neg,pthresh_pos,pthresh_neg;
481   PC_MG             *mg          = (PC_MG*)pc->data;
482   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
483   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
484 
485   PetscFunctionBegin;
486   /* trim and rescale with reallocation */
487   ierr = MatGetOwnershipRange(*P,&ps,&pf);CHKERRQ(ierr);
488   ierr = MatGetOwnershipRangeColumn(*P,&pcs,&pcf);CHKERRQ(ierr);
489   pn = pf-ps;
490   pcn = pcf-pcs;
491   ierr = PetscMalloc1(pn,&lsparse);CHKERRQ(ierr);
492   ierr = PetscMalloc1(pn,&gsparse);CHKERRQ(ierr);
493   /* allocate */
494   cmax = 0;
495   for (i=ps;i<pf;i++) {
496     lsparse[i-ps] = 0;
497     gsparse[i-ps] = 0;
498     ierr = MatGetRow(*P,i,&ncols,&pcol,&pval);CHKERRQ(ierr);
499     if (ncols > cmax) {
500       cmax = ncols;
501     }
502     pmax_pos = 0.;
503     pmax_neg = 0.;
504     for (j=0;j<ncols;j++) {
505       if (PetscRealPart(pval[j]) > pmax_pos) {
506         pmax_pos = PetscRealPart(pval[j]);
507       } else if (PetscRealPart(pval[j]) < pmax_neg) {
508         pmax_neg = PetscRealPart(pval[j]);
509       }
510     }
511     for (j=0;j<ncols;j++) {
512       if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
513         if (pcol[j] >= pcs && pcol[j] < pcf) {
514           lsparse[i-ps]++;
515         } else {
516           gsparse[i-ps]++;
517         }
518       }
519     }
520     ierr = MatRestoreRow(*P,i,&ncols,&pcol,&pval);CHKERRQ(ierr);
521   }
522 
523   ierr = PetscMalloc1(cmax,&pnval);CHKERRQ(ierr);
524   ierr = PetscMalloc1(cmax,&pncol);CHKERRQ(ierr);
525 
526   ierr = MatCreate(PetscObjectComm((PetscObject)*P),&Pnew);CHKERRQ(ierr);
527   ierr = MatSetType(Pnew, MATAIJ);CHKERRQ(ierr);
528   ierr = MatSetSizes(Pnew,pn,pcn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
529   ierr = MatSeqAIJSetPreallocation(Pnew,0,lsparse);CHKERRQ(ierr);
530   ierr = MatMPIAIJSetPreallocation(Pnew,0,lsparse,0,gsparse);CHKERRQ(ierr);
531 
532   for (i=ps;i<pf;i++) {
533     ierr = MatGetRow(*P,i,&ncols,&pcol,&pval);CHKERRQ(ierr);
534     pmax_pos = 0.;
535     pmax_neg = 0.;
536     for (j=0;j<ncols;j++) {
537       if (PetscRealPart(pval[j]) > pmax_pos) {
538         pmax_pos = PetscRealPart(pval[j]);
539       } else if (PetscRealPart(pval[j]) < pmax_neg) {
540         pmax_neg = PetscRealPart(pval[j]);
541       }
542     }
543     pthresh_pos = 0.;
544     pthresh_neg = 0.;
545     ptot_pos = 0.;
546     ptot_neg = 0.;
547     for (j=0;j<ncols;j++) {
548       if (PetscRealPart(pval[j]) >= cls->interp_threshold*pmax_pos) {
549         pthresh_pos += PetscRealPart(pval[j]);
550       } else if (PetscRealPart(pval[j]) <= cls->interp_threshold*pmax_neg) {
551         pthresh_neg += PetscRealPart(pval[j]);
552       }
553       if (PetscRealPart(pval[j]) > 0.) {
554         ptot_pos += PetscRealPart(pval[j]);
555       } else {
556         ptot_neg += PetscRealPart(pval[j]);
557       }
558     }
559     if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
560     if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
561     idx=0;
562     for (j=0;j<ncols;j++) {
563       if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold) {
564         pnval[idx] = ptot_pos*pval[j];
565         pncol[idx] = pcol[j];
566         idx++;
567       } else if (PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
568         pnval[idx] = ptot_neg*pval[j];
569         pncol[idx] = pcol[j];
570         idx++;
571       }
572     }
573     ierr = MatRestoreRow(*P,i,&ncols,&pcol,&pval);CHKERRQ(ierr);
574     ierr = MatSetValues(Pnew,1,&i,idx,pncol,pnval,INSERT_VALUES);CHKERRQ(ierr);
575   }
576 
577   ierr = MatAssemblyBegin(Pnew, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
578   ierr = MatAssemblyEnd(Pnew, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
579   ierr = MatDestroy(P);CHKERRQ(ierr);
580 
581   *P = Pnew;
582   ierr = PetscFree(lsparse);CHKERRQ(ierr);
583   ierr = PetscFree(gsparse);CHKERRQ(ierr);
584   ierr = PetscFree(pncol);CHKERRQ(ierr);
585   ierr = PetscFree(pnval);CHKERRQ(ierr);
586   PetscFunctionReturn(0);
587 }
588 
589 #undef __FUNCT__
590 #define __FUNCT__ "PCGAMGProlongator_Classical_Standard"
591 PetscErrorCode PCGAMGProlongator_Classical_Standard(PC pc, const Mat A, const Mat G, PetscCoarsenData *agg_lists,Mat *P)
592 {
593   PetscErrorCode    ierr;
594   Mat               *lA;
595   Vec               lv,v,cv;
596   PetscScalar       *lcid;
597   IS                lis;
598   PetscInt          fs,fe,cs,ce,nl,i,j,k,li,lni,ci;
599   VecScatter        lscat;
600   PetscInt          fn,cn,cid,c_indx;
601   PetscBool         iscoarse;
602   PetscScalar       c_scalar;
603   const PetscScalar *vcol;
604   const PetscInt    *icol;
605   const PetscInt    *gidx;
606   PetscInt          ncols;
607   PetscInt          *lsparse,*gsparse;
608   MatType           mtype;
609   PetscInt          maxcols;
610   PetscReal         diag,jdiag,jwttotal;
611   PetscScalar       *pvcol,vi;
612   PetscInt          *picol;
613   PetscInt          pncols;
614   PetscScalar       *pcontrib,pentry,pjentry;
615   /* PC_MG             *mg          = (PC_MG*)pc->data; */
616   /* PC_GAMG           *gamg        = (PC_GAMG*)mg->innerctx; */
617 
618   PetscFunctionBegin;
619 
620   ierr = MatGetOwnershipRange(A,&fs,&fe);CHKERRQ(ierr);
621   fn = fe-fs;
622   ierr = MatGetVecs(A,NULL,&v);CHKERRQ(ierr);
623   ierr = ISCreateStride(PETSC_COMM_SELF,fe-fs,fs,1,&lis);CHKERRQ(ierr);
624   /* increase the overlap by two to get neighbors of neighbors */
625   ierr = MatIncreaseOverlap(A,1,&lis,2);CHKERRQ(ierr);
626   ierr = ISSort(lis);CHKERRQ(ierr);
627   /* create a communication structure for the overlapped portion and transmit coarse indices */
628 
629   /* get the local part of A */
630   ierr = MatGetSubMatrices(A,1,&lis,&lis,MAT_INITIAL_MATRIX,&lA);CHKERRQ(ierr);
631   /* build the scatter out of it */
632   ierr = ISGetLocalSize(lis,&nl);CHKERRQ(ierr);
633   ierr = VecCreateSeq(PETSC_COMM_SELF,nl,&lv);CHKERRQ(ierr);
634   ierr = VecScatterCreate(v,lis,lv,NULL,&lscat);CHKERRQ(ierr);
635 
636   ierr = PetscMalloc1(fn,&lsparse);CHKERRQ(ierr);
637   ierr = PetscMalloc1(fn,&gsparse);CHKERRQ(ierr);
638   ierr = PetscMalloc1(nl,&pcontrib);CHKERRQ(ierr);
639 
640   /* create coarse vector */
641   cn = 0;
642   for (i=0;i<fn;i++) {
643     ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse);CHKERRQ(ierr);
644     if (!iscoarse) {
645       cn++;
646     }
647   }
648   ierr = VecCreateMPI(PetscObjectComm((PetscObject)A),cn,PETSC_DECIDE,&cv);CHKERRQ(ierr);
649   ierr = VecGetOwnershipRange(cv,&cs,&ce);CHKERRQ(ierr);
650   cn = 0;
651   for (i=0;i<fn;i++) {
652     ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr);
653     if (!iscoarse) {
654       cid = cs+cn;
655       cn++;
656     } else {
657       cid = -1;
658     }
659     *(PetscInt*)&c_scalar = cid;
660     c_indx = fs+i;
661     ierr = VecSetValues(v,1,&c_indx,&c_scalar,INSERT_VALUES);CHKERRQ(ierr);
662   }
663   ierr = VecScatterBegin(lscat,v,lv,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
664   ierr = VecScatterEnd(lscat,v,lv,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
665   /* count to preallocate the prolongator */
666   ierr = ISGetIndices(lis,&gidx);CHKERRQ(ierr);
667   ierr = VecGetArray(lv,&lcid);CHKERRQ(ierr);
668   maxcols = 0;
669   /* count the number of unique contributing coarse cells for each fine */
670   for (i=0;i<nl;i++) {
671     pcontrib[i] = 0.;
672     ierr = MatGetRow(lA[0],i,&ncols,&icol,NULL);CHKERRQ(ierr);
673     if (gidx[i] >= fs && gidx[i] < fe) {
674       li = gidx[i] - fs;
675       lsparse[li] = 0;
676       gsparse[li] = 0;
677       cid = *(PetscInt*)&(lcid[i]);
678       if (cid >= 0) {
679         lsparse[li] = 1;
680       } else {
681         for (j=0;j<ncols;j++) {
682           if (*(PetscInt*)&(lcid[icol[j]]) >= 0) {
683             pcontrib[icol[j]] = 1.;
684           } else {
685             ci = icol[j];
686             ierr = MatRestoreRow(lA[0],i,&ncols,&icol,NULL);CHKERRQ(ierr);
687             ierr = MatGetRow(lA[0],ci,&ncols,&icol,NULL);CHKERRQ(ierr);
688             for (k=0;k<ncols;k++) {
689               if (*(PetscInt*)&(lcid[icol[k]]) >= 0) {
690                 pcontrib[icol[k]] = 1.;
691               }
692             }
693             ierr = MatRestoreRow(lA[0],ci,&ncols,&icol,NULL);CHKERRQ(ierr);
694             ierr = MatGetRow(lA[0],i,&ncols,&icol,NULL);CHKERRQ(ierr);
695           }
696         }
697         for (j=0;j<ncols;j++) {
698           if (*(PetscInt*)&(lcid[icol[j]]) >= 0 && pcontrib[icol[j]] != 0.) {
699             lni = *(PetscInt*)&(lcid[icol[j]]);
700             if (lni >= cs && lni < ce) {
701               lsparse[li]++;
702             } else {
703               gsparse[li]++;
704             }
705             pcontrib[icol[j]] = 0.;
706           } else {
707             ci = icol[j];
708             ierr = MatRestoreRow(lA[0],i,&ncols,&icol,NULL);CHKERRQ(ierr);
709             ierr = MatGetRow(lA[0],ci,&ncols,&icol,NULL);CHKERRQ(ierr);
710             for (k=0;k<ncols;k++) {
711               if (*(PetscInt*)&(lcid[icol[k]]) >= 0 && pcontrib[icol[k]] != 0.) {
712                 lni = *(PetscInt*)&(lcid[icol[k]]);
713                 if (lni >= cs && lni < ce) {
714                   lsparse[li]++;
715                 } else {
716                   gsparse[li]++;
717                 }
718                 pcontrib[icol[k]] = 0.;
719               }
720             }
721             ierr = MatRestoreRow(lA[0],ci,&ncols,&icol,NULL);CHKERRQ(ierr);
722             ierr = MatGetRow(lA[0],i,&ncols,&icol,NULL);CHKERRQ(ierr);
723           }
724         }
725       }
726       if (lsparse[li] + gsparse[li] > maxcols) maxcols = lsparse[li]+gsparse[li];
727     }
728     ierr = MatRestoreRow(lA[0],i,&ncols,&icol,&vcol);CHKERRQ(ierr);
729   }
730   ierr = PetscMalloc1(maxcols,&picol);CHKERRQ(ierr);
731   ierr = PetscMalloc1(maxcols,&pvcol);CHKERRQ(ierr);
732   ierr = MatCreate(PetscObjectComm((PetscObject)A),P);CHKERRQ(ierr);
733   ierr = MatGetType(A,&mtype);CHKERRQ(ierr);
734   ierr = MatSetType(*P,mtype);CHKERRQ(ierr);
735   ierr = MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
736   ierr = MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);CHKERRQ(ierr);
737   ierr = MatSeqAIJSetPreallocation(*P,0,lsparse);CHKERRQ(ierr);
738   for (i=0;i<nl;i++) {
739     diag = 0.;
740     if (gidx[i] >= fs && gidx[i] < fe) {
741       li = gidx[i] - fs;
742       pncols=0;
743       cid = *(PetscInt*)&(lcid[i]);
744       if (cid >= 0) {
745         pncols = 1;
746         picol[0] = cid;
747         pvcol[0] = 1.;
748       } else {
749         ierr = MatGetRow(lA[0],i,&ncols,&icol,&vcol);CHKERRQ(ierr);
750         for (j=0;j<ncols;j++) {
751           pentry = vcol[j];
752           if (*(PetscInt*)&(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             ierr = MatRestoreRow(lA[0],i,&ncols,&icol,&vcol);CHKERRQ(ierr);
760             ierr = MatGetRow(lA[0],ci,&ncols,&icol,&vcol);CHKERRQ(ierr);
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 (*(PetscInt*)&(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 (*(PetscInt*)&(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             ierr = MatRestoreRow(lA[0],ci,&ncols,&icol,&vcol);CHKERRQ(ierr);
786             ierr = MatGetRow(lA[0],i,&ncols,&icol,&vcol);CHKERRQ(ierr);
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 (*(PetscInt*)&(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 = *(PetscInt*)&(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               ierr = MatRestoreRow(lA[0],i,&ncols,&icol,&vcol);CHKERRQ(ierr);
807               ierr = MatGetRow(lA[0],ci,&ncols,&icol,&vcol);CHKERRQ(ierr);
808               for (k=0;k<ncols;k++) {
809                 if (*(PetscInt*)&(lcid[icol[k]]) >= 0 && pcontrib[icol[k]] != 0.) {
810                   if (PetscAbsScalar(pcontrib[icol[k]]) > 0.0) {
811                     lni = *(PetscInt*)&(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               ierr = MatRestoreRow(lA[0],ci,&ncols,&icol,&vcol);CHKERRQ(ierr);
820               ierr = MatGetRow(lA[0],i,&ncols,&icol,&vcol);CHKERRQ(ierr);
821             }
822             pcontrib[icol[j]] = 0.;
823           }
824           ierr = MatRestoreRow(lA[0],i,&ncols,&icol,&vcol);CHKERRQ(ierr);
825         }
826       }
827       ci = gidx[i];
828       li = gidx[i] - fs;
829       if (pncols > 0) {
830         ierr = MatSetValues(*P,1,&ci,pncols,picol,pvcol,INSERT_VALUES);CHKERRQ(ierr);
831       }
832     }
833   }
834   ierr = ISRestoreIndices(lis,&gidx);CHKERRQ(ierr);
835   ierr = VecRestoreArray(lv,&lcid);CHKERRQ(ierr);
836 
837   ierr = PetscFree(pcontrib);CHKERRQ(ierr);
838   ierr = PetscFree(picol);CHKERRQ(ierr);
839   ierr = PetscFree(pvcol);CHKERRQ(ierr);
840   ierr = PetscFree(lsparse);CHKERRQ(ierr);
841   ierr = PetscFree(gsparse);CHKERRQ(ierr);
842   ierr = ISDestroy(&lis);CHKERRQ(ierr);
843   ierr = MatDestroyMatrices(1,&lA);CHKERRQ(ierr);
844   ierr = VecDestroy(&lv);CHKERRQ(ierr);
845   ierr = VecDestroy(&cv);CHKERRQ(ierr);
846   ierr = VecDestroy(&v);CHKERRQ(ierr);
847   ierr = VecScatterDestroy(&lscat);CHKERRQ(ierr);
848 
849   ierr = MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
850   ierr = MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
851 
852   /*
853   Mat Pold;
854   ierr = PCGAMGProlongator_Classical(pc,A,G,agg_lists,&Pold);CHKERRQ(ierr);
855   ierr = MatView(Pold,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
856   ierr = MatView(*P,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
857   ierr = MatDestroy(&Pold);CHKERRQ(ierr);
858    */
859   PetscFunctionReturn(0);
860 }
861 
862 #undef __FUNCT__
863 #define __FUNCT__ "PCGAMGOptProl_Classical_Jacobi"
864 PetscErrorCode PCGAMGOptProl_Classical_Jacobi(PC pc,const Mat A,Mat *P)
865 {
866 
867   PetscErrorCode    ierr;
868   PetscInt          f,s,n,cf,cs,i,idx;
869   PetscInt          *coarserows;
870   PetscInt          ncols;
871   const PetscInt    *pcols;
872   const PetscScalar *pvals;
873   Mat               Pnew;
874   Vec               diag;
875   PC_MG             *mg          = (PC_MG*)pc->data;
876   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
877   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
878 
879   PetscFunctionBegin;
880   if (cls->nsmooths == 0) {
881     ierr = PCGAMGTruncateProlongator_Private(pc,P);CHKERRQ(ierr);
882     PetscFunctionReturn(0);
883   }
884   ierr = MatGetOwnershipRange(*P,&s,&f);CHKERRQ(ierr);
885   n = f-s;
886   ierr = MatGetOwnershipRangeColumn(*P,&cs,&cf);CHKERRQ(ierr);
887   ierr = PetscMalloc1(n,&coarserows);CHKERRQ(ierr);
888   /* identify the rows corresponding to coarse unknowns */
889   idx = 0;
890   for (i=s;i<f;i++) {
891     ierr = MatGetRow(*P,i,&ncols,&pcols,&pvals);CHKERRQ(ierr);
892     /* assume, for now, that it's a coarse unknown if it has a single unit entry */
893     if (ncols == 1) {
894       if (pvals[0] == 1.) {
895         coarserows[idx] = i;
896         idx++;
897       }
898     }
899     ierr = MatRestoreRow(*P,i,&ncols,&pcols,&pvals);CHKERRQ(ierr);
900   }
901   ierr = MatGetVecs(A,&diag,0);CHKERRQ(ierr);
902   ierr = MatGetDiagonal(A,diag);CHKERRQ(ierr);
903   ierr = VecReciprocal(diag);CHKERRQ(ierr);
904   for (i=0;i<cls->nsmooths;i++) {
905     ierr = MatMatMult(A,*P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Pnew);CHKERRQ(ierr);
906     ierr = MatZeroRows(Pnew,idx,coarserows,0.,NULL,NULL);CHKERRQ(ierr);
907     ierr = MatDiagonalScale(Pnew,diag,0);CHKERRQ(ierr);
908     ierr = MatAYPX(Pnew,-1.0,*P,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
909     ierr = MatDestroy(P);CHKERRQ(ierr);
910     *P  = Pnew;
911     Pnew = NULL;
912   }
913   ierr = VecDestroy(&diag);CHKERRQ(ierr);
914   ierr = PetscFree(coarserows);CHKERRQ(ierr);
915   ierr = PCGAMGTruncateProlongator_Private(pc,P);CHKERRQ(ierr);
916   PetscFunctionReturn(0);
917 }
918 
919 #undef __FUNCT__
920 #define __FUNCT__ "PCGAMGProlongator_Classical"
921 PetscErrorCode PCGAMGProlongator_Classical(PC pc, const Mat A, const Mat G, PetscCoarsenData *agg_lists,Mat *P)
922 {
923   PetscErrorCode    ierr;
924   PetscErrorCode    (*f)(PC,Mat,Mat,PetscCoarsenData*,Mat*);
925   PC_MG             *mg          = (PC_MG*)pc->data;
926   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
927   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
928 
929   PetscFunctionBegin;
930   ierr = PetscFunctionListFind(PCGAMGClassicalProlongatorList,cls->prolongtype,&f);CHKERRQ(ierr);
931   if (!f)SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Cannot find PCGAMG Classical prolongator type");
932   ierr = (*f)(pc,A,G,agg_lists,P);CHKERRQ(ierr);
933   PetscFunctionReturn(0);
934 }
935 
936 #undef __FUNCT__
937 #define __FUNCT__ "PCGAMGDestroy_Classical"
938 PetscErrorCode PCGAMGDestroy_Classical(PC pc)
939 {
940   PetscErrorCode ierr;
941   PC_MG          *mg          = (PC_MG*)pc->data;
942   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
943 
944   PetscFunctionBegin;
945   ierr = PetscFree(pc_gamg->subctx);CHKERRQ(ierr);
946   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",NULL);CHKERRQ(ierr);
947   PetscFunctionReturn(0);
948 }
949 
950 #undef __FUNCT__
951 #define __FUNCT__ "PCGAMGSetFromOptions_Classical"
952 PetscErrorCode PCGAMGSetFromOptions_Classical(PC pc)
953 {
954   PC_MG             *mg          = (PC_MG*)pc->data;
955   PC_GAMG           *pc_gamg     = (PC_GAMG*)mg->innerctx;
956   PC_GAMG_Classical *cls         = (PC_GAMG_Classical*)pc_gamg->subctx;
957   char              tname[256];
958   PetscErrorCode    ierr;
959   PetscBool         flg;
960 
961   PetscFunctionBegin;
962   ierr = PetscOptionsHead("GAMG-Classical options");CHKERRQ(ierr);
963   ierr = PetscOptionsFList("-pc_gamg_classical_type","Type of Classical AMG prolongation",
964                           "PCGAMGClassicalSetType",PCGAMGClassicalProlongatorList,cls->prolongtype, tname, sizeof(tname), &flg);CHKERRQ(ierr);
965   if (flg) {
966     ierr = PCGAMGClassicalSetType(pc,tname);CHKERRQ(ierr);
967   }
968   ierr = PetscOptionsReal("-pc_gamg_classical_interp_threshold","Threshold for classical interpolator entries","",cls->interp_threshold,&cls->interp_threshold,NULL);CHKERRQ(ierr);
969   ierr = PetscOptionsInt("-pc_gamg_classical_nsmooths","Threshold for classical interpolator entries","",cls->nsmooths,&cls->nsmooths,NULL);CHKERRQ(ierr);
970   ierr = PetscOptionsTail();CHKERRQ(ierr);
971   PetscFunctionReturn(0);
972 }
973 
974 #undef __FUNCT__
975 #define __FUNCT__ "PCGAMGSetData_Classical"
976 PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
977 {
978   PC_MG          *mg      = (PC_MG*)pc->data;
979   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
980 
981   PetscFunctionBegin;
982   /* no data for classical AMG */
983   pc_gamg->data = NULL;
984   pc_gamg->data_cell_cols = 0;
985   pc_gamg->data_cell_rows = 0;
986   pc_gamg->data_sz        = 0;
987   PetscFunctionReturn(0);
988 }
989 
990 
991 #undef __FUNCT__
992 #define __FUNCT__ "PCGAMGClassicalFinalizePackage"
993 PetscErrorCode PCGAMGClassicalFinalizePackage(void)
994 {
995   PetscErrorCode ierr;
996 
997   PetscFunctionBegin;
998   PCGAMGClassicalPackageInitialized = PETSC_FALSE;
999   ierr = PetscFunctionListDestroy(&PCGAMGClassicalProlongatorList);CHKERRQ(ierr);
1000   PetscFunctionReturn(0);
1001 }
1002 
1003 #undef __FUNCT__
1004 #define __FUNCT__ "PCGAMGClassicalInitializePackage"
1005 PetscErrorCode PCGAMGClassicalInitializePackage(void)
1006 {
1007   PetscErrorCode ierr;
1008 
1009   PetscFunctionBegin;
1010   if (PCGAMGClassicalPackageInitialized) PetscFunctionReturn(0);
1011   ierr = PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALDIRECT,PCGAMGProlongator_Classical_Direct);CHKERRQ(ierr);
1012   ierr = PetscFunctionListAdd(&PCGAMGClassicalProlongatorList,PCGAMGCLASSICALSTANDARD,PCGAMGProlongator_Classical_Standard);CHKERRQ(ierr);
1013   ierr = PetscRegisterFinalize(PCGAMGClassicalFinalizePackage);CHKERRQ(ierr);
1014   PetscFunctionReturn(0);
1015 }
1016 
1017 /* -------------------------------------------------------------------------- */
1018 /*
1019    PCCreateGAMG_Classical
1020 
1021 */
1022 #undef __FUNCT__
1023 #define __FUNCT__ "PCCreateGAMG_Classical"
1024 PetscErrorCode  PCCreateGAMG_Classical(PC pc)
1025 {
1026   PetscErrorCode ierr;
1027   PC_MG             *mg      = (PC_MG*)pc->data;
1028   PC_GAMG           *pc_gamg = (PC_GAMG*)mg->innerctx;
1029   PC_GAMG_Classical *pc_gamg_classical;
1030 
1031   PetscFunctionBegin;
1032   ierr = PCGAMGClassicalInitializePackage();
1033   if (pc_gamg->subctx) {
1034     /* call base class */
1035     ierr = PCDestroy_GAMG(pc);CHKERRQ(ierr);
1036   }
1037 
1038   /* create sub context for SA */
1039   ierr = PetscNewLog(pc,&pc_gamg_classical);CHKERRQ(ierr);
1040   pc_gamg->subctx = pc_gamg_classical;
1041   pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
1042   /* reset does not do anything; setup not virtual */
1043 
1044   /* set internal function pointers */
1045   pc_gamg->ops->destroy        = PCGAMGDestroy_Classical;
1046   pc_gamg->ops->graph          = PCGAMGGraph_Classical;
1047   pc_gamg->ops->coarsen        = PCGAMGCoarsen_Classical;
1048   pc_gamg->ops->prolongator    = PCGAMGProlongator_Classical;
1049   pc_gamg->ops->optprol        = PCGAMGOptProl_Classical_Jacobi;
1050   pc_gamg->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
1051 
1052   pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
1053   pc_gamg_classical->interp_threshold = 0.2;
1054   pc_gamg_classical->nsmooths         = 0;
1055   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGClassicalSetType_C",PCGAMGClassicalSetType_GAMG);CHKERRQ(ierr);
1056   ierr = PCGAMGClassicalSetType(pc,PCGAMGCLASSICALSTANDARD);CHKERRQ(ierr);
1057   PetscFunctionReturn(0);
1058 }
1059