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