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