Actual source code: gamg.c
1: /*
2: GAMG geometric-algebric multigrid PC - Mark Adams 2011
3: */
4: #include <petsc/private/matimpl.h>
5: #include <../src/ksp/pc/impls/gamg/gamg.h>
6: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>
8: #if defined(PETSC_HAVE_CUDA)
9: #include <cuda_runtime.h>
10: #endif
12: #if defined(PETSC_HAVE_HIP)
13: #include <hip/hip_runtime.h>
14: #endif
16: PetscLogEvent petsc_gamg_setup_events[GAMG_NUM_SET];
17: PetscLogEvent petsc_gamg_setup_matmat_events[PETSC_MG_MAXLEVELS][3];
19: // #define GAMG_STAGES
20: #if defined(GAMG_STAGES)
21: static PetscLogStage gamg_stages[PETSC_MG_MAXLEVELS];
22: #endif
24: static PetscFunctionList GAMGList = NULL;
25: static PetscBool PCGAMGPackageInitialized;
27: PetscErrorCode PCReset_GAMG(PC pc)
28: {
29: PC_MG *mg = (PC_MG *)pc->data;
30: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
32: PetscFunctionBegin;
33: PetscCall(PetscFree(pc_gamg->data));
34: pc_gamg->data_sz = 0;
35: PetscCall(PetscFree(pc_gamg->orig_data));
36: for (PetscInt level = 0; level < PETSC_MG_MAXLEVELS; level++) {
37: mg->min_eigen_DinvA[level] = 0;
38: mg->max_eigen_DinvA[level] = 0;
39: }
40: pc_gamg->emin = 0;
41: pc_gamg->emax = 0;
42: PetscCall(PCReset_MG(pc));
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*
47: PCGAMGCreateLevel_GAMG: create coarse op with RAP. repartition and/or reduce number
48: of active processors.
50: Input Parameter:
51: . pc - parameters + side effect: coarse data in 'pc_gamg->data' and
52: 'pc_gamg->data_sz' are changed via repartitioning/reduction.
53: . Amat_fine - matrix on this fine (k) level
54: . cr_bs - coarse block size
55: In/Output Parameter:
56: . a_P_inout - prolongation operator to the next level (k-->k-1)
57: . a_nactive_proc - number of active procs
58: Output Parameter:
59: . a_Amat_crs - coarse matrix that is created (k-1)
60: */
61: static PetscErrorCode PCGAMGCreateLevel_GAMG(PC pc, Mat Amat_fine, PetscInt cr_bs, Mat *a_P_inout, Mat *a_Amat_crs, PetscMPIInt *a_nactive_proc, IS *Pcolumnperm, PetscBool is_last)
62: {
63: PC_MG *mg = (PC_MG *)pc->data;
64: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
65: Mat Cmat, Pold = *a_P_inout;
66: MPI_Comm comm;
67: PetscMPIInt rank, size, new_size, nactive = *a_nactive_proc;
68: PetscInt ncrs_eq, ncrs, f_bs;
70: PetscFunctionBegin;
71: PetscCall(PetscObjectGetComm((PetscObject)Amat_fine, &comm));
72: PetscCallMPI(MPI_Comm_rank(comm, &rank));
73: PetscCallMPI(MPI_Comm_size(comm, &size));
74: PetscCall(MatGetBlockSize(Amat_fine, &f_bs));
75: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
76: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
77: PetscCall(MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat));
78: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
79: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
81: if (Pcolumnperm) *Pcolumnperm = NULL;
83: /* set 'ncrs' (nodes), 'ncrs_eq' (equations)*/
84: PetscCall(MatGetLocalSize(Cmat, &ncrs_eq, NULL));
85: if (pc_gamg->data_cell_rows > 0) {
86: ncrs = pc_gamg->data_sz / pc_gamg->data_cell_cols / pc_gamg->data_cell_rows;
87: } else {
88: PetscInt bs;
89: PetscCall(MatGetBlockSize(Cmat, &bs));
90: ncrs = ncrs_eq / bs;
91: }
92: /* get number of PEs to make active 'new_size', reduce, can be any integer 1-P */
93: if (pc_gamg->level_reduction_factors[pc_gamg->current_level] == 0 && PetscDefined(HAVE_CUDA) && pc_gamg->current_level == 0) { /* 0 turns reducing to 1 process/device on; do for HIP, etc. */
94: #if defined(PETSC_HAVE_CUDA)
95: PetscShmComm pshmcomm;
96: PetscMPIInt locrank;
97: MPI_Comm loccomm;
98: PetscInt s_nnodes, r_nnodes, new_new_size;
99: cudaError_t cerr;
100: int devCount;
101: PetscCall(PetscShmCommGet(comm, &pshmcomm));
102: PetscCall(PetscShmCommGetMpiShmComm(pshmcomm, &loccomm));
103: PetscCallMPI(MPI_Comm_rank(loccomm, &locrank));
104: s_nnodes = !locrank;
105: PetscCall(MPIU_Allreduce(&s_nnodes, &r_nnodes, 1, MPIU_INT, MPI_SUM, comm));
106: PetscCheck((size % r_nnodes) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "odd number of nodes np=%d nnodes%" PetscInt_FMT, size, r_nnodes);
107: devCount = 0;
108: cerr = cudaGetDeviceCount(&devCount);
109: cudaGetLastError(); /* Reset the last error */
110: if (cerr == cudaSuccess && devCount >= 1) { /* There are devices, else go to heuristic */
111: new_new_size = r_nnodes * devCount;
112: new_size = new_new_size;
113: PetscCall(PetscInfo(pc, "%s: Fine grid with Cuda. %" PetscInt_FMT " nodes. Change new active set size %d --> %d (devCount=%d #nodes=%" PetscInt_FMT ")\n", ((PetscObject)pc)->prefix, r_nnodes, nactive, new_size, devCount, r_nnodes));
114: } else {
115: PetscCall(PetscInfo(pc, "%s: With Cuda but no device. Use heuristics.", ((PetscObject)pc)->prefix));
116: goto HEURISTIC;
117: }
118: #else
119: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "should not be here");
120: #endif
121: } else if (pc_gamg->level_reduction_factors[pc_gamg->current_level] > 0) {
122: PetscCheck(nactive % pc_gamg->level_reduction_factors[pc_gamg->current_level] == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "odd number of active process %d wrt reduction factor %" PetscInt_FMT, nactive, pc_gamg->level_reduction_factors[pc_gamg->current_level]);
123: new_size = nactive / pc_gamg->level_reduction_factors[pc_gamg->current_level];
124: PetscCall(PetscInfo(pc, "%s: Manually setting reduction to %d active processes (%d/%" PetscInt_FMT ")\n", ((PetscObject)pc)->prefix, new_size, nactive, pc_gamg->level_reduction_factors[pc_gamg->current_level]));
125: } else if (is_last && !pc_gamg->use_parallel_coarse_grid_solver) {
126: new_size = 1;
127: PetscCall(PetscInfo(pc, "%s: Force coarsest grid reduction to %d active processes\n", ((PetscObject)pc)->prefix, new_size));
128: } else {
129: PetscInt ncrs_eq_glob;
130: #if defined(PETSC_HAVE_CUDA)
131: HEURISTIC:
132: #endif
133: PetscCall(MatGetSize(Cmat, &ncrs_eq_glob, NULL));
134: new_size = (PetscMPIInt)((float)ncrs_eq_glob / (float)pc_gamg->min_eq_proc + 0.5); /* hardwire min. number of eq/proc */
135: if (!new_size) new_size = 1; /* not likely, possible? */
136: else if (new_size >= nactive) new_size = nactive; /* no change, rare */
137: PetscCall(PetscInfo(pc, "%s: Coarse grid reduction from %d to %d active processes\n", ((PetscObject)pc)->prefix, nactive, new_size));
138: }
139: if (new_size == nactive) {
140: *a_Amat_crs = Cmat; /* output - no repartitioning or reduction - could bail here */
141: if (new_size < size) {
142: /* odd case where multiple coarse grids are on one processor or no coarsening ... */
143: PetscCall(PetscInfo(pc, "%s: reduced grid using same number of processors (%d) as last grid (use larger coarse grid)\n", ((PetscObject)pc)->prefix, nactive));
144: if (pc_gamg->cpu_pin_coarse_grids) {
145: PetscCall(MatBindToCPU(*a_Amat_crs, PETSC_TRUE));
146: PetscCall(MatBindToCPU(*a_P_inout, PETSC_TRUE));
147: }
148: }
149: /* we know that the grid structure can be reused in MatPtAP */
150: } else { /* reduce active processors - we know that the grid structure can NOT be reused in MatPtAP */
151: PetscInt *counts, *newproc_idx, ii, jj, kk, strideNew, *tidx, ncrs_new, ncrs_eq_new, nloc_old, expand_factor = 1, rfactor = 1;
152: IS is_eq_newproc, is_eq_num, is_eq_num_prim, new_eq_indices;
153: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
154: nloc_old = ncrs_eq / cr_bs;
155: PetscCheck(ncrs_eq % cr_bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ncrs_eq %" PetscInt_FMT " not divisible by cr_bs %" PetscInt_FMT, ncrs_eq, cr_bs);
156: /* get new_size and rfactor */
157: if (pc_gamg->layout_type == PCGAMG_LAYOUT_SPREAD || !pc_gamg->repart) {
158: /* find factor */
159: if (new_size == 1) rfactor = size; /* don't modify */
160: else {
161: PetscReal best_fact = 0.;
162: jj = -1;
163: for (kk = 1; kk <= size; kk++) {
164: if (!(size % kk)) { /* a candidate */
165: PetscReal nactpe = (PetscReal)size / (PetscReal)kk, fact = nactpe / (PetscReal)new_size;
166: if (fact > 1.0) fact = 1. / fact; /* keep fact < 1 */
167: if (fact > best_fact) {
168: best_fact = fact;
169: jj = kk;
170: }
171: }
172: }
173: if (jj != -1) rfactor = jj;
174: else rfactor = 1; /* a prime */
175: if (pc_gamg->layout_type == PCGAMG_LAYOUT_COMPACT) expand_factor = 1;
176: else expand_factor = rfactor;
177: }
178: new_size = size / rfactor; /* make new size one that is factor */
179: if (new_size == nactive) { /* no repartitioning or reduction, bail out because nested here (rare) */
180: *a_Amat_crs = Cmat;
181: PetscCall(PetscInfo(pc, "%s: Finding factorable processor set stopped reduction: new_size=%d, neq(loc)=%" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, new_size, ncrs_eq));
182: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
185: }
186: /* make 'is_eq_newproc' */
187: PetscCall(PetscMalloc1(size, &counts));
188: if (pc_gamg->repart) { /* Repartition Cmat_{k} and move columns of P^{k}_{k-1} and coordinates of primal part accordingly */
189: Mat adj;
190: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_REPART], 0, 0, 0, 0));
191: PetscCall(PetscInfo(pc, "%s: Repartition: size (active): %d --> %d, %" PetscInt_FMT " local equations, using %s process layout\n", ((PetscObject)pc)->prefix, *a_nactive_proc, new_size, ncrs_eq, (pc_gamg->layout_type == PCGAMG_LAYOUT_COMPACT) ? "compact" : "spread"));
192: /* get 'adj' */
193: if (cr_bs == 1) {
194: PetscCall(MatConvert(Cmat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj));
195: } else {
196: /* make a scalar matrix to partition (no Stokes here) */
197: Mat tMat;
198: PetscInt Istart_crs, Iend_crs, ncols, jj, Ii;
199: const PetscScalar *vals;
200: const PetscInt *idx;
201: PetscInt *d_nnz, *o_nnz, M, N;
202: static PetscInt llev = 0; /* ugly but just used for debugging */
203: MatType mtype;
205: PetscCall(PetscMalloc2(ncrs, &d_nnz, ncrs, &o_nnz));
206: PetscCall(MatGetOwnershipRange(Cmat, &Istart_crs, &Iend_crs));
207: PetscCall(MatGetSize(Cmat, &M, &N));
208: for (Ii = Istart_crs, jj = 0; Ii < Iend_crs; Ii += cr_bs, jj++) {
209: PetscCall(MatGetRow(Cmat, Ii, &ncols, NULL, NULL));
210: d_nnz[jj] = ncols / cr_bs;
211: o_nnz[jj] = ncols / cr_bs;
212: PetscCall(MatRestoreRow(Cmat, Ii, &ncols, NULL, NULL));
213: if (d_nnz[jj] > ncrs) d_nnz[jj] = ncrs;
214: if (o_nnz[jj] > (M / cr_bs - ncrs)) o_nnz[jj] = M / cr_bs - ncrs;
215: }
217: PetscCall(MatGetType(Amat_fine, &mtype));
218: PetscCall(MatCreate(comm, &tMat));
219: PetscCall(MatSetSizes(tMat, ncrs, ncrs, PETSC_DETERMINE, PETSC_DETERMINE));
220: PetscCall(MatSetType(tMat, mtype));
221: PetscCall(MatSeqAIJSetPreallocation(tMat, 0, d_nnz));
222: PetscCall(MatMPIAIJSetPreallocation(tMat, 0, d_nnz, 0, o_nnz));
223: PetscCall(PetscFree2(d_nnz, o_nnz));
225: for (ii = Istart_crs; ii < Iend_crs; ii++) {
226: PetscInt dest_row = ii / cr_bs;
227: PetscCall(MatGetRow(Cmat, ii, &ncols, &idx, &vals));
228: for (jj = 0; jj < ncols; jj++) {
229: PetscInt dest_col = idx[jj] / cr_bs;
230: PetscScalar v = 1.0;
231: PetscCall(MatSetValues(tMat, 1, &dest_row, 1, &dest_col, &v, ADD_VALUES));
232: }
233: PetscCall(MatRestoreRow(Cmat, ii, &ncols, &idx, &vals));
234: }
235: PetscCall(MatAssemblyBegin(tMat, MAT_FINAL_ASSEMBLY));
236: PetscCall(MatAssemblyEnd(tMat, MAT_FINAL_ASSEMBLY));
238: if (llev++ == -1) {
239: PetscViewer viewer;
240: char fname[32];
241: PetscCall(PetscSNPrintf(fname, sizeof(fname), "part_mat_%" PetscInt_FMT ".mat", llev));
242: PetscCall(PetscViewerBinaryOpen(comm, fname, FILE_MODE_WRITE, &viewer));
243: PetscCall(MatView(tMat, viewer));
244: PetscCall(PetscViewerDestroy(&viewer));
245: }
246: PetscCall(MatConvert(tMat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj));
247: PetscCall(MatDestroy(&tMat));
248: } /* create 'adj' */
250: { /* partition: get newproc_idx */
251: char prefix[256];
252: const char *pcpre;
253: const PetscInt *is_idx;
254: MatPartitioning mpart;
255: IS proc_is;
257: PetscCall(MatPartitioningCreate(comm, &mpart));
258: PetscCall(MatPartitioningSetAdjacency(mpart, adj));
259: PetscCall(PCGetOptionsPrefix(pc, &pcpre));
260: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : ""));
261: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix));
262: PetscCall(MatPartitioningSetFromOptions(mpart));
263: PetscCall(MatPartitioningSetNParts(mpart, new_size));
264: PetscCall(MatPartitioningApply(mpart, &proc_is));
265: PetscCall(MatPartitioningDestroy(&mpart));
267: /* collect IS info */
268: PetscCall(PetscMalloc1(ncrs_eq, &newproc_idx));
269: PetscCall(ISGetIndices(proc_is, &is_idx));
270: for (kk = jj = 0; kk < nloc_old; kk++) {
271: for (ii = 0; ii < cr_bs; ii++, jj++) { newproc_idx[jj] = is_idx[kk] * expand_factor; /* distribution */ }
272: }
273: PetscCall(ISRestoreIndices(proc_is, &is_idx));
274: PetscCall(ISDestroy(&proc_is));
275: }
276: PetscCall(MatDestroy(&adj));
278: PetscCall(ISCreateGeneral(comm, ncrs_eq, newproc_idx, PETSC_COPY_VALUES, &is_eq_newproc));
279: PetscCall(PetscFree(newproc_idx));
280: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REPART], 0, 0, 0, 0));
281: } else { /* simple aggregation of parts -- 'is_eq_newproc' */
282: PetscInt targetPE;
283: PetscCheck(new_size != nactive, PETSC_COMM_SELF, PETSC_ERR_PLIB, "new_size==nactive. Should not happen");
284: PetscCall(PetscInfo(pc, "%s: Number of equations (loc) %" PetscInt_FMT " with simple aggregation\n", ((PetscObject)pc)->prefix, ncrs_eq));
285: targetPE = (rank / rfactor) * expand_factor;
286: PetscCall(ISCreateStride(comm, ncrs_eq, targetPE, 0, &is_eq_newproc));
287: } /* end simple 'is_eq_newproc' */
289: /*
290: Create an index set from the is_eq_newproc index set to indicate the mapping TO
291: */
292: PetscCall(ISPartitioningToNumbering(is_eq_newproc, &is_eq_num));
293: is_eq_num_prim = is_eq_num;
294: /*
295: Determine how many equations/vertices are assigned to each processor
296: */
297: PetscCall(ISPartitioningCount(is_eq_newproc, size, counts));
298: ncrs_eq_new = counts[rank];
299: PetscCall(ISDestroy(&is_eq_newproc));
300: ncrs_new = ncrs_eq_new / cr_bs;
302: PetscCall(PetscFree(counts));
303: /* data movement scope -- this could be moved to subclasses so that we don't try to cram all auxiliary data into some complex abstracted thing */
304: {
305: Vec src_crd, dest_crd;
306: const PetscInt *idx, ndata_rows = pc_gamg->data_cell_rows, ndata_cols = pc_gamg->data_cell_cols, node_data_sz = ndata_rows * ndata_cols;
307: VecScatter vecscat;
308: PetscScalar *array;
309: IS isscat;
310: /* move data (for primal equations only) */
311: /* Create a vector to contain the newly ordered element information */
312: PetscCall(VecCreate(comm, &dest_crd));
313: PetscCall(VecSetSizes(dest_crd, node_data_sz * ncrs_new, PETSC_DECIDE));
314: PetscCall(VecSetType(dest_crd, VECSTANDARD)); /* this is needed! */
315: /*
316: There are 'ndata_rows*ndata_cols' data items per node, (one can think of the vectors of having
317: a block size of ...). Note, ISs are expanded into equation space by 'cr_bs'.
318: */
319: PetscCall(PetscMalloc1(ncrs * node_data_sz, &tidx));
320: PetscCall(ISGetIndices(is_eq_num_prim, &idx));
321: for (ii = 0, jj = 0; ii < ncrs; ii++) {
322: PetscInt id = idx[ii * cr_bs] / cr_bs; /* get node back */
323: for (kk = 0; kk < node_data_sz; kk++, jj++) tidx[jj] = id * node_data_sz + kk;
324: }
325: PetscCall(ISRestoreIndices(is_eq_num_prim, &idx));
326: PetscCall(ISCreateGeneral(comm, node_data_sz * ncrs, tidx, PETSC_COPY_VALUES, &isscat));
327: PetscCall(PetscFree(tidx));
328: /*
329: Create a vector to contain the original vertex information for each element
330: */
331: PetscCall(VecCreateSeq(PETSC_COMM_SELF, node_data_sz * ncrs, &src_crd));
332: for (jj = 0; jj < ndata_cols; jj++) {
333: const PetscInt stride0 = ncrs * pc_gamg->data_cell_rows;
334: for (ii = 0; ii < ncrs; ii++) {
335: for (kk = 0; kk < ndata_rows; kk++) {
336: PetscInt ix = ii * ndata_rows + kk + jj * stride0, jx = ii * node_data_sz + kk * ndata_cols + jj;
337: PetscScalar tt = (PetscScalar)pc_gamg->data[ix];
338: PetscCall(VecSetValues(src_crd, 1, &jx, &tt, INSERT_VALUES));
339: }
340: }
341: }
342: PetscCall(VecAssemblyBegin(src_crd));
343: PetscCall(VecAssemblyEnd(src_crd));
344: /*
345: Scatter the element vertex information (still in the original vertex ordering)
346: to the correct processor
347: */
348: PetscCall(VecScatterCreate(src_crd, NULL, dest_crd, isscat, &vecscat));
349: PetscCall(ISDestroy(&isscat));
350: PetscCall(VecScatterBegin(vecscat, src_crd, dest_crd, INSERT_VALUES, SCATTER_FORWARD));
351: PetscCall(VecScatterEnd(vecscat, src_crd, dest_crd, INSERT_VALUES, SCATTER_FORWARD));
352: PetscCall(VecScatterDestroy(&vecscat));
353: PetscCall(VecDestroy(&src_crd));
354: /*
355: Put the element vertex data into a new allocation of the gdata->ele
356: */
357: PetscCall(PetscFree(pc_gamg->data));
358: PetscCall(PetscMalloc1(node_data_sz * ncrs_new, &pc_gamg->data));
360: pc_gamg->data_sz = node_data_sz * ncrs_new;
361: strideNew = ncrs_new * ndata_rows;
363: PetscCall(VecGetArray(dest_crd, &array));
364: for (jj = 0; jj < ndata_cols; jj++) {
365: for (ii = 0; ii < ncrs_new; ii++) {
366: for (kk = 0; kk < ndata_rows; kk++) {
367: PetscInt ix = ii * ndata_rows + kk + jj * strideNew, jx = ii * node_data_sz + kk * ndata_cols + jj;
368: pc_gamg->data[ix] = PetscRealPart(array[jx]);
369: }
370: }
371: }
372: PetscCall(VecRestoreArray(dest_crd, &array));
373: PetscCall(VecDestroy(&dest_crd));
374: }
375: /* move A and P (columns) with new layout */
376: /* PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[SET13],0,0,0,0)); */
377: /*
378: Invert for MatCreateSubMatrix
379: */
380: PetscCall(ISInvertPermutation(is_eq_num, ncrs_eq_new, &new_eq_indices));
381: PetscCall(ISSort(new_eq_indices)); /* is this needed? */
382: PetscCall(ISSetBlockSize(new_eq_indices, cr_bs));
383: if (is_eq_num != is_eq_num_prim) { PetscCall(ISDestroy(&is_eq_num_prim)); /* could be same as 'is_eq_num' */ }
384: if (Pcolumnperm) {
385: PetscCall(PetscObjectReference((PetscObject)new_eq_indices));
386: *Pcolumnperm = new_eq_indices;
387: }
388: PetscCall(ISDestroy(&is_eq_num));
389: /* PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[SET13],0,0,0,0)); */
390: /* PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[SET14],0,0,0,0)); */
392: /* 'a_Amat_crs' output */
393: {
394: Mat mat;
395: PetscBool isset, isspd, isher;
396: #if !defined(PETSC_USE_COMPLEX)
397: PetscBool issym;
398: #endif
400: PetscCall(MatCreateSubMatrix(Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat));
401: PetscCall(MatIsSPDKnown(Cmat, &isset, &isspd)); // like MatPropagateSymmetryOptions, but should set MAT_STRUCTURALLY_SYMMETRIC ?
402: if (isset) PetscCall(MatSetOption(mat, MAT_SPD, isspd));
403: else {
404: PetscCall(MatIsHermitianKnown(Cmat, &isset, &isher));
405: if (isset) PetscCall(MatSetOption(mat, MAT_HERMITIAN, isher));
406: else {
407: #if !defined(PETSC_USE_COMPLEX)
408: PetscCall(MatIsSymmetricKnown(Cmat, &isset, &issym));
409: if (isset) PetscCall(MatSetOption(mat, MAT_SYMMETRIC, issym));
410: #endif
411: }
412: }
413: *a_Amat_crs = mat;
414: }
415: PetscCall(MatDestroy(&Cmat));
417: /* PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[SET14],0,0,0,0)); */
418: /* prolongator */
419: {
420: IS findices;
421: PetscInt Istart, Iend;
422: Mat Pnew;
424: PetscCall(MatGetOwnershipRange(Pold, &Istart, &Iend));
425: /* PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[SET15],0,0,0,0)); */
426: PetscCall(ISCreateStride(comm, Iend - Istart, Istart, 1, &findices));
427: PetscCall(ISSetBlockSize(findices, f_bs));
428: PetscCall(MatCreateSubMatrix(Pold, findices, new_eq_indices, MAT_INITIAL_MATRIX, &Pnew));
429: PetscCall(ISDestroy(&findices));
430: PetscCall(MatSetOption(Pnew, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
432: /* PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[SET15],0,0,0,0)); */
433: PetscCall(MatDestroy(a_P_inout));
435: /* output - repartitioned */
436: *a_P_inout = Pnew;
437: }
438: PetscCall(ISDestroy(&new_eq_indices));
440: *a_nactive_proc = new_size; /* output */
442: /* pinning on reduced grids, not a bad heuristic and optimization gets folded into process reduction optimization */
443: if (pc_gamg->cpu_pin_coarse_grids) {
444: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
445: static PetscInt llev = 2;
446: PetscCall(PetscInfo(pc, "%s: Pinning level %" PetscInt_FMT " to the CPU\n", ((PetscObject)pc)->prefix, llev++));
447: #endif
448: PetscCall(MatBindToCPU(*a_Amat_crs, PETSC_TRUE));
449: PetscCall(MatBindToCPU(*a_P_inout, PETSC_TRUE));
450: if (1) { /* HACK: move this to MatBindCPU_MPIAIJXXX; lvec is created, need to pin it, this is done in MatSetUpMultiply_MPIAIJ. Hack */
451: Mat A = *a_Amat_crs, P = *a_P_inout;
452: PetscMPIInt size;
453: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
454: if (size > 1) {
455: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
456: PetscCall(VecBindToCPU(a->lvec, PETSC_TRUE));
457: PetscCall(VecBindToCPU(p->lvec, PETSC_TRUE));
458: }
459: }
460: }
461: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
462: } // processor reduce
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: // used in GEO
467: PetscErrorCode PCGAMGSquareGraph_GAMG(PC a_pc, Mat Gmat1, Mat *Gmat2)
468: {
469: const char *prefix;
470: char addp[32];
471: PC_MG *mg = (PC_MG *)a_pc->data;
472: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
474: PetscFunctionBegin;
475: PetscCall(PCGetOptionsPrefix(a_pc, &prefix));
476: PetscCall(PetscInfo(a_pc, "%s: Square Graph on level %" PetscInt_FMT "\n", ((PetscObject)a_pc)->prefix, pc_gamg->current_level + 1));
477: PetscCall(MatProductCreate(Gmat1, Gmat1, NULL, Gmat2));
478: PetscCall(MatSetOptionsPrefix(*Gmat2, prefix));
479: PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_square_%" PetscInt_FMT "_", pc_gamg->current_level));
480: PetscCall(MatAppendOptionsPrefix(*Gmat2, addp));
481: if ((*Gmat2)->structurally_symmetric == PETSC_BOOL3_TRUE) {
482: PetscCall(MatProductSetType(*Gmat2, MATPRODUCT_AB));
483: } else {
484: PetscCall(MatSetOption(Gmat1, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
485: PetscCall(MatProductSetType(*Gmat2, MATPRODUCT_AtB));
486: }
487: PetscCall(MatProductSetFromOptions(*Gmat2));
488: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0], 0, 0, 0, 0));
489: PetscCall(MatProductSymbolic(*Gmat2));
490: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0], 0, 0, 0, 0));
491: PetscCall(MatProductClear(*Gmat2));
492: /* we only need the sparsity, cheat and tell PETSc the matrix has been assembled */
493: (*Gmat2)->assembled = PETSC_TRUE;
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: /*
498: PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner
499: by setting data structures and options.
501: Input Parameter:
502: . pc - the preconditioner context
504: */
505: PetscErrorCode PCSetUp_GAMG(PC pc)
506: {
507: PC_MG *mg = (PC_MG *)pc->data;
508: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
509: Mat Pmat = pc->pmat;
510: PetscInt fine_level, level, level1, bs, M, N, qq, lidx, nASMBlocksArr[PETSC_MG_MAXLEVELS];
511: MPI_Comm comm;
512: PetscMPIInt rank, size, nactivepe;
513: Mat Aarr[PETSC_MG_MAXLEVELS], Parr[PETSC_MG_MAXLEVELS];
514: IS *ASMLocalIDsArr[PETSC_MG_MAXLEVELS];
515: PetscBool is_last = PETSC_FALSE;
516: #if defined(PETSC_USE_INFO)
517: PetscLogDouble nnz0 = 0., nnztot = 0.;
518: MatInfo info;
519: #endif
521: PetscFunctionBegin;
522: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
523: PetscCallMPI(MPI_Comm_rank(comm, &rank));
524: PetscCallMPI(MPI_Comm_size(comm, &size));
525: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
526: if (pc->setupcalled) {
527: if (!pc_gamg->reuse_prol || pc->flag == DIFFERENT_NONZERO_PATTERN) {
528: /* reset everything */
529: PetscCall(PCReset_MG(pc));
530: pc->setupcalled = 0;
531: } else {
532: PC_MG_Levels **mglevels = mg->levels;
533: /* just do Galerkin grids */
534: Mat B, dA, dB;
535: if (pc_gamg->Nlevels > 1) {
536: PetscInt gl;
537: /* currently only handle case where mat and pmat are the same on coarser levels */
538: PetscCall(KSPGetOperators(mglevels[pc_gamg->Nlevels - 1]->smoothd, &dA, &dB));
539: /* (re)set to get dirty flag */
540: PetscCall(KSPSetOperators(mglevels[pc_gamg->Nlevels - 1]->smoothd, dA, dB));
542: for (level = pc_gamg->Nlevels - 2, gl = 0; level >= 0; level--, gl++) {
543: MatReuse reuse = MAT_INITIAL_MATRIX;
544: #if defined(GAMG_STAGES)
545: PetscCall(PetscLogStagePush(gamg_stages[gl]));
546: #endif
547: /* matrix structure can change from repartitioning or process reduction but don't know if we have process reduction here. Should fix */
548: PetscCall(KSPGetOperators(mglevels[level]->smoothd, NULL, &B));
549: if (B->product) {
550: if (B->product->A == dB && B->product->B == mglevels[level + 1]->interpolate) reuse = MAT_REUSE_MATRIX;
551: }
552: if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDestroy(&mglevels[level]->A));
553: if (reuse == MAT_REUSE_MATRIX) {
554: PetscCall(PetscInfo(pc, "%s: RAP after first solve, reuse matrix level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
555: } else {
556: PetscCall(PetscInfo(pc, "%s: RAP after first solve, new matrix level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
557: }
558: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[gl][1], 0, 0, 0, 0));
559: PetscCall(MatPtAP(dB, mglevels[level + 1]->interpolate, reuse, PETSC_DEFAULT, &B));
560: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[gl][1], 0, 0, 0, 0));
561: if (reuse == MAT_INITIAL_MATRIX) mglevels[level]->A = B;
562: PetscCall(KSPSetOperators(mglevels[level]->smoothd, B, B));
563: dB = B;
564: #if defined(GAMG_STAGES)
565: PetscCall(PetscLogStagePop());
566: #endif
567: }
568: }
570: PetscCall(PCSetUp_MG(pc));
571: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
574: }
576: if (!pc_gamg->data) {
577: if (pc_gamg->orig_data) {
578: PetscCall(MatGetBlockSize(Pmat, &bs));
579: PetscCall(MatGetLocalSize(Pmat, &qq, NULL));
581: pc_gamg->data_sz = (qq / bs) * pc_gamg->orig_data_cell_rows * pc_gamg->orig_data_cell_cols;
582: pc_gamg->data_cell_rows = pc_gamg->orig_data_cell_rows;
583: pc_gamg->data_cell_cols = pc_gamg->orig_data_cell_cols;
585: PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data));
586: for (qq = 0; qq < pc_gamg->data_sz; qq++) pc_gamg->data[qq] = pc_gamg->orig_data[qq];
587: } else {
588: PetscCheck(pc_gamg->ops->createdefaultdata, comm, PETSC_ERR_PLIB, "'createdefaultdata' not set(?) need to support NULL data");
589: PetscCall(pc_gamg->ops->createdefaultdata(pc, Pmat));
590: }
591: }
593: /* cache original data for reuse */
594: if (!pc_gamg->orig_data && (PetscBool)(!pc_gamg->reuse_prol)) {
595: PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->orig_data));
596: for (qq = 0; qq < pc_gamg->data_sz; qq++) pc_gamg->orig_data[qq] = pc_gamg->data[qq];
597: pc_gamg->orig_data_cell_rows = pc_gamg->data_cell_rows;
598: pc_gamg->orig_data_cell_cols = pc_gamg->data_cell_cols;
599: }
601: /* get basic dims */
602: PetscCall(MatGetBlockSize(Pmat, &bs));
603: PetscCall(MatGetSize(Pmat, &M, &N));
605: #if defined(PETSC_USE_INFO)
606: PetscCall(MatGetInfo(Pmat, MAT_GLOBAL_SUM, &info)); /* global reduction */
607: nnz0 = info.nz_used;
608: nnztot = info.nz_used;
609: #endif
610: PetscCall(PetscInfo(pc, "%s: level %d) N=%" PetscInt_FMT ", n data rows=%" PetscInt_FMT ", n data cols=%" PetscInt_FMT ", nnz/row (ave)=%" PetscInt_FMT ", np=%d\n", ((PetscObject)pc)->prefix, 0, M, pc_gamg->data_cell_rows, pc_gamg->data_cell_cols, (PetscInt)(nnz0 / (PetscReal)M + 0.5), size));
612: /* Get A_i and R_i */
613: for (level = 0, Aarr[0] = Pmat, nactivepe = size; level < (pc_gamg->Nlevels - 1) && (!level || M > pc_gamg->coarse_eq_limit); level++) {
614: pc_gamg->current_level = level;
615: PetscCheck(level < PETSC_MG_MAXLEVELS, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Too many levels %" PetscInt_FMT, level);
616: level1 = level + 1;
617: #if defined(GAMG_STAGES)
618: if (!gamg_stages[level]) {
619: char str[32];
620: PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", (int)level));
621: PetscCall(PetscLogStageRegister(str, &gamg_stages[level]));
622: }
623: PetscCall(PetscLogStagePush(gamg_stages[level]));
624: #endif
625: { /* construct prolongator */
626: Mat Gmat;
627: PetscCoarsenData *agg_lists;
628: Mat Prol11;
630: PetscCall(PCGAMGCreateGraph(pc, Aarr[level], &Gmat));
631: PetscCall(pc_gamg->ops->coarsen(pc, &Gmat, &agg_lists));
632: PetscCall(pc_gamg->ops->prolongator(pc, Aarr[level], Gmat, agg_lists, &Prol11));
634: /* could have failed to create new level */
635: if (Prol11) {
636: const char *prefix;
637: char addp[32];
639: /* get new block size of coarse matrices */
640: PetscCall(MatGetBlockSizes(Prol11, NULL, &bs));
642: if (pc_gamg->ops->optprolongator) {
643: /* smooth */
644: PetscCall(pc_gamg->ops->optprolongator(pc, Aarr[level], &Prol11));
645: }
647: if (pc_gamg->use_aggs_in_asm) {
648: PetscInt bs;
649: PetscCall(MatGetBlockSizes(Prol11, &bs, NULL)); // not timed directly, ugly, could remove, but good ASM method
650: PetscCall(PetscCDGetASMBlocks(agg_lists, bs, Gmat, &nASMBlocksArr[level], &ASMLocalIDsArr[level]));
651: }
653: PetscCall(PCGetOptionsPrefix(pc, &prefix));
654: PetscCall(MatSetOptionsPrefix(Prol11, prefix));
655: PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_prolongator_%d_", (int)level));
656: PetscCall(MatAppendOptionsPrefix(Prol11, addp));
657: /* Always generate the transpose with CUDA
658: Such behaviour can be adapted with -pc_gamg_prolongator_ prefixed options */
659: PetscCall(MatSetOption(Prol11, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
660: PetscCall(MatSetFromOptions(Prol11));
661: Parr[level1] = Prol11;
662: } else Parr[level1] = NULL; /* failed to coarsen */
664: PetscCall(MatDestroy(&Gmat));
665: PetscCall(PetscCDDestroy(agg_lists));
666: } /* construct prolongator scope */
667: if (!level) Aarr[0] = Pmat; /* use Pmat for finest level setup */
668: if (!Parr[level1]) { /* failed to coarsen */
669: PetscCall(PetscInfo(pc, "%s: Stop gridding, level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
670: #if defined(GAMG_STAGES)
671: PetscCall(PetscLogStagePop());
672: #endif
673: break;
674: }
675: PetscCall(MatGetSize(Parr[level1], &M, &N)); /* N is next M, a loop test variables */
676: PetscCheck(!is_last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Is last ?");
677: if (N <= pc_gamg->coarse_eq_limit) is_last = PETSC_TRUE;
678: if (level1 == pc_gamg->Nlevels - 1) is_last = PETSC_TRUE;
679: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0));
680: PetscCall(pc_gamg->ops->createlevel(pc, Aarr[level], bs, &Parr[level1], &Aarr[level1], &nactivepe, NULL, is_last));
681: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0));
683: PetscCall(MatGetSize(Aarr[level1], &M, &N)); /* M is loop test variables */
684: #if defined(PETSC_USE_INFO)
685: PetscCall(MatGetInfo(Aarr[level1], MAT_GLOBAL_SUM, &info));
686: nnztot += info.nz_used;
687: #endif
688: PetscCall(PetscInfo(pc, "%s: %d) N=%" PetscInt_FMT ", n data cols=%" PetscInt_FMT ", nnz/row (ave)=%" PetscInt_FMT ", %d active pes\n", ((PetscObject)pc)->prefix, (int)level1, M, pc_gamg->data_cell_cols, (PetscInt)(info.nz_used / (PetscReal)M), nactivepe));
690: #if defined(GAMG_STAGES)
691: PetscCall(PetscLogStagePop());
692: #endif
693: /* stop if one node or one proc -- could pull back for singular problems */
694: if ((pc_gamg->data_cell_cols && M / pc_gamg->data_cell_cols < 2) || (!pc_gamg->data_cell_cols && M / bs < 2)) {
695: PetscCall(PetscInfo(pc, "%s: HARD stop of coarsening on level %" PetscInt_FMT ". Grid too small: %" PetscInt_FMT " block nodes\n", ((PetscObject)pc)->prefix, level, M / bs));
696: level++;
697: break;
698: }
699: } /* levels */
700: PetscCall(PetscFree(pc_gamg->data));
702: PetscCall(PetscInfo(pc, "%s: %" PetscInt_FMT " levels, operator complexity = %g\n", ((PetscObject)pc)->prefix, level + 1, nnztot / nnz0));
703: pc_gamg->Nlevels = level + 1;
704: fine_level = level;
705: PetscCall(PCMGSetLevels(pc, pc_gamg->Nlevels, NULL));
707: if (pc_gamg->Nlevels > 1) { /* don't setup MG if one level */
709: /* set default smoothers & set operators */
710: for (lidx = 1, level = pc_gamg->Nlevels - 2; lidx <= fine_level; lidx++, level--) {
711: KSP smoother;
712: PC subpc;
714: PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
715: PetscCall(KSPGetPC(smoother, &subpc));
717: PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE));
718: /* set ops */
719: PetscCall(KSPSetOperators(smoother, Aarr[level], Aarr[level]));
720: PetscCall(PCMGSetInterpolation(pc, lidx, Parr[level + 1]));
722: /* set defaults */
723: PetscCall(KSPSetType(smoother, KSPCHEBYSHEV));
725: /* set blocks for ASM smoother that uses the 'aggregates' */
726: if (pc_gamg->use_aggs_in_asm) {
727: PetscInt sz;
728: IS *iss;
730: sz = nASMBlocksArr[level];
731: iss = ASMLocalIDsArr[level];
732: PetscCall(PCSetType(subpc, PCASM));
733: PetscCall(PCASMSetOverlap(subpc, 0));
734: PetscCall(PCASMSetType(subpc, PC_ASM_BASIC));
735: if (!sz) {
736: IS is;
737: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_COPY_VALUES, &is));
738: PetscCall(PCASMSetLocalSubdomains(subpc, 1, NULL, &is));
739: PetscCall(ISDestroy(&is));
740: } else {
741: PetscInt kk;
742: PetscCall(PCASMSetLocalSubdomains(subpc, sz, iss, NULL));
743: for (kk = 0; kk < sz; kk++) PetscCall(ISDestroy(&iss[kk]));
744: PetscCall(PetscFree(iss));
745: }
746: ASMLocalIDsArr[level] = NULL;
747: nASMBlocksArr[level] = 0;
748: } else {
749: PetscCall(PCSetType(subpc, PCJACOBI));
750: }
751: }
752: {
753: /* coarse grid */
754: KSP smoother, *k2;
755: PC subpc, pc2;
756: PetscInt ii, first;
757: Mat Lmat = Aarr[(level = pc_gamg->Nlevels - 1)];
758: lidx = 0;
760: PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
761: PetscCall(KSPSetOperators(smoother, Lmat, Lmat));
762: if (!pc_gamg->use_parallel_coarse_grid_solver) {
763: PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE));
764: PetscCall(KSPGetPC(smoother, &subpc));
765: PetscCall(PCSetType(subpc, PCBJACOBI));
766: PetscCall(PCSetUp(subpc));
767: PetscCall(PCBJacobiGetSubKSP(subpc, &ii, &first, &k2));
768: PetscCheck(ii == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " is not one", ii);
769: PetscCall(KSPGetPC(k2[0], &pc2));
770: PetscCall(PCSetType(pc2, PCLU));
771: PetscCall(PCFactorSetShiftType(pc2, MAT_SHIFT_INBLOCKS));
772: PetscCall(KSPSetTolerances(k2[0], PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1));
773: PetscCall(KSPSetType(k2[0], KSPPREONLY));
774: }
775: }
777: /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */
778: PetscObjectOptionsBegin((PetscObject)pc);
779: PetscCall(PCSetFromOptions_MG(pc, PetscOptionsObject));
780: PetscOptionsEnd();
781: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
783: /* set cheby eigen estimates from SA to use in the solver */
784: if (pc_gamg->use_sa_esteig) {
785: for (lidx = 1, level = pc_gamg->Nlevels - 2; level >= 0; lidx++, level--) {
786: KSP smoother;
787: PetscBool ischeb;
789: PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
790: PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb));
791: if (ischeb) {
792: KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data;
794: // The command line will override these settings because KSPSetFromOptions is called in PCSetUp_MG
795: if (mg->max_eigen_DinvA[level] > 0) {
796: // SA uses Jacobi for P; we use SA estimates if the smoother is also Jacobi or if the user explicitly requested it.
797: // TODO: This should test whether it's the same Jacobi variant (DIAG, ROWSUM, etc.)
798: PetscReal emax, emin;
800: emin = mg->min_eigen_DinvA[level];
801: emax = mg->max_eigen_DinvA[level];
802: PetscCall(PetscInfo(pc, "%s: PCSetUp_GAMG: call KSPChebyshevSetEigenvalues on level %" PetscInt_FMT " (N=%" PetscInt_FMT ") with emax = %g emin = %g\n", ((PetscObject)pc)->prefix, level, Aarr[level]->rmap->N, (double)emax, (double)emin));
803: cheb->emin_provided = emin;
804: cheb->emax_provided = emax;
805: }
806: }
807: }
808: }
810: PetscCall(PCSetUp_MG(pc));
812: /* clean up */
813: for (level = 1; level < pc_gamg->Nlevels; level++) {
814: PetscCall(MatDestroy(&Parr[level]));
815: PetscCall(MatDestroy(&Aarr[level]));
816: }
817: } else {
818: KSP smoother;
820: PetscCall(PetscInfo(pc, "%s: One level solver used (system is seen as DD). Using default solver.\n", ((PetscObject)pc)->prefix));
821: PetscCall(PCMGGetSmoother(pc, 0, &smoother));
822: PetscCall(KSPSetOperators(smoother, Aarr[0], Aarr[0]));
823: PetscCall(KSPSetType(smoother, KSPPREONLY));
824: PetscCall(PCSetUp_MG(pc));
825: }
826: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: /*
831: PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner
832: that was created with PCCreate_GAMG().
834: Input Parameter:
835: . pc - the preconditioner context
837: Application Interface Routine: PCDestroy()
838: */
839: PetscErrorCode PCDestroy_GAMG(PC pc)
840: {
841: PC_MG *mg = (PC_MG *)pc->data;
842: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
844: PetscFunctionBegin;
845: PetscCall(PCReset_GAMG(pc));
846: if (pc_gamg->ops->destroy) PetscCall((*pc_gamg->ops->destroy)(pc));
847: PetscCall(PetscFree(pc_gamg->ops));
848: PetscCall(PetscFree(pc_gamg->gamg_type_name));
849: PetscCall(PetscFree(pc_gamg));
850: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", NULL));
851: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", NULL));
852: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", NULL));
853: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", NULL));
854: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", NULL));
855: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", NULL));
856: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", NULL));
857: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseParallelCoarseGridSolve_C", NULL));
858: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", NULL));
859: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", NULL));
860: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", NULL));
861: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", NULL));
862: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", NULL));
863: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", NULL));
864: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", NULL));
865: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", NULL));
866: PetscCall(PCDestroy_MG(pc));
867: PetscFunctionReturn(PETSC_SUCCESS);
868: }
870: /*@
871: PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction in `PCGAMG`
873: Logically Collective
875: Input Parameters:
876: + pc - the preconditioner context
877: - n - the number of equations
879: Options Database Key:
880: . -pc_gamg_process_eq_limit <limit> - set the limit
882: Level: intermediate
884: Note:
885: `PCGAMG` will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process
886: that has degrees of freedom
888: .seealso: `PCGAMG`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()`
889: @*/
890: PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n)
891: {
892: PetscFunctionBegin;
894: PetscTryMethod(pc, "PCGAMGSetProcEqLim_C", (PC, PetscInt), (pc, n));
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n)
899: {
900: PC_MG *mg = (PC_MG *)pc->data;
901: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
903: PetscFunctionBegin;
904: if (n > 0) pc_gamg->min_eq_proc = n;
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: /*@
909: PCGAMGSetCoarseEqLim - Set maximum number of equations on the coarsest grid of `PCGAMG`
911: Collective
913: Input Parameters:
914: + pc - the preconditioner context
915: - n - maximum number of equations to aim for
917: Options Database Key:
918: . -pc_gamg_coarse_eq_limit <limit> - set the limit
920: Level: intermediate
922: Note:
923: For example -pc_gamg_coarse_eq_limit 1000 will stop coarsening once the coarse grid
924: has less than 1000 unknowns.
926: .seealso: `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()`
927: @*/
928: PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n)
929: {
930: PetscFunctionBegin;
932: PetscTryMethod(pc, "PCGAMGSetCoarseEqLim_C", (PC, PetscInt), (pc, n));
933: PetscFunctionReturn(PETSC_SUCCESS);
934: }
936: static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n)
937: {
938: PC_MG *mg = (PC_MG *)pc->data;
939: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
941: PetscFunctionBegin;
942: if (n > 0) pc_gamg->coarse_eq_limit = n;
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: /*@
947: PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids when reducing the number of MPI ranks to use
949: Collective
951: Input Parameters:
952: + pc - the preconditioner context
953: - n - `PETSC_TRUE` or `PETSC_FALSE`
955: Options Database Key:
956: . -pc_gamg_repartition <true,false> - turn on the repartitioning
958: Level: intermediate
960: Note:
961: This will generally improve the loading balancing of the work on each level
963: .seealso: `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`
964: @*/
965: PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n)
966: {
967: PetscFunctionBegin;
969: PetscTryMethod(pc, "PCGAMGSetRepartition_C", (PC, PetscBool), (pc, n));
970: PetscFunctionReturn(PETSC_SUCCESS);
971: }
973: static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n)
974: {
975: PC_MG *mg = (PC_MG *)pc->data;
976: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
978: PetscFunctionBegin;
979: pc_gamg->repart = n;
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: /*@
984: PCGAMGSetUseSAEstEig - Use the eigen estimate from smoothed aggregation for the Chebyshev smoother during the solution process
986: Collective
988: Input Parameters:
989: + pc - the preconditioner context
990: - n - number of its
992: Options Database Key:
993: . -pc_gamg_use_sa_esteig <true,false> - use the eigen estimate
995: Level: advanced
997: Notes:
998: Smoothed aggregation constructs the smoothed prolongator $P = (I - \omega D^{-1} A) T$ where $T$ is the tentative prolongator and $D$ is the diagonal of $A$.
999: Eigenvalue estimates (based on a few `PCCG` or `PCGMRES` iterations) are computed to choose $\omega$ so that this is a stable smoothing operation.
1000: If Chebyshev with Jacobi (diagonal) preconditioning is used for smoothing, then the eigenvalue estimates can be reused during the solution process
1001: This option is only used when the smoother uses Jacobi, and should be turned off if a different `PCJacobiType` is used.
1002: It became default in PETSc 3.17.
1004: .seealso: `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`
1005: @*/
1006: PetscErrorCode PCGAMGSetUseSAEstEig(PC pc, PetscBool n)
1007: {
1008: PetscFunctionBegin;
1010: PetscTryMethod(pc, "PCGAMGSetUseSAEstEig_C", (PC, PetscBool), (pc, n));
1011: PetscFunctionReturn(PETSC_SUCCESS);
1012: }
1014: static PetscErrorCode PCGAMGSetUseSAEstEig_GAMG(PC pc, PetscBool n)
1015: {
1016: PC_MG *mg = (PC_MG *)pc->data;
1017: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1019: PetscFunctionBegin;
1020: pc_gamg->use_sa_esteig = n;
1021: PetscFunctionReturn(PETSC_SUCCESS);
1022: }
1024: /*@
1025: PCGAMGSetEigenvalues - Set WHAT eigenvalues WHY?
1027: Collective
1029: Input Parameters:
1030: + pc - the preconditioner context
1031: - emax - max eigenvalue
1032: . emin - min eigenvalue
1034: Options Database Key:
1035: . -pc_gamg_eigenvalues <emin,emax> - estimates of the eigenvalues
1037: Level: intermediate
1039: .seealso: `PCGAMG`, `PCGAMGSetUseSAEstEig()`
1040: @*/
1041: PetscErrorCode PCGAMGSetEigenvalues(PC pc, PetscReal emax, PetscReal emin)
1042: {
1043: PetscFunctionBegin;
1045: PetscTryMethod(pc, "PCGAMGSetEigenvalues_C", (PC, PetscReal, PetscReal), (pc, emax, emin));
1046: PetscFunctionReturn(PETSC_SUCCESS);
1047: }
1049: static PetscErrorCode PCGAMGSetEigenvalues_GAMG(PC pc, PetscReal emax, PetscReal emin)
1050: {
1051: PC_MG *mg = (PC_MG *)pc->data;
1052: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1054: PetscFunctionBegin;
1055: PetscCheck(emax > emin, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Maximum eigenvalue must be larger than minimum: max %g min %g", (double)emax, (double)emin);
1056: PetscCheck(emax * emin > 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Both eigenvalues must be of the same sign: max %g min %g", (double)emax, (double)emin);
1057: pc_gamg->emax = emax;
1058: pc_gamg->emin = emin;
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: /*@
1063: PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding a `PCGAMG` algebraic multigrid preconditioner
1065: Collective
1067: Input Parameters:
1068: + pc - the preconditioner context
1069: - n - `PETSC_TRUE` or `PETSC_FALSE`
1071: Options Database Key:
1072: . -pc_gamg_reuse_interpolation <true,false> - reuse the previous interpolation
1074: Level: intermediate
1076: Note:
1077: May negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows
1078: rebuilding the preconditioner quicker.
1080: .seealso: `PCGAMG`
1081: @*/
1082: PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n)
1083: {
1084: PetscFunctionBegin;
1086: PetscTryMethod(pc, "PCGAMGSetReuseInterpolation_C", (PC, PetscBool), (pc, n));
1087: PetscFunctionReturn(PETSC_SUCCESS);
1088: }
1090: static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n)
1091: {
1092: PC_MG *mg = (PC_MG *)pc->data;
1093: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1095: PetscFunctionBegin;
1096: pc_gamg->reuse_prol = n;
1097: PetscFunctionReturn(PETSC_SUCCESS);
1098: }
1100: /*@
1101: PCGAMGASMSetUseAggs - Have the `PCGAMG` smoother on each level use the aggregates defined by the coarsening process as the subdomains for the additive Schwarz preconditioner
1102: used as the smoother
1104: Collective
1106: Input Parameters:
1107: + pc - the preconditioner context
1108: - flg - `PETSC_TRUE` to use aggregates, `PETSC_FALSE` to not
1110: Options Database Key:
1111: . -pc_gamg_asm_use_agg <true,false> - use aggregates to define the additive Schwarz subdomains
1113: Level: intermediate
1115: .seealso: `PCGAMG`, `PCASM`, `PCSetType`
1116: @*/
1117: PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg)
1118: {
1119: PetscFunctionBegin;
1121: PetscTryMethod(pc, "PCGAMGASMSetUseAggs_C", (PC, PetscBool), (pc, flg));
1122: PetscFunctionReturn(PETSC_SUCCESS);
1123: }
1125: static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg)
1126: {
1127: PC_MG *mg = (PC_MG *)pc->data;
1128: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1130: PetscFunctionBegin;
1131: pc_gamg->use_aggs_in_asm = flg;
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: /*@
1136: PCGAMGSetUseParallelCoarseGridSolve - allow a parallel coarse grid solver
1138: Collective
1140: Input Parameters:
1141: + pc - the preconditioner context
1142: - flg - `PETSC_TRUE` to not force coarse grid onto one processor
1144: Options Database Key:
1145: . -pc_gamg_use_parallel_coarse_grid_solver - use a parallel coarse grid direct solver
1147: Level: intermediate
1149: .seealso: `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetCpuPinCoarseGrids()`
1150: @*/
1151: PetscErrorCode PCGAMGSetUseParallelCoarseGridSolve(PC pc, PetscBool flg)
1152: {
1153: PetscFunctionBegin;
1155: PetscTryMethod(pc, "PCGAMGSetUseParallelCoarseGridSolve_C", (PC, PetscBool), (pc, flg));
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }
1159: static PetscErrorCode PCGAMGSetUseParallelCoarseGridSolve_GAMG(PC pc, PetscBool flg)
1160: {
1161: PC_MG *mg = (PC_MG *)pc->data;
1162: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1164: PetscFunctionBegin;
1165: pc_gamg->use_parallel_coarse_grid_solver = flg;
1166: PetscFunctionReturn(PETSC_SUCCESS);
1167: }
1169: /*@
1170: PCGAMGSetCpuPinCoarseGrids - pin the coarse grids created in `PCGAMG` to run only on the CPU since the problems may be too small to run efficiently on the GPUs
1172: Collective
1174: Input Parameters:
1175: + pc - the preconditioner context
1176: - flg - `PETSC_TRUE` to pin coarse grids to the CPU
1178: Options Database Key:
1179: . -pc_gamg_cpu_pin_coarse_grids - pin the coarse grids to the CPU
1181: Level: intermediate
1183: .seealso: `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetUseParallelCoarseGridSolve()`
1184: @*/
1185: PetscErrorCode PCGAMGSetCpuPinCoarseGrids(PC pc, PetscBool flg)
1186: {
1187: PetscFunctionBegin;
1189: PetscTryMethod(pc, "PCGAMGSetCpuPinCoarseGrids_C", (PC, PetscBool), (pc, flg));
1190: PetscFunctionReturn(PETSC_SUCCESS);
1191: }
1193: static PetscErrorCode PCGAMGSetCpuPinCoarseGrids_GAMG(PC pc, PetscBool flg)
1194: {
1195: PC_MG *mg = (PC_MG *)pc->data;
1196: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1198: PetscFunctionBegin;
1199: pc_gamg->cpu_pin_coarse_grids = flg;
1200: PetscFunctionReturn(PETSC_SUCCESS);
1201: }
1203: /*@
1204: PCGAMGSetCoarseGridLayoutType - place coarse grids on processors with natural order (compact type)
1206: Collective
1208: Input Parameters:
1209: + pc - the preconditioner context
1210: - flg - `PCGAMGLayoutType` type, either `PCGAMG_LAYOUT_COMPACT` or `PCGAMG_LAYOUT_SPREAD`
1212: Options Database Key:
1213: . -pc_gamg_coarse_grid_layout_type - place the coarse grids with natural ordering
1215: Level: intermediate
1217: .seealso: `PCGAMG`, `PCGAMGSetUseParallelCoarseGridSolve()`, `PCGAMGSetCpuPinCoarseGrids()`, `PCGAMGLayoutType`, `PCGAMG_LAYOUT_COMPACT`, `PCGAMG_LAYOUT_SPREAD`
1218: @*/
1219: PetscErrorCode PCGAMGSetCoarseGridLayoutType(PC pc, PCGAMGLayoutType flg)
1220: {
1221: PetscFunctionBegin;
1223: PetscTryMethod(pc, "PCGAMGSetCoarseGridLayoutType_C", (PC, PCGAMGLayoutType), (pc, flg));
1224: PetscFunctionReturn(PETSC_SUCCESS);
1225: }
1227: static PetscErrorCode PCGAMGSetCoarseGridLayoutType_GAMG(PC pc, PCGAMGLayoutType flg)
1228: {
1229: PC_MG *mg = (PC_MG *)pc->data;
1230: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1232: PetscFunctionBegin;
1233: pc_gamg->layout_type = flg;
1234: PetscFunctionReturn(PETSC_SUCCESS);
1235: }
1237: /*@
1238: PCGAMGSetNlevels - Sets the maximum number of levels `PCGAMG` will use
1240: Not Collective
1242: Input Parameters:
1243: + pc - the preconditioner
1244: - n - the maximum number of levels to use
1246: Options Database Key:
1247: . -pc_mg_levels <n> - set the maximum number of levels to allow
1249: Level: intermediate
1251: Developer Note:
1252: Should be called `PCGAMGSetMaximumNumberlevels()` and possible be shared with `PCMG`
1254: .seealso: `PCGAMG`
1255: @*/
1256: PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n)
1257: {
1258: PetscFunctionBegin;
1260: PetscTryMethod(pc, "PCGAMGSetNlevels_C", (PC, PetscInt), (pc, n));
1261: PetscFunctionReturn(PETSC_SUCCESS);
1262: }
1264: static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n)
1265: {
1266: PC_MG *mg = (PC_MG *)pc->data;
1267: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1269: PetscFunctionBegin;
1270: pc_gamg->Nlevels = n;
1271: PetscFunctionReturn(PETSC_SUCCESS);
1272: }
1274: /*@
1275: PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph
1277: Not Collective
1279: Input Parameters:
1280: + pc - the preconditioner context
1281: . v - array of threshold values for finest n levels; 0.0 means keep all nonzero entries in the graph; negative means keep even zero entries in the graph
1282: - n - number of threshold values provided in array
1284: Options Database Key:
1285: . -pc_gamg_threshold <threshold> - the threshold to drop edges
1287: Level: intermediate
1289: Notes:
1290: Increasing the threshold decreases the rate of coarsening. Conversely reducing the threshold increases the rate of coarsening (aggressive coarsening) and thereby reduces the complexity of the coarse grids, and generally results in slower solver converge rates. Reducing coarse grid complexity reduced the complexity of Galerkin coarse grid construction considerably.
1291: Before coarsening or aggregating the graph, `PCGAMG` removes small values from the graph with this threshold, and thus reducing the coupling in the graph and a different (perhaps better) coarser set of points.
1293: If `n` is less than the total number of coarsenings (see `PCGAMGSetNlevels()`), then threshold scaling (see `PCGAMGSetThresholdScale()`) is used for each successive coarsening.
1294: In this case, `PCGAMGSetThresholdScale()` must be called before `PCGAMGSetThreshold()`.
1295: If `n` is greater than the total number of levels, the excess entries in threshold will not be used.
1297: .seealso: `PCGAMG`, `PCGAMGSetAggressiveLevels()`, `PCGAMGSetThresholdScale()`
1298: @*/
1299: PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal v[], PetscInt n)
1300: {
1301: PetscFunctionBegin;
1304: PetscTryMethod(pc, "PCGAMGSetThreshold_C", (PC, PetscReal[], PetscInt), (pc, v, n));
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal v[], PetscInt n)
1309: {
1310: PC_MG *mg = (PC_MG *)pc->data;
1311: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1312: PetscInt i;
1313: PetscFunctionBegin;
1314: for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->threshold[i] = v[i];
1315: for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale;
1316: PetscFunctionReturn(PETSC_SUCCESS);
1317: }
1319: /*@
1320: PCGAMGSetRankReductionFactors - Set a manual schedule for MPI rank reduction on coarse grids
1322: Collective
1324: Input Parameters:
1325: + pc - the preconditioner context
1326: . v - array of reduction factors. 0 for first value forces a reduction to one process/device on first level in CUDA
1327: - n - number of values provided in array
1329: Options Database Key:
1330: . -pc_gamg_rank_reduction_factors <factors> - provide the schedule
1332: Level: intermediate
1334: .seealso: `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`
1335: @*/
1336: PetscErrorCode PCGAMGSetRankReductionFactors(PC pc, PetscInt v[], PetscInt n)
1337: {
1338: PetscFunctionBegin;
1341: PetscTryMethod(pc, "PCGAMGSetRankReductionFactors_C", (PC, PetscInt[], PetscInt), (pc, v, n));
1342: PetscFunctionReturn(PETSC_SUCCESS);
1343: }
1345: static PetscErrorCode PCGAMGSetRankReductionFactors_GAMG(PC pc, PetscInt v[], PetscInt n)
1346: {
1347: PC_MG *mg = (PC_MG *)pc->data;
1348: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1349: PetscInt i;
1350: PetscFunctionBegin;
1351: for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->level_reduction_factors[i] = v[i];
1352: for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->level_reduction_factors[i] = -1; /* 0 stop putting one process/device on first level */
1353: PetscFunctionReturn(PETSC_SUCCESS);
1354: }
1356: /*@
1357: PCGAMGSetThresholdScale - Relative threshold reduction at each level
1359: Not Collective
1361: Input Parameters:
1362: + pc - the preconditioner context
1363: - scale - the threshold value reduction, usually < 1.0
1365: Options Database Key:
1366: . -pc_gamg_threshold_scale <v> - set the relative threshold reduction on each level
1368: Level: advanced
1370: Note:
1371: The initial threshold (for an arbitrary number of levels starting from the finest) can be set with `PCGAMGSetThreshold()`.
1372: This scaling is used for each subsequent coarsening, but must be called before `PCGAMGSetThreshold()`.
1374: .seealso: `PCGAMG`, `PCGAMGSetThreshold()`
1375: @*/
1376: PetscErrorCode PCGAMGSetThresholdScale(PC pc, PetscReal v)
1377: {
1378: PetscFunctionBegin;
1380: PetscTryMethod(pc, "PCGAMGSetThresholdScale_C", (PC, PetscReal), (pc, v));
1381: PetscFunctionReturn(PETSC_SUCCESS);
1382: }
1384: static PetscErrorCode PCGAMGSetThresholdScale_GAMG(PC pc, PetscReal v)
1385: {
1386: PC_MG *mg = (PC_MG *)pc->data;
1387: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1388: PetscFunctionBegin;
1389: pc_gamg->threshold_scale = v;
1390: PetscFunctionReturn(PETSC_SUCCESS);
1391: }
1393: /*@C
1394: PCGAMGSetType - Set the type of algorithm `PCGAMG` should use
1396: Collective
1398: Input Parameters:
1399: + pc - the preconditioner context
1400: - type - `PCGAMGAGG`, `PCGAMGGEO`, or `PCGAMGCLASSICAL`
1402: Options Database Key:
1403: . -pc_gamg_type <agg,geo,classical> - type of algebraic multigrid to apply
1405: Level: intermediate
1407: .seealso: `PCGAMGGetType()`, `PCGAMG`, `PCGAMGType`
1408: @*/
1409: PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type)
1410: {
1411: PetscFunctionBegin;
1413: PetscTryMethod(pc, "PCGAMGSetType_C", (PC, PCGAMGType), (pc, type));
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: /*@C
1418: PCGAMGGetType - Get the type of algorithm `PCGAMG` will use
1420: Collective
1422: Input Parameter:
1423: . pc - the preconditioner context
1425: Output Parameter:
1426: . type - the type of algorithm used
1428: Level: intermediate
1430: .seealso: `PCGAMG`, `PCGAMGSetType()`, `PCGAMGType`
1431: @*/
1432: PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type)
1433: {
1434: PetscFunctionBegin;
1436: PetscUseMethod(pc, "PCGAMGGetType_C", (PC, PCGAMGType *), (pc, type));
1437: PetscFunctionReturn(PETSC_SUCCESS);
1438: }
1440: static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type)
1441: {
1442: PC_MG *mg = (PC_MG *)pc->data;
1443: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1445: PetscFunctionBegin;
1446: *type = pc_gamg->type;
1447: PetscFunctionReturn(PETSC_SUCCESS);
1448: }
1450: static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type)
1451: {
1452: PC_MG *mg = (PC_MG *)pc->data;
1453: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1454: PetscErrorCode (*r)(PC);
1456: PetscFunctionBegin;
1457: pc_gamg->type = type;
1458: PetscCall(PetscFunctionListFind(GAMGList, type, &r));
1459: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown GAMG type %s given", type);
1460: if (pc_gamg->ops->destroy) {
1461: PetscCall((*pc_gamg->ops->destroy)(pc));
1462: PetscCall(PetscMemzero(pc_gamg->ops, sizeof(struct _PCGAMGOps)));
1463: pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1464: /* cleaning up common data in pc_gamg - this should disappear someday */
1465: pc_gamg->data_cell_cols = 0;
1466: pc_gamg->data_cell_rows = 0;
1467: pc_gamg->orig_data_cell_cols = 0;
1468: pc_gamg->orig_data_cell_rows = 0;
1469: PetscCall(PetscFree(pc_gamg->data));
1470: pc_gamg->data_sz = 0;
1471: }
1472: PetscCall(PetscFree(pc_gamg->gamg_type_name));
1473: PetscCall(PetscStrallocpy(type, &pc_gamg->gamg_type_name));
1474: PetscCall((*r)(pc));
1475: PetscFunctionReturn(PETSC_SUCCESS);
1476: }
1478: static PetscErrorCode PCView_GAMG(PC pc, PetscViewer viewer)
1479: {
1480: PC_MG *mg = (PC_MG *)pc->data;
1481: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1482: PetscReal gc = 0, oc = 0;
1484: PetscFunctionBegin;
1485: PetscCall(PetscViewerASCIIPrintf(viewer, " GAMG specific options\n"));
1486: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for dropping small values in graph on each level ="));
1487: for (PetscInt i = 0; i < mg->nlevels; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)pc_gamg->threshold[i]));
1488: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1489: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold scaling factor for each level not specified = %g\n", (double)pc_gamg->threshold_scale));
1490: if (pc_gamg->use_aggs_in_asm) PetscCall(PetscViewerASCIIPrintf(viewer, " Using aggregates from coarsening process to define subdomains for PCASM\n"));
1491: if (pc_gamg->use_parallel_coarse_grid_solver) PetscCall(PetscViewerASCIIPrintf(viewer, " Using parallel coarse grid solver (all coarse grid equations not put on one process)\n"));
1492: if (pc_gamg->ops->view) PetscCall((*pc_gamg->ops->view)(pc, viewer));
1493: PetscCall(PCMGGetGridComplexity(pc, &gc, &oc));
1494: PetscCall(PetscViewerASCIIPrintf(viewer, " Complexity: grid = %g operator = %g\n", (double)gc, (double)oc));
1495: PetscFunctionReturn(PETSC_SUCCESS);
1496: }
1498: PetscErrorCode PCSetFromOptions_GAMG(PC pc, PetscOptionItems *PetscOptionsObject)
1499: {
1500: PC_MG *mg = (PC_MG *)pc->data;
1501: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1502: PetscBool flag;
1503: MPI_Comm comm;
1504: char prefix[256], tname[32];
1505: PetscInt i, n;
1506: const char *pcpre;
1507: static const char *LayoutTypes[] = {"compact", "spread", "PCGAMGLayoutType", "PC_GAMG_LAYOUT", NULL};
1508: PetscFunctionBegin;
1509: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1510: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG options");
1511: PetscCall(PetscOptionsFList("-pc_gamg_type", "Type of AMG method", "PCGAMGSetType", GAMGList, pc_gamg->gamg_type_name, tname, sizeof(tname), &flag));
1512: if (flag) PetscCall(PCGAMGSetType(pc, tname));
1513: PetscCall(PetscOptionsBool("-pc_gamg_repartition", "Repartion coarse grids", "PCGAMGSetRepartition", pc_gamg->repart, &pc_gamg->repart, NULL));
1514: PetscCall(PetscOptionsBool("-pc_gamg_use_sa_esteig", "Use eigen estimate from smoothed aggregation for smoother", "PCGAMGSetUseSAEstEig", pc_gamg->use_sa_esteig, &pc_gamg->use_sa_esteig, NULL));
1515: PetscCall(PetscOptionsBool("-pc_gamg_reuse_interpolation", "Reuse prolongation operator", "PCGAMGReuseInterpolation", pc_gamg->reuse_prol, &pc_gamg->reuse_prol, NULL));
1516: PetscCall(PetscOptionsBool("-pc_gamg_asm_use_agg", "Use aggregation aggregates for ASM smoother", "PCGAMGASMSetUseAggs", pc_gamg->use_aggs_in_asm, &pc_gamg->use_aggs_in_asm, NULL));
1517: PetscCall(PetscOptionsBool("-pc_gamg_use_parallel_coarse_grid_solver", "Use parallel coarse grid solver (otherwise put last grid on one process)", "PCGAMGSetUseParallelCoarseGridSolve", pc_gamg->use_parallel_coarse_grid_solver, &pc_gamg->use_parallel_coarse_grid_solver, NULL));
1518: PetscCall(PetscOptionsBool("-pc_gamg_cpu_pin_coarse_grids", "Pin coarse grids to the CPU", "PCGAMGSetCpuPinCoarseGrids", pc_gamg->cpu_pin_coarse_grids, &pc_gamg->cpu_pin_coarse_grids, NULL));
1519: PetscCall(PetscOptionsEnum("-pc_gamg_coarse_grid_layout_type", "compact: place reduced grids on processes in natural order; spread: distribute to whole machine for more memory bandwidth", "PCGAMGSetCoarseGridLayoutType", LayoutTypes,
1520: (PetscEnum)pc_gamg->layout_type, (PetscEnum *)&pc_gamg->layout_type, NULL));
1521: PetscCall(PetscOptionsInt("-pc_gamg_process_eq_limit", "Limit (goal) on number of equations per process on coarse grids", "PCGAMGSetProcEqLim", pc_gamg->min_eq_proc, &pc_gamg->min_eq_proc, NULL));
1522: PetscCall(PetscOptionsInt("-pc_gamg_coarse_eq_limit", "Limit on number of equations for the coarse grid", "PCGAMGSetCoarseEqLim", pc_gamg->coarse_eq_limit, &pc_gamg->coarse_eq_limit, NULL));
1523: PetscCall(PetscOptionsReal("-pc_gamg_threshold_scale", "Scaling of threshold for each level not specified", "PCGAMGSetThresholdScale", pc_gamg->threshold_scale, &pc_gamg->threshold_scale, NULL));
1524: n = PETSC_MG_MAXLEVELS;
1525: PetscCall(PetscOptionsRealArray("-pc_gamg_threshold", "Relative threshold to use for dropping edges in aggregation graph", "PCGAMGSetThreshold", pc_gamg->threshold, &n, &flag));
1526: if (!flag || n < PETSC_MG_MAXLEVELS) {
1527: if (!flag) n = 1;
1528: i = n;
1529: do {
1530: pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale;
1531: } while (++i < PETSC_MG_MAXLEVELS);
1532: }
1533: n = PETSC_MG_MAXLEVELS;
1534: PetscCall(PetscOptionsIntArray("-pc_gamg_rank_reduction_factors", "Manual schedule of coarse grid reduction factors that overrides internal heuristics (0 for first reduction puts one process/device)", "PCGAMGSetRankReductionFactors", pc_gamg->level_reduction_factors, &n, &flag));
1535: if (!flag) i = 0;
1536: else i = n;
1537: do {
1538: pc_gamg->level_reduction_factors[i] = -1;
1539: } while (++i < PETSC_MG_MAXLEVELS);
1540: PetscCall(PetscOptionsInt("-pc_mg_levels", "Set number of MG levels", "PCGAMGSetNlevels", pc_gamg->Nlevels, &pc_gamg->Nlevels, NULL));
1541: {
1542: PetscReal eminmax[2] = {0., 0.};
1543: n = 2;
1544: PetscCall(PetscOptionsRealArray("-pc_gamg_eigenvalues", "extreme eigenvalues for smoothed aggregation", "PCGAMGSetEigenvalues", eminmax, &n, &flag));
1545: if (flag) {
1546: PetscCheck(n == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_gamg_eigenvalues: must specify 2 parameters, min and max eigenvalues");
1547: PetscCall(PCGAMGSetEigenvalues(pc, eminmax[1], eminmax[0]));
1548: }
1549: }
1550: /* set options for subtype */
1551: PetscCall((*pc_gamg->ops->setfromoptions)(pc, PetscOptionsObject));
1553: PetscCall(PCGetOptionsPrefix(pc, &pcpre));
1554: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : ""));
1555: PetscOptionsHeadEnd();
1556: PetscFunctionReturn(PETSC_SUCCESS);
1557: }
1559: /*MC
1560: PCGAMG - Geometric algebraic multigrid (AMG) preconditioner
1562: Options Database Keys:
1563: + -pc_gamg_type <type,default=agg> - one of agg, geo, or classical
1564: . -pc_gamg_repartition <bool,default=false> - repartition the degrees of freedom across the coarse grids as they are determined
1565: . -pc_gamg_asm_use_agg <bool,default=false> - use the aggregates from the coasening process to defined the subdomains on each level for the PCASM smoother
1566: . -pc_gamg_process_eq_limit <limit, default=50> - `PCGAMG` will reduce the number of MPI ranks used directly on the coarse grids so that there are around <limit>
1567: equations on each process that has degrees of freedom
1568: . -pc_gamg_coarse_eq_limit <limit, default=50> - Set maximum number of equations on coarsest grid to aim for.
1569: . -pc_gamg_reuse_interpolation <bool,default=true> - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations (should always be true)
1570: . -pc_gamg_threshold[] <thresh,default=[-1,...]> - Before aggregating the graph `PCGAMG` will remove small values from the graph on each level (< 0 does no filtering)
1571: - -pc_gamg_threshold_scale <scale,default=1> - Scaling of threshold on each coarser grid if not specified
1573: Options Database Keys for Aggregation:
1574: + -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation
1575: . -pc_gamg_square_graph <n,default=1> - alias for pc_gamg_aggressive_coarsening (deprecated)
1576: - -pc_gamg_aggressive_coarsening <n,default=1> - number of aggressive coarsening (MIS-2) levels from finest.
1578: Options Database Keys for Multigrid:
1579: + -pc_mg_cycle_type <v> - v or w, see `PCMGSetCycleType()`
1580: . -pc_mg_distinct_smoothup - configure the up and down (pre and post) smoothers separately, see PCMGSetDistinctSmoothUp()
1581: . -pc_mg_type <multiplicative> - (one of) additive multiplicative full kascade
1582: - -pc_mg_levels <levels> - Number of levels of multigrid to use. GAMG has a heuristic so pc_mg_levels is not usually used with GAMG
1584: Level: intermediate
1586: Notes:
1587: To obtain good performance for `PCGAMG` for vector valued problems you must
1588: call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point
1589: call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator
1591: See [the Users Manual section on PCGAMG](sec_amg) for more details.
1593: .seealso: `PCCreate()`, `PCSetType()`, `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`, `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`,
1594: `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`, `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetUseParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`, `PCGAMGGetType()`, `PCGAMGSetReuseInterpolation()`, `PCGAMGSetUseSAEstEig()`
1595: M*/
1597: PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc)
1598: {
1599: PC_GAMG *pc_gamg;
1600: PC_MG *mg;
1602: PetscFunctionBegin;
1603: /* register AMG type */
1604: PetscCall(PCGAMGInitializePackage());
1606: /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */
1607: PetscCall(PCSetType(pc, PCMG));
1608: PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCGAMG));
1610: /* create a supporting struct and attach it to pc */
1611: PetscCall(PetscNew(&pc_gamg));
1612: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
1613: mg = (PC_MG *)pc->data;
1614: mg->innerctx = pc_gamg;
1616: PetscCall(PetscNew(&pc_gamg->ops));
1618: /* these should be in subctx but repartitioning needs simple arrays */
1619: pc_gamg->data_sz = 0;
1620: pc_gamg->data = NULL;
1622: /* overwrite the pointers of PCMG by the functions of base class PCGAMG */
1623: pc->ops->setfromoptions = PCSetFromOptions_GAMG;
1624: pc->ops->setup = PCSetUp_GAMG;
1625: pc->ops->reset = PCReset_GAMG;
1626: pc->ops->destroy = PCDestroy_GAMG;
1627: mg->view = PCView_GAMG;
1629: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", PCGAMGSetProcEqLim_GAMG));
1630: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", PCGAMGSetCoarseEqLim_GAMG));
1631: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", PCGAMGSetRepartition_GAMG));
1632: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", PCGAMGSetEigenvalues_GAMG));
1633: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", PCGAMGSetUseSAEstEig_GAMG));
1634: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", PCGAMGSetReuseInterpolation_GAMG));
1635: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", PCGAMGASMSetUseAggs_GAMG));
1636: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseParallelCoarseGridSolve_C", PCGAMGSetUseParallelCoarseGridSolve_GAMG));
1637: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", PCGAMGSetCpuPinCoarseGrids_GAMG));
1638: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", PCGAMGSetCoarseGridLayoutType_GAMG));
1639: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", PCGAMGSetThreshold_GAMG));
1640: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", PCGAMGSetRankReductionFactors_GAMG));
1641: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", PCGAMGSetThresholdScale_GAMG));
1642: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", PCGAMGSetType_GAMG));
1643: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", PCGAMGGetType_GAMG));
1644: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", PCGAMGSetNlevels_GAMG));
1645: pc_gamg->repart = PETSC_FALSE;
1646: pc_gamg->reuse_prol = PETSC_TRUE;
1647: pc_gamg->use_aggs_in_asm = PETSC_FALSE;
1648: pc_gamg->use_parallel_coarse_grid_solver = PETSC_FALSE;
1649: pc_gamg->cpu_pin_coarse_grids = PETSC_FALSE;
1650: pc_gamg->layout_type = PCGAMG_LAYOUT_SPREAD;
1651: pc_gamg->min_eq_proc = 50;
1652: pc_gamg->coarse_eq_limit = 50;
1653: for (int i = 0; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = -1;
1654: pc_gamg->threshold_scale = 1.;
1655: pc_gamg->Nlevels = PETSC_MG_MAXLEVELS;
1656: pc_gamg->current_level = 0; /* don't need to init really */
1657: pc_gamg->use_sa_esteig = PETSC_TRUE;
1658: pc_gamg->emin = 0;
1659: pc_gamg->emax = 0;
1661: pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1663: /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */
1664: PetscCall(PCGAMGSetType(pc, PCGAMGAGG));
1665: PetscFunctionReturn(PETSC_SUCCESS);
1666: }
1668: /*@C
1669: PCGAMGInitializePackage - This function initializes everything in the `PCGAMG` package. It is called
1670: from `PCInitializePackage()`.
1672: Level: developer
1674: .seealso: `PetscInitialize()`
1675: @*/
1676: PetscErrorCode PCGAMGInitializePackage(void)
1677: {
1678: PetscInt l;
1680: PetscFunctionBegin;
1681: if (PCGAMGPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
1682: PCGAMGPackageInitialized = PETSC_TRUE;
1683: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGGEO, PCCreateGAMG_GEO));
1684: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGAGG, PCCreateGAMG_AGG));
1685: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGCLASSICAL, PCCreateGAMG_Classical));
1686: PetscCall(PetscRegisterFinalize(PCGAMGFinalizePackage));
1688: /* general events */
1689: PetscCall(PetscLogEventRegister("PCSetUp_GAMG+", PC_CLASSID, &petsc_gamg_setup_events[GAMG_SETUP]));
1690: PetscCall(PetscLogEventRegister(" PCGAMGCreateG", PC_CLASSID, &petsc_gamg_setup_events[GAMG_GRAPH]));
1691: PetscCall(PetscLogEventRegister(" GAMG Coarsen", PC_CLASSID, &petsc_gamg_setup_events[GAMG_COARSEN]));
1692: PetscCall(PetscLogEventRegister(" GAMG MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[GAMG_MIS]));
1693: PetscCall(PetscLogEventRegister(" PCGAMGProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROL]));
1694: PetscCall(PetscLogEventRegister(" GAMG Prol-col", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLA]));
1695: PetscCall(PetscLogEventRegister(" GAMG Prol-lift", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLB]));
1696: PetscCall(PetscLogEventRegister(" PCGAMGOptProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPT]));
1697: PetscCall(PetscLogEventRegister(" GAMG smooth", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPTSM]));
1698: PetscCall(PetscLogEventRegister(" PCGAMGCreateL", PC_CLASSID, &petsc_gamg_setup_events[GAMG_LEVEL]));
1699: PetscCall(PetscLogEventRegister(" GAMG PtAP", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PTAP]));
1700: PetscCall(PetscLogEventRegister(" GAMG Reduce", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REDUCE]));
1701: PetscCall(PetscLogEventRegister(" GAMG Repart", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REPART]));
1702: /* PetscCall(PetscLogEventRegister(" GAMG Inv-Srt", PC_CLASSID, &petsc_gamg_setup_events[SET13])); */
1703: /* PetscCall(PetscLogEventRegister(" GAMG Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14])); */
1704: /* PetscCall(PetscLogEventRegister(" GAMG Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15])); */
1705: for (l = 0; l < PETSC_MG_MAXLEVELS; l++) {
1706: char ename[32];
1708: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Squ l%02" PetscInt_FMT, l));
1709: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][0]));
1710: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Gal l%02" PetscInt_FMT, l));
1711: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][1]));
1712: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Opt l%02" PetscInt_FMT, l));
1713: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][2]));
1714: }
1715: #if defined(GAMG_STAGES)
1716: { /* create timer stages */
1717: char str[32];
1718: PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", 0));
1719: PetscCall(PetscLogStageRegister(str, &gamg_stages[0]));
1720: }
1721: #endif
1722: PetscFunctionReturn(PETSC_SUCCESS);
1723: }
1725: /*@C
1726: PCGAMGFinalizePackage - This function frees everything from the `PCGAMG` package. It is
1727: called from `PetscFinalize()` automatically.
1729: Level: developer
1731: .seealso: `PetscFinalize()`
1732: @*/
1733: PetscErrorCode PCGAMGFinalizePackage(void)
1734: {
1735: PetscFunctionBegin;
1736: PCGAMGPackageInitialized = PETSC_FALSE;
1737: PetscCall(PetscFunctionListDestroy(&GAMGList));
1738: PetscFunctionReturn(PETSC_SUCCESS);
1739: }
1741: /*@C
1742: PCGAMGRegister - Register a `PCGAMG` implementation.
1744: Input Parameters:
1745: + type - string that will be used as the name of the `PCGAMG` type.
1746: - create - function for creating the gamg context.
1748: Level: developer
1750: .seealso: `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()`
1751: @*/
1752: PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC))
1753: {
1754: PetscFunctionBegin;
1755: PetscCall(PCGAMGInitializePackage());
1756: PetscCall(PetscFunctionListAdd(&GAMGList, type, create));
1757: PetscFunctionReturn(PETSC_SUCCESS);
1758: }
1760: /*@C
1761: PCGAMGCreateGraph - Creates a graph that is used by the ``PCGAMGType`` in the coarsening process
1763: Input Parameters:
1764: + pc - the `PCGAMG`
1765: - A - the matrix, for any level
1767: Output Parameter:
1768: . G - the graph
1770: Level: advanced
1772: .seealso: `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()`
1773: @*/
1774: PetscErrorCode PCGAMGCreateGraph(PC pc, Mat A, Mat *G)
1775: {
1776: PC_MG *mg = (PC_MG *)pc->data;
1777: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1779: PetscFunctionBegin;
1780: PetscCall(pc_gamg->ops->creategraph(pc, A, G));
1781: PetscFunctionReturn(PETSC_SUCCESS);
1782: }