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