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