xref: /petsc/src/ksp/pc/impls/gamg/classical.c (revision b375234f3613e90cb42b059dc2aa19764ec63904)
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,Mat A,Mat *G)
56 {
57   PetscInt          s,f,idx;
58   PetscInt          r,c,ncols,*gcol;
59   const PetscInt    *rcol;
60   const PetscScalar *rval;
61   PetscScalar       *gval;
62   PetscReal         rmax,Amax;
63   PetscInt          cmax;
64   PC_MG             *mg;
65   PC_GAMG           *gamg;
66   PetscErrorCode    ierr;
67   PetscInt          *gsparse,*lsparse;
68   Mat               lA,gA;
69   MatType           mtype;
70 
71   PetscFunctionBegin;
72   mg   = (PC_MG *)pc->data;
73   gamg = (PC_GAMG *)mg->innerctx;
74 
75   ierr = MatGetOwnershipRange(A,&s,&f);CHKERRQ(ierr);
76 
77   ierr = PCGAMGClassicalGraphSplitting_Private(A,&lA,&gA);CHKERRQ(ierr);
78 
79   ierr = PetscMalloc(sizeof(PetscInt)*(f - s),&lsparse);CHKERRQ(ierr);
80   if (gA) {ierr = PetscMalloc(sizeof(PetscInt)*(f - s),&gsparse);CHKERRQ(ierr);}
81   else {
82     gsparse = NULL;
83   }
84 
85   /* find the maximum off-diagonal entry in the matrix */
86   rmax = 0.;
87   cmax = 0;
88   for (r = s;r < f;r++) {
89     ierr = MatGetRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
90     for (c = 0; c < ncols; c++) {
91       if (rcol[c] != r)
92         if (PetscAbsScalar(rval[c]) > rmax) rmax = PetscAbsScalar(rval[c]);
93     }
94     if (ncols > cmax) cmax = ncols;
95     ierr = MatRestoreRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
96   }
97   ierr = PetscMalloc(sizeof(PetscScalar)*cmax,&gval);CHKERRQ(ierr);
98   ierr = PetscMalloc(sizeof(PetscInt)*cmax,&gcol);CHKERRQ(ierr);
99   ierr = MPI_Allreduce(&rmax,&Amax,1,MPI_DOUBLE,MPI_MAX,PetscObjectComm((PetscObject)pc));
100 
101   ierr = PetscInfo2(pc,"Maximum off-diagonal value in classical AMG graph: %f threshold: %f \n",rmax,gamg->threshold);CHKERRQ(ierr);
102 
103   for (r = 0;r < f-s;r++) {
104     lsparse[r] = 0;
105     if (gsparse) gsparse[r] = 0;
106   }
107 
108   /* for now this recreates the entire matrix due to a bug in MatCoarsen */
109   for (r = 0;r < f-s;r++) {
110     ierr = MatGetRow(lA,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
111     idx = 0;
112     for (c = 0; c < ncols; c++) {
113       if (PetscAbsScalar(rval[c]) > gamg->threshold*Amax) {
114         idx++;
115       }
116     }
117     ierr = MatRestoreRow(lA,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
118     lsparse[r] = idx;
119   }
120   if (gA) {
121     for (r = 0;r < f-s;r++) {
122       ierr = MatGetRow(gA,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
123       idx = 0;
124       for (c = 0; c < ncols; c++) {
125         if (PetscAbsScalar(rval[c]) > gamg->threshold*Amax) {
126           idx++;
127         }
128       }
129       ierr = MatRestoreRow(gA,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
130       gsparse[r] = idx;
131     }
132   }
133   ierr = MatCreate(PetscObjectComm((PetscObject)A),G); CHKERRQ(ierr);
134   ierr = MatGetType(A,&mtype);CHKERRQ(ierr);
135   ierr = MatSetType(*G,mtype);CHKERRQ(ierr);
136   ierr = MatSetSizes(*G,f-s,f-s,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
137   ierr = MatMPIAIJSetPreallocation(*G,0,lsparse,0,gsparse);CHKERRQ(ierr);
138   ierr = MatSeqAIJSetPreallocation(*G,0,lsparse);CHKERRQ(ierr);
139   for (r = s;r < f;r++) {
140     ierr = MatGetRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
141     idx = 0;
142     for (c = 0; c < ncols; c++) {
143       /* classical strength of connection */
144       if (PetscAbsScalar(rval[c]) > gamg->threshold*Amax) {
145         gcol[idx] = rcol[c];
146         gval[idx] = rval[c];
147         idx++;
148       }
149     }
150     ierr = MatSetValues(*G,1,&r,idx,gcol,gval,INSERT_VALUES);CHKERRQ(ierr);
151     ierr = MatRestoreRow(A,r,&ncols,&rcol,&rval);CHKERRQ(ierr);
152   }
153   ierr = MatAssemblyBegin(*G, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
154   ierr = MatAssemblyEnd(*G, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
155 
156   ierr = PetscFree(gval);CHKERRQ(ierr);
157   ierr = PetscFree(gcol);CHKERRQ(ierr);
158   ierr = PetscFree(lsparse);CHKERRQ(ierr);
159   ierr = PetscFree(gsparse);CHKERRQ(ierr);
160   PetscFunctionReturn(0);
161 }
162 
163 
164 #undef __FUNCT__
165 #define __FUNCT__ "PCGAMGCoarsen_Classical"
166 PetscErrorCode PCGAMGCoarsen_Classical(PC pc,Mat *G,PetscCoarsenData **agg_lists)
167 {
168   PetscErrorCode   ierr;
169   MatCoarsen       crs;
170   MPI_Comm         fcomm = ((PetscObject)pc)->comm;
171 
172   PetscFunctionBegin;
173 
174 
175   /* construct the graph if necessary */
176   if (!G) {
177     SETERRQ(fcomm,PETSC_ERR_ARG_WRONGSTATE,"Must set Graph in PC in PCGAMG before coarsening");
178   }
179 
180   ierr = MatCoarsenCreate(fcomm,&crs);CHKERRQ(ierr);
181   ierr = MatCoarsenSetFromOptions(crs);CHKERRQ(ierr);
182   ierr = MatCoarsenSetAdjacency(crs,*G);CHKERRQ(ierr);
183   ierr = MatCoarsenSetStrictAggs(crs,PETSC_TRUE);CHKERRQ(ierr);
184   ierr = MatCoarsenApply(crs);CHKERRQ(ierr);
185   ierr = MatCoarsenGetData(crs,agg_lists);CHKERRQ(ierr);
186   ierr = MatCoarsenDestroy(&crs);CHKERRQ(ierr);
187 
188   PetscFunctionReturn(0);
189 }
190 
191 #undef __FUNCT__
192 #define __FUNCT__ "PCGAMGClassicalGhost_Private"
193 /*
194  Find all ghost nodes that are coarse and output the fine/coarse splitting for those as well
195 
196  Input:
197  G - graph;
198  gvec - Global Vector
199  avec - Local part of the scattered vec
200  bvec - Global part of the scattered vec
201 
202  Output:
203  findx - indirection t
204 
205  */
206 PetscErrorCode PCGAMGClassicalGhost_Private(Mat G,Vec v,Vec gv)
207 {
208   PetscErrorCode ierr;
209   Mat_MPIAIJ     *aij = (Mat_MPIAIJ*)G->data;
210   PetscBool      isMPIAIJ;
211 
212   PetscFunctionBegin;
213   ierr = PetscObjectTypeCompare((PetscObject)G, MATMPIAIJ, &isMPIAIJ ); CHKERRQ(ierr);
214   if (isMPIAIJ) {
215     ierr = VecScatterBegin(aij->Mvctx,v,gv,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
216     ierr = VecScatterEnd(aij->Mvctx,v,gv,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
217   }
218   PetscFunctionReturn(0);
219 }
220 
221 #undef __FUNCT__
222 #define __FUNCT__ "PCGAMGProlongator_Classical"
223 PetscErrorCode PCGAMGProlongator_Classical(PC pc, Mat A, Mat G, PetscCoarsenData *agg_lists,Mat *P)
224 {
225   PetscErrorCode    ierr;
226   MPI_Comm          comm;
227   Mat               lG,gG,lA,gA;     /* on and off diagonal matrices */
228   PetscInt          fn;                        /* fine local blocked sizes */
229   PetscInt          cn;                        /* coarse local blocked sizes */
230   PetscInt          gn;                        /* size of the off-diagonal fine vector */
231   PetscInt          fs,fe;                     /* fine (row) ownership range*/
232   PetscInt          cs,ce;                     /* coarse (column) ownership range */
233   PetscInt          i,j,k;                     /* indices! */
234   PetscBool         iscoarse;                  /* flag for determining if a node is coarse */
235   PetscInt          *lcid,*gcid;               /* on and off-processor coarse unknown IDs */
236   PetscInt          *lsparse,*gsparse;         /* on and off-processor sparsity patterns for prolongator */
237   PetscScalar       pij;
238   const PetscScalar *rval;
239   const PetscInt    *rcol;
240   PetscScalar       g_pos,g_neg,a_pos,a_neg,diag,invdiag,alpha,beta;
241   Vec               F;   /* vec of coarse size */
242   Vec               C;   /* vec of fine size */
243   Vec               gF;  /* vec of off-diagonal fine size */
244   MatType           mtype;
245   PetscInt          c_indx;
246   const PetscScalar *vcols;
247   const PetscInt    *icols;
248   PetscScalar       c_scalar;
249   PetscInt          ncols,col;
250   PetscInt          row_f,row_c;
251   PetscInt          cmax=0,ncolstotal,idx;
252   PetscScalar       *pvals;
253   PetscInt          *pcols;
254 
255   PetscFunctionBegin;
256   comm = ((PetscObject)pc)->comm;
257   ierr = MatGetOwnershipRange(A,&fs,&fe); CHKERRQ(ierr);
258   fn = (fe - fs);
259 
260   ierr = MatGetVecs(A,&F,NULL);CHKERRQ(ierr);
261 
262   /* get the number of local unknowns and the indices of the local unknowns */
263 
264   ierr = PetscMalloc(sizeof(PetscInt)*fn,&lsparse);CHKERRQ(ierr);
265   ierr = PetscMalloc(sizeof(PetscInt)*fn,&gsparse);CHKERRQ(ierr);
266   ierr = PetscMalloc(sizeof(PetscInt)*fn,&lcid);CHKERRQ(ierr);
267 
268   /* count the number of coarse unknowns */
269   cn = 0;
270   for (i=0;i<fn;i++) {
271     /* filter out singletons */
272     ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr);
273     lcid[i] = -1;
274     if (!iscoarse) {
275       cn++;
276     }
277   }
278 
279    /* create the coarse vector */
280   ierr = VecCreateMPI(comm,cn,PETSC_DECIDE,&C);CHKERRQ(ierr);
281   ierr = VecGetOwnershipRange(C,&cs,&ce);CHKERRQ(ierr);
282 
283   /* construct a global vector indicating the global indices of the coarse unknowns */
284   cn = 0;
285   for (i=0;i<fn;i++) {
286     ierr = PetscCDEmptyAt(agg_lists,i,&iscoarse); CHKERRQ(ierr);
287     if (!iscoarse) {
288       lcid[i] = cs+cn;
289       cn++;
290     } else {
291       lcid[i] = -1;
292     }
293     c_scalar = (PetscScalar)lcid[i];
294     c_indx = fs+i;
295     ierr = VecSetValues(F,1,&c_indx,&c_scalar,INSERT_VALUES);CHKERRQ(ierr);
296   }
297 
298   ierr = VecAssemblyBegin(F);CHKERRQ(ierr);
299   ierr = VecAssemblyEnd(F);CHKERRQ(ierr);
300 
301   /* split the graph into two */
302   ierr = PCGAMGClassicalGraphSplitting_Private(G,&lG,&gG);CHKERRQ(ierr);
303   ierr = PCGAMGClassicalGraphSplitting_Private(A,&lA,&gA);CHKERRQ(ierr);
304 
305   /* scatter to the ghost vector */
306   ierr = PCGAMGClassicalCreateGhostVector_Private(G,&gF,NULL);CHKERRQ(ierr);
307   ierr = PCGAMGClassicalGhost_Private(G,F,gF);CHKERRQ(ierr);
308 
309   if (gG) {
310     ierr = VecGetSize(gF,&gn);CHKERRQ(ierr);
311     ierr = PetscMalloc(sizeof(PetscInt)*gn,&gcid);CHKERRQ(ierr);
312     for (i=0;i<gn;i++) {
313       ierr = VecGetValues(gF,1,&i,&c_scalar);CHKERRQ(ierr);
314       gcid[i] = (PetscInt)PetscRealPart(c_scalar);
315     }
316   }
317 
318   ierr = VecDestroy(&F);CHKERRQ(ierr);
319   ierr = VecDestroy(&gF);CHKERRQ(ierr);
320   ierr = VecDestroy(&C);CHKERRQ(ierr);
321 
322   /* count the on and off processor sparsity patterns for the prolongator */
323 
324   for (i=0;i<fn;i++) {
325     /* on */
326     ierr = MatGetRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
327     ncolstotal = ncols;
328     lsparse[i] = 0;
329     if (lcid[i] >= 0) {
330       lsparse[i] = 1;
331       gsparse[i] = 0;
332     } else {
333       for (j = 0;j < ncols;j++) {
334         col = icols[j];
335         if (lcid[col] >= 0 && vcols[j] != 0.) {
336           lsparse[i] += 1;
337         }
338       }
339       ierr = MatRestoreRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
340       ncolstotal += ncols;
341       /* off */
342       gsparse[i] = 0;
343       if (gG) {
344         ierr = MatGetRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
345         for (j = 0; j < ncols; j++) {
346           col = icols[j];
347           if (gcid[col] >= 0 && vcols[j] != 0.) {
348             gsparse[i] += 1;
349           }
350         }
351         ierr = MatRestoreRow(gG,i,&ncols,NULL,NULL);CHKERRQ(ierr);
352       }
353       if (ncolstotal > cmax) cmax = ncolstotal;
354     }
355   }
356 
357   ierr = PetscMalloc(sizeof(PetscInt)*cmax,&pcols);CHKERRQ(ierr);
358   ierr = PetscMalloc(sizeof(PetscScalar)*cmax,&pvals);CHKERRQ(ierr);
359 
360   /* preallocate and create the prolongator */
361   ierr = MatCreate(comm,P); CHKERRQ(ierr);
362   ierr = MatGetType(G,&mtype);CHKERRQ(ierr);
363   ierr = MatSetType(*P,mtype);CHKERRQ(ierr);
364 
365   ierr = MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
366   ierr = MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);CHKERRQ(ierr);
367   ierr = MatSeqAIJSetPreallocation(*P,0,lsparse);CHKERRQ(ierr);
368 
369   /* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
370   for (i = 0;i < fn;i++) {
371     /* determine on or off */
372     row_f = i + fs;
373     row_c = lcid[i];
374     if (row_c >= 0) {
375       pij = 1.;
376       ierr = MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);CHKERRQ(ierr);
377     } else {
378       g_pos = 0.;
379       g_neg = 0.;
380       a_pos = 0.;
381       a_neg = 0.;
382       diag = 0.;
383 
384       /* local strong connections */
385       ierr = MatGetRow(lG,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
386       for (k = 0; k < ncols; k++) {
387         if (PetscRealPart(rval[k]) > 0) {
388           g_pos += rval[k];
389         } else {
390           g_neg += rval[k];
391         }
392       }
393       ierr = MatRestoreRow(lG,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
394 
395       /* ghosted strong connections */
396       if (gG) {
397         ierr = MatGetRow(gG,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
398         for (k = 0; k < ncols; k++) {
399           if (PetscRealPart(rval[k]) > 0.) {
400             g_pos += rval[k];
401           } else {
402             g_neg += rval[k];
403           }
404         }
405         ierr = MatRestoreRow(gG,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
406       }
407 
408       /* local all connections */
409       ierr = MatGetRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
410       for (k = 0; k < ncols; k++) {
411         if (PetscRealPart(rval[k]) > 0) {
412           a_pos += rval[k];
413         } else {
414           a_neg += rval[k];
415         }
416         if (rcol[k] == i) {
417           diag = rval[k];
418         }
419       }
420       ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
421 
422       /* ghosted all connections */
423       if (gA) {
424         ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
425         for (k = 0; k < ncols; k++) {
426           if (PetscRealPart(rval[k]) > 0.) {
427             a_pos += PetscRealPart(rval[k]);
428           } else {
429             a_neg += PetscRealPart(rval[k]);
430           }
431         }
432         ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
433       }
434 
435       if (g_neg == 0.) {
436         alpha = 0.;
437       } else {
438         alpha = -a_neg/g_neg;
439       }
440 
441       if (g_pos == 0.) {
442         diag += a_pos;
443         beta = 0.;
444       } else {
445         beta = -a_pos/g_pos;
446       }
447 
448       invdiag = 1. / diag;
449       /* on */
450       ierr = MatGetRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
451       idx = 0;
452       for (j = 0;j < ncols;j++) {
453         col = icols[j];
454         if (lcid[col] >= 0 && vcols[j] != 0.) {
455           row_f = i + fs;
456           row_c = lcid[col];
457           /* set the values for on-processor ones */
458           if (PetscRealPart(vcols[j]) < 0.) {
459             pij = vcols[j]*alpha*invdiag;
460           } else {
461             pij = vcols[j]*beta*invdiag;
462           }
463           if (PetscAbsScalar(pij) != 0.) {
464             pvals[idx] = pij;
465             pcols[idx] = row_c;
466             idx++;
467           }
468         }
469       }
470       ierr = MatRestoreRow(lG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
471       /* off */
472       if (gG) {
473         ierr = MatGetRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
474         for (j = 0; j < ncols; j++) {
475           col = icols[j];
476           if (gcid[col] >= 0 && vcols[j] != 0.) {
477             row_f = i + fs;
478             row_c = gcid[col];
479             /* set the values for on-processor ones */
480             if (PetscRealPart(vcols[j]) < 0.) {
481               pij = vcols[j]*alpha*invdiag;
482             } else {
483               pij = vcols[j]*beta*invdiag;
484             }
485             if (PetscAbsScalar(pij) != 0.) {
486               pvals[idx] = pij;
487               pcols[idx] = row_c;
488               idx++;
489             }
490           }
491         }
492         ierr = MatRestoreRow(gG,i,&ncols,&icols,&vcols);CHKERRQ(ierr);
493       }
494       ierr = MatSetValues(*P,1,&row_f,idx,pcols,pvals,INSERT_VALUES);CHKERRQ(ierr);
495     }
496   }
497 
498   ierr = MatAssemblyBegin(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
499   ierr = MatAssemblyEnd(*P, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
500 
501   ierr = PetscFree(lsparse);CHKERRQ(ierr);
502   ierr = PetscFree(gsparse);CHKERRQ(ierr);
503   ierr = PetscFree(pcols);CHKERRQ(ierr);
504   ierr = PetscFree(pvals);CHKERRQ(ierr);
505   ierr = PetscFree(lcid);CHKERRQ(ierr);
506   if (gG) {ierr = PetscFree(gcid);CHKERRQ(ierr);}
507 
508   PetscFunctionReturn(0);
509 }
510 
511 #undef __FUNCT__
512 #define __FUNCT__ "PCGAMGDestroy_Classical"
513 PetscErrorCode PCGAMGDestroy_Classical(PC pc)
514 {
515   PetscErrorCode ierr;
516   PC_MG          *mg          = (PC_MG*)pc->data;
517   PC_GAMG        *pc_gamg     = (PC_GAMG*)mg->innerctx;
518 
519   PetscFunctionBegin;
520   ierr = PetscFree(pc_gamg->subctx);CHKERRQ(ierr);
521   PetscFunctionReturn(0);
522 }
523 
524 #undef __FUNCT__
525 #define __FUNCT__ "PCGAMGSetFromOptions_Classical"
526 PetscErrorCode PCGAMGSetFromOptions_Classical(PC pc)
527 {
528   PetscFunctionBegin;
529   PetscFunctionReturn(0);
530 }
531 
532 #undef __FUNCT__
533 #define __FUNCT__ "PCGAMGSetData_Classical"
534 PetscErrorCode PCGAMGSetData_Classical(PC pc, Mat A)
535 {
536   PC_MG          *mg      = (PC_MG*)pc->data;
537   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
538 
539   PetscFunctionBegin;
540   /* no data for classical AMG */
541   pc_gamg->data           = NULL;
542   pc_gamg->data_cell_cols = 1;
543   pc_gamg->data_cell_rows = 1;
544   pc_gamg->data_sz = 0;
545   PetscFunctionReturn(0);
546 }
547 
548 /* -------------------------------------------------------------------------- */
549 /*
550    PCCreateGAMG_Classical
551 
552 */
553 #undef __FUNCT__
554 #define __FUNCT__ "PCCreateGAMG_Classical"
555 PetscErrorCode  PCCreateGAMG_Classical(PC pc)
556 {
557   PetscErrorCode ierr;
558   PC_MG             *mg      = (PC_MG*)pc->data;
559   PC_GAMG           *pc_gamg = (PC_GAMG*)mg->innerctx;
560   PC_GAMG_Classical *pc_gamg_classical;
561 
562   PetscFunctionBegin;
563   if (pc_gamg->subctx) {
564     /* call base class */
565     ierr = PCDestroy_GAMG(pc);CHKERRQ(ierr);
566   }
567 
568   /* create sub context for SA */
569   ierr = PetscNewLog(pc, PC_GAMG_Classical, &pc_gamg_classical);CHKERRQ(ierr);
570   pc_gamg->subctx = pc_gamg_classical;
571   pc->ops->setfromoptions = PCGAMGSetFromOptions_Classical;
572   /* reset does not do anything; setup not virtual */
573 
574   /* set internal function pointers */
575   pc_gamg->ops->destroy     = PCGAMGDestroy_Classical;
576   pc_gamg->ops->graph       = PCGAMGGraph_Classical;
577   pc_gamg->ops->coarsen     = PCGAMGCoarsen_Classical;
578   pc_gamg->ops->prolongator = PCGAMGProlongator_Classical;
579   pc_gamg->ops->optprol     = NULL;
580   pc_gamg->ops->formkktprol = NULL;
581 
582   pc_gamg->ops->createdefaultdata = PCGAMGSetData_Classical;
583   PetscFunctionReturn(0);
584 }
585