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