xref: /petsc/src/ksp/pc/impls/gamg/classical.c (revision 167fb786fdb731ab03cfde3016f1ed5bc75cc9ed) 
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 dummy; /* empty struct; save for later */
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 (PetscAbsScalar(rval[c]) > rmax && rcol[c] != r) {
93         rmax = PetscAbsScalar(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 (PetscAbsScalar(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 (PetscAbsScalar(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   Mat               lG,gG,lA,gA;     /* on and off diagonal matrices */
213   PetscInt          fn;                        /* fine local blocked sizes */
214   PetscInt          cn;                        /* coarse local blocked sizes */
215   PetscInt          gn;                        /* size of the off-diagonal fine vector */
216   PetscInt          fs,fe;                     /* fine (row) ownership range*/
217   PetscInt          cs,ce;                     /* coarse (column) ownership range */
218   PetscInt          i,j,k;                     /* indices! */
219   PetscBool         iscoarse;                  /* flag for determining if a node is coarse */
220   PetscInt          *lcid,*gcid;               /* on and off-processor coarse unknown IDs */
221   PetscInt          *lsparse,*gsparse;         /* on and off-processor sparsity patterns for prolongator */
222   PetscScalar       pij;
223   const PetscScalar *rval;
224   const PetscInt    *rcol;
225   PetscScalar       g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta;
226   Vec               F;   /* vec of coarse size */
227   Vec               C;   /* vec of fine size */
228   Vec               gF;  /* vec of off-diagonal fine size */
229   MatType           mtype;
230   PetscInt          c_indx;
231   const PetscScalar *vcols;
232   const PetscInt    *icols;
233   PetscScalar       c_scalar;
234   PetscInt          ncols,col;
235   PetscInt          row_f,row_c;
236   PetscInt          cmax=0,ncolstotal,idx;
237   PetscScalar       *pvals;
238   PetscInt          *pcols;
239 
240   PetscFunctionBegin;
241   comm = ((PetscObject)pc)->comm;
242   ierr = MatGetOwnershipRange(A,&fs,&fe); CHKERRQ(ierr);
243   fn = (fe - fs);
244 
245   ierr = MatGetVecs(A,&F,NULL);CHKERRQ(ierr);
246 
247   /* get the number of local unknowns and the indices of the local unknowns */
248 
249   ierr = PetscMalloc(sizeof(PetscInt)*fn,&lsparse);CHKERRQ(ierr);
250   ierr = PetscMalloc(sizeof(PetscInt)*fn,&gsparse);CHKERRQ(ierr);
251   ierr = PetscMalloc(sizeof(PetscInt)*fn,&lcid);CHKERRQ(ierr);
252 
253   /* count the number of coarse unknowns */
254   cn = 0;
255   for (i=0;i<fn;i++) {
256     /* filter out singletons */
257     ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr);
258     lcid[i] = -1;
259     if (!iscoarse) {
260       cn++;
261     }
262   }
263 
264    /* create the coarse vector */
265   ierr = VecCreateMPI(comm,cn,PETSC_DECIDE,&C);CHKERRQ(ierr);
266   ierr = VecGetOwnershipRange(C,&cs,&ce);CHKERRQ(ierr);
267 
268   /* construct a global vector indicating the global indices of the coarse unknowns */
269   cn = 0;
270   for (i=0;i<fn;i++) {
271     ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr);
272     if (!iscoarse) {
273       lcid[i] = cs+cn;
274       cn++;
275     } else {
276       lcid[i] = -1;
277     }
278     *((PetscInt *)&c_scalar) = lcid[i];
279     c_indx = fs+i;
280     ierr = VecSetValues(F,1,&c_indx,&c_scalar,INSERT_VALUES);CHKERRQ(ierr);
281   }
282 
283   ierr = VecAssemblyBegin(F);CHKERRQ(ierr);
284   ierr = VecAssemblyEnd(F);CHKERRQ(ierr);
285 
286   /* split the operator into two */
287   ierr = PCGAMGClassicalGraphSplitting_Private(G,&lG,&gG);CHKERRQ(ierr);
288   ierr = PCGAMGClassicalGraphSplitting_Private(A,&lA,&gA);CHKERRQ(ierr);
289 
290   /* scatter to the ghost vector */
291   ierr = PCGAMGClassicalCreateGhostVector_Private(G,&gF,NULL);CHKERRQ(ierr);
292   ierr = PCGAMGClassicalGhost_Private(G,F,gF);CHKERRQ(ierr);
293 
294   if (gG) {
295     ierr = VecGetSize(gF,&gn);CHKERRQ(ierr);
296     ierr = PetscMalloc(sizeof(PetscInt)*gn,&gcid);CHKERRQ(ierr);
297     for (i=0;i<gn;i++) {
298       ierr = VecGetValues(gF,1,&i,&c_scalar);CHKERRQ(ierr);
299       gcid[i] = *((PetscInt *)&c_scalar);
300     }
301   }
302 
303   ierr = VecDestroy(&F);CHKERRQ(ierr);
304   ierr = VecDestroy(&gF);CHKERRQ(ierr);
305   ierr = VecDestroy(&C);CHKERRQ(ierr);
306 
307   /* count the on and off processor sparsity patterns for the prolongator */
308 
309   for (i=0;i<fn;i++) {
310     /* on */
311     ncolstotal = ncols;
312     lsparse[i] = 0;
313     gsparse[i] = 0;
314     if (lcid[i] >= 0) {
315       lsparse[i] = 1;
316       gsparse[i] = 0;
317     } else {
318       ierr = MatGetRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
319       for (j = 0;j < ncols;j++) {
320         col = icols[j];
321         if (lcid[col] >= 0 && vcols[j] != 0.) {
322           lsparse[i] += 1;
323         }
324       }
325       ierr = MatRestoreRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
326       ncolstotal += ncols;
327       /* off */
328       if (gG) {
329         ierr = MatGetRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
330         for (j = 0; j < ncols; j++) {
331           col = icols[j];
332           if (gcid[col] >= 0 && vcols[j] != 0.) {
333             gsparse[i] += 1;
334           }
335         }
336         ierr = MatRestoreRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
337       }
338       if (ncolstotal > cmax) cmax = ncolstotal;
339     }
340   }
341 
342   ierr = PetscMalloc(sizeof(PetscInt)*cmax,&pcols);CHKERRQ(ierr);
343   ierr = PetscMalloc(sizeof(PetscScalar)*cmax,&pvals);CHKERRQ(ierr);
344 
345   /* preallocate and create the prolongator */
346   ierr = MatCreate(comm,P); CHKERRQ(ierr);
347   ierr = MatGetType(G,&mtype);CHKERRQ(ierr);
348   ierr = MatSetType(*P,mtype);CHKERRQ(ierr);
349 
350   ierr = MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
351   ierr = MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);CHKERRQ(ierr);
352   ierr = MatSeqAIJSetPreallocation(*P,0,lsparse);CHKERRQ(ierr);
353 
354   /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
355   for (i = 0;i < fn;i++) {
356     /* determine on or off */
357     row_f = i + fs;
358     row_c = lcid[i];
359     if (row_c >= 0) {
360       pij = 1.;
361       ierr = MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);CHKERRQ(ierr);
362     } else {
363       PetscInt nstrong=0,ntotal=0;
364       g_pos = 0.;
365       g_neg = 0.;
366       a_pos = 0.;
367       a_neg = 0.;
368       diag = 0.;
369 
370       /* local strong connections */
371       ierr = MatGetRow(lG,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
372       for (k = 0; k < ncols; k++) {
373         if (lcid[rcol[k]] >= 0) {
374           if (PetscRealPart(rval[k]) > 0) {
375             g_pos += rval[k];
376           } else {
377             g_neg += rval[k];
378           }
379           nstrong++;
380         }
381       }
382       ierr = MatRestoreRow(lG,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
383 
384       /* ghosted strong connections */
385       if (gG) {
386         ierr = MatGetRow(gG,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
387         for (k = 0; k < ncols; k++) {
388           if (gcid[rcol[k]] >= 0) {
389             if (PetscRealPart(rval[k]) > 0.) {
390               g_pos += rval[k];
391             } else {
392               g_neg += rval[k];
393             }
394             nstrong++;
395           }
396         }
397         ierr = MatRestoreRow(gG,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
398       }
399 
400       /* local all connections */
401       ierr = MatGetRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
402       for (k = 0; k < ncols; k++) {
403         if (rcol[k] != i) {
404           if (PetscRealPart(rval[k]) > 0) {
405             a_pos += rval[k];
406           } else {
407             a_neg += rval[k];
408           }
409           ntotal++;
410         } else diag = rval[k];
411       }
412       ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
413 
414       /* ghosted all connections */
415       if (gA) {
416         ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
417         for (k = 0; k < ncols; k++) {
418           if (PetscRealPart(rval[k]) > 0.) {
419             a_pos += PetscRealPart(rval[k]);
420           } else {
421             a_neg += PetscRealPart(rval[k]);
422           }
423           ntotal++;
424         }
425         ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
426       }
427 
428       if (g_neg == 0.) {
429         alpha = 0.;
430       } else {
431         alpha = -a_neg/g_neg;
432       }
433 
434       if (g_pos == 0.) {
435         diag += a_pos;
436         beta = 0.;
437       } else {
438         beta = -a_pos/g_pos;
439       }
440       if (diag == 0.) {
441         invdiag = 0.;
442       } else invdiag = 1. / diag;
443       /* on */
444       ierr = MatGetRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
445       idx = 0;
446       for (j = 0;j < ncols;j++) {
447         col = icols[j];
448         if (lcid[col] >= 0 && vcols[j] != 0.) {
449           row_f = i + fs;
450           row_c = lcid[col];
451           /* set the values for on-processor ones */
452           if (PetscRealPart(vcols[j]) < 0.) {
453             pij = vcols[j]*alpha*invdiag;
454           } else {
455             pij = vcols[j]*beta*invdiag;
456           }
457           if (PetscAbsScalar(pij) != 0.) {
458             pvals[idx] = pij;
459             pcols[idx] = row_c;
460             idx++;
461           }
462         }
463       }
464       ierr = MatRestoreRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
465       /* off */
466       if (gG) {
467         ierr = MatGetRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
468         for (j = 0; j < ncols; j++) {
469           col = icols[j];
470           if (gcid[col] >= 0 && vcols[j] != 0.) {
471             row_f = i + fs;
472             row_c = gcid[col];
473             /* set the values for on-processor ones */
474             if (PetscRealPart(vcols[j]) < 0.) {
475               pij = vcols[j]*alpha*invdiag;
476             } else {
477               pij = vcols[j]*beta*invdiag;
478             }
479             if (PetscAbsScalar(pij) != 0.) {
480               pvals[idx] = pij;
481               pcols[idx] = row_c;
482               idx++;
483             }
484           }
485         }
486         ierr = MatRestoreRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
487       }
488       ierr = MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);CHKERRQ(ierr);
489     }
490   }
491 
492   ierr = MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
493   ierr = MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
494 
495   ierr = PetscFree(lsparse);CHKERRQ(ierr);
496   ierr = PetscFree(gsparse);CHKERRQ(ierr);
497   ierr = PetscFree(pcols);CHKERRQ(ierr);
498   ierr = PetscFree(pvals);CHKERRQ(ierr);
499   ierr = PetscFree(lcid);CHKERRQ(ierr);
500   if (gG) {ierr = PetscFree(gcid);CHKERRQ(ierr);}
501 
502   PetscFunctionReturn(0);
503 }
504 
505 #undef __FUNCT__
506 #define __FUNCT__ "PCGAMGDestroy_Classical"
507 PetscErrorCode PCGAMGDestroy_Classical(PC pc)
508 {
509   PetscErrorCode ierr;
510   PC_MG          *mg          = (PC_MG*)pc->data;
511   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
512 
513   PetscFunctionBegin;
514   ierr = PetscFree(pc_gamg->subctx);CHKERRQ(ierr);
515   PetscFunctionReturn(0);
516 }
517 
518 #undef __FUNCT__
519 #define __FUNCT__ "PCGAMGSetFromOptions_Classical"
520 PetscErrorCode PCGAMGSetFromOptions_Classical(PC pc)
521 {
522   PetscFunctionBegin;
523   PetscFunctionReturn(0);
524 }
525 
526 #undef __FUNCT__
527 #define __FUNCT__ "PCGAMGSetData_Classical"
528 PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
529 {
530   PC_MG          *mg      = (PC_MG*)pc->data;
531   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
532 
533   PetscFunctionBegin;
534   /* no data for classical AMG */
535   pc_gamg->data           = NULL;
536   pc_gamg->data_cell_cols = 1;
537   pc_gamg->data_cell_rows = 1;
538   pc_gamg->data_sz = 0;
539   PetscFunctionReturn(0);
540 }
541 
542 /* -------------------------------------------------------------------------- */
543 /*
544    PCCreateGAMG_Classical
545 
546 */
547 #undef __FUNCT__
548 #define __FUNCT__ "PCCreateGAMG_Classical"
549 PetscErrorCode  PCCreateGAMG_Classical(PC pc)
550 {
551   PetscErrorCode ierr;
552   PC_MG             *mg      = (PC_MG*)pc->data;
553   PC_GAMG           *pc_gamg = (PC_GAMG*)mg->innerctx;
554   PC_GAMG_Classical *pc_gamg_classical;
555 
556   PetscFunctionBegin;
557   if (pc_gamg->subctx) {
558     /* call base class */
559     ierr = PCDestroy_GAMG(pc);CHKERRQ(ierr);
560   }
561 
562   /* create sub context for SA */
563   ierr = PetscNewLog(pc, PC_GAMG_Classical, &pc_gamg_classical);CHKERRQ(ierr);
564   pc_gamg->subctx = pc_gamg_classical;
565   pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
566   /* reset does not do anything; setup not virtual */
567 
568   /* set internal function pointers */
569   pc_gamg->ops->destroy     = PCGAMGDestroy_Classical;
570   pc_gamg->ops->graph       = PCGAMGGraph_Classical;
571   pc_gamg->ops->coarsen     = PCGAMGCoarsen_Classical;
572   pc_gamg->ops->prolongator = PCGAMGProlongator_Classical;
573   pc_gamg->ops->optprol     = NULL;
574 
575   pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
576   PetscFunctionReturn(0);
577 }
578