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