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