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