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