Actual source code: gamg.c
1: /*
2: GAMG geometric-algebric multigrid PC - Mark Adams 2011
3: */
4: #include <../src/ksp/pc/impls/gamg/gamg.h>
5: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>
7: #if defined(PETSC_HAVE_CUDA)
8: #include <petscdevice_cuda.h>
9: #endif
11: #if defined(PETSC_HAVE_HIP)
12: #include <petscdevice_hip.h>
13: #endif
15: PetscLogEvent petsc_gamg_setup_events[GAMG_NUM_SET];
16: PetscLogEvent petsc_gamg_setup_matmat_events[PETSC_MG_MAXLEVELS][3];
18: // #define GAMG_STAGES
19: #if defined(GAMG_STAGES)
20: static PetscLogStage gamg_stages[PETSC_MG_MAXLEVELS];
21: #endif
23: static PetscFunctionList GAMGList = NULL;
24: static PetscBool PCGAMGPackageInitialized;
26: static PetscErrorCode PCReset_GAMG(PC pc)
27: {
28: PC_MG *mg = (PC_MG *)pc->data;
29: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
31: PetscFunctionBegin;
32: PetscCall(PetscFree(pc_gamg->data));
33: pc_gamg->data_sz = 0;
34: PetscCall(PetscFree(pc_gamg->orig_data));
35: for (PetscInt level = 0; level < PETSC_MG_MAXLEVELS; level++) {
36: mg->min_eigen_DinvA[level] = 0;
37: mg->max_eigen_DinvA[level] = 0;
38: }
39: pc_gamg->emin = 0;
40: pc_gamg->emax = 0;
41: PetscCall(PCReset_MG(pc));
42: PetscCall(MatCoarsenDestroy(&pc_gamg->asm_crs));
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 = NULL, 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));
76: if (Pcolumnperm) *Pcolumnperm = NULL;
78: /* set 'ncrs' (nodes), 'ncrs_eq' (equations)*/
79: PetscCall(MatGetLocalSize(Pold, NULL, &ncrs_eq));
80: if (pc_gamg->data_cell_rows > 0) {
81: ncrs = pc_gamg->data_sz / pc_gamg->data_cell_cols / pc_gamg->data_cell_rows;
82: } else {
83: PetscInt bs;
84: PetscCall(MatGetBlockSizes(Pold, NULL, &bs));
85: ncrs = ncrs_eq / bs;
86: }
87: /* get number of PEs to make active 'new_size', reduce, can be any integer 1-P */
88: 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. */
89: #if defined(PETSC_HAVE_CUDA)
90: PetscShmComm pshmcomm;
91: PetscMPIInt locrank;
92: MPI_Comm loccomm;
93: PetscInt s_nnodes, r_nnodes, new_new_size;
94: cudaError_t cerr;
95: int devCount;
96: PetscCall(PetscShmCommGet(comm, &pshmcomm));
97: PetscCall(PetscShmCommGetMpiShmComm(pshmcomm, &loccomm));
98: PetscCallMPI(MPI_Comm_rank(loccomm, &locrank));
99: s_nnodes = !locrank;
100: PetscCallMPI(MPIU_Allreduce(&s_nnodes, &r_nnodes, 1, MPIU_INT, MPI_SUM, comm));
101: PetscCheck((size % r_nnodes) == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "odd number of nodes np=%d nnodes%" PetscInt_FMT, size, r_nnodes);
102: devCount = 0;
103: cerr = cudaGetDeviceCount(&devCount);
104: cudaGetLastError(); /* Reset the last error */
105: if (cerr == cudaSuccess && devCount >= 1) { /* There are devices, else go to heuristic */
106: new_new_size = r_nnodes * devCount;
107: new_size = new_new_size;
108: 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));
109: } else {
110: PetscCall(PetscInfo(pc, "%s: With Cuda but no device. Use heuristics.\n", ((PetscObject)pc)->prefix));
111: goto HEURISTIC;
112: }
113: #else
114: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "should not be here");
115: #endif
116: } else if (pc_gamg->level_reduction_factors[pc_gamg->current_level] > 0) {
117: if (nactive < pc_gamg->level_reduction_factors[pc_gamg->current_level]) {
118: new_size = 1;
119: PetscCall(PetscInfo(pc, "%s: reduction factor too small for %d active processes: reduce to one process\n", ((PetscObject)pc)->prefix, new_size));
120: } else {
121: 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]);
122: new_size = nactive / pc_gamg->level_reduction_factors[pc_gamg->current_level];
123: 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]));
124: }
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(Pold, NULL, &ncrs_eq_glob));
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: /* output - no repartitioning or reduction - could bail here
141: we know that the grid structure can be reused in MatPtAP */
142: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
143: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
144: PetscCall(MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, a_Amat_crs));
145: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
146: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
147: if (new_size < size) {
148: /* odd case where multiple coarse grids are on one processor or no coarsening ... */
149: PetscCall(PetscInfo(pc, "%s: reduced grid using same number of processors (%d) as last grid (use larger coarse grid)\n", ((PetscObject)pc)->prefix, nactive));
150: if (pc_gamg->cpu_pin_coarse_grids) {
151: PetscCall(MatBindToCPU(*a_Amat_crs, PETSC_TRUE));
152: PetscCall(MatBindToCPU(*a_P_inout, PETSC_TRUE));
153: }
154: }
155: } else { /* reduce active processors - we know that the grid structure can NOT be reused in MatPtAP */
156: PetscInt *counts, *newproc_idx, ii, jj, kk, strideNew, *tidx, ncrs_new, ncrs_eq_new, nloc_old, expand_factor = 1, rfactor = 1;
157: IS is_eq_newproc, is_eq_num, new_eq_indices;
158: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
159: nloc_old = ncrs_eq / cr_bs;
160: 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);
161: /* get new_size and rfactor */
162: if (pc_gamg->layout_type == PCGAMG_LAYOUT_SPREAD || !pc_gamg->repart) {
163: /* find factor */
164: if (new_size == 1) rfactor = size; /* don't modify */
165: else {
166: PetscReal best_fact = 0.;
167: jj = -1;
168: for (kk = 1; kk <= size; kk++) {
169: if (!(size % kk)) { /* a candidate */
170: PetscReal nactpe = (PetscReal)size / (PetscReal)kk, fact = nactpe / (PetscReal)new_size;
171: if (fact > 1.0) fact = 1. / fact; /* keep fact < 1 */
172: if (fact > best_fact) {
173: best_fact = fact;
174: jj = kk;
175: }
176: }
177: }
178: if (jj != -1) rfactor = jj;
179: else rfactor = 1; /* a prime */
180: if (pc_gamg->layout_type == PCGAMG_LAYOUT_COMPACT) expand_factor = 1;
181: else expand_factor = rfactor;
182: }
183: new_size = size / rfactor; /* make new size one that is factor */
184: if (new_size == nactive) { /* no repartitioning or reduction, bail out because nested here (rare) */
185: 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));
186: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
187: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
188: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
189: PetscCall(MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, a_Amat_crs));
190: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
191: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
194: }
195: /* make 'is_eq_newproc' */
196: if (pc_gamg->repart) { /* Repartition Cmat_{k} and move columns of P^{k}_{k-1} and coordinates of primal part accordingly */
197: Mat adj;
199: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
200: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
201: PetscCall(MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat));
202: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
203: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
204: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_REPART], 0, 0, 0, 0));
205: 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"));
206: /* get 'adj' */
207: if (cr_bs == 1) {
208: PetscCall(MatConvert(Cmat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj));
209: } else {
210: /* make a scalar matrix to partition (no Stokes here) */
211: Mat tMat;
212: PetscInt Istart_crs, Iend_crs, ncols, jj, Ii;
213: const PetscScalar *vals;
214: const PetscInt *idx;
215: PetscInt *d_nnz, *o_nnz, M, N, maxnnz = 0, *j_buf = NULL;
216: PetscScalar *v_buff = NULL;
217: static PetscInt llev = 0; /* ugly but just used for debugging */
218: MatType mtype;
220: PetscCall(PetscMalloc2(ncrs, &d_nnz, ncrs, &o_nnz));
221: PetscCall(MatGetOwnershipRange(Cmat, &Istart_crs, &Iend_crs));
222: PetscCall(MatGetSize(Cmat, &M, &N));
223: for (Ii = Istart_crs, jj = 0; Ii < Iend_crs; Ii += cr_bs, jj++) {
224: PetscCall(MatGetRow(Cmat, Ii, &ncols, NULL, NULL));
225: d_nnz[jj] = ncols / cr_bs;
226: o_nnz[jj] = ncols / cr_bs;
227: if (ncols > maxnnz) maxnnz = ncols;
228: PetscCall(MatRestoreRow(Cmat, Ii, &ncols, NULL, NULL));
229: if (d_nnz[jj] > ncrs) d_nnz[jj] = ncrs;
230: if (o_nnz[jj] > (M / cr_bs - ncrs)) o_nnz[jj] = M / cr_bs - ncrs;
231: }
233: PetscCall(MatGetType(Amat_fine, &mtype));
234: PetscCall(MatCreate(comm, &tMat));
235: PetscCall(MatSetSizes(tMat, ncrs, ncrs, PETSC_DETERMINE, PETSC_DETERMINE));
236: PetscCall(MatSetType(tMat, mtype));
237: PetscCall(MatSeqAIJSetPreallocation(tMat, 0, d_nnz));
238: PetscCall(MatMPIAIJSetPreallocation(tMat, 0, d_nnz, 0, o_nnz));
239: PetscCall(PetscFree2(d_nnz, o_nnz));
240: PetscCall(PetscMalloc2(maxnnz, &j_buf, maxnnz, &v_buff));
241: for (ii = 0; ii < maxnnz; ii++) v_buff[ii] = 1.;
243: for (ii = Istart_crs; ii < Iend_crs; ii++) {
244: PetscInt dest_row = ii / cr_bs;
245: PetscCall(MatGetRow(Cmat, ii, &ncols, &idx, &vals));
246: for (jj = 0; jj < ncols; jj++) j_buf[jj] = idx[jj] / cr_bs;
247: PetscCall(MatSetValues(tMat, 1, &dest_row, ncols, j_buf, v_buff, ADD_VALUES));
248: PetscCall(MatRestoreRow(Cmat, ii, &ncols, &idx, &vals));
249: }
250: PetscCall(MatAssemblyBegin(tMat, MAT_FINAL_ASSEMBLY));
251: PetscCall(MatAssemblyEnd(tMat, MAT_FINAL_ASSEMBLY));
252: PetscCall(PetscFree2(j_buf, v_buff));
254: if (llev++ == -1) {
255: PetscViewer viewer;
256: char fname[32];
257: PetscCall(PetscSNPrintf(fname, sizeof(fname), "part_mat_%" PetscInt_FMT ".mat", llev));
258: PetscCall(PetscViewerBinaryOpen(comm, fname, FILE_MODE_WRITE, &viewer));
259: PetscCall(MatView(tMat, viewer));
260: PetscCall(PetscViewerDestroy(&viewer));
261: }
262: PetscCall(MatConvert(tMat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj));
263: PetscCall(MatDestroy(&tMat));
264: } /* create 'adj' */
266: { /* partition: get newproc_idx */
267: char prefix[256];
268: const char *pcpre;
269: const PetscInt *is_idx;
270: MatPartitioning mpart;
271: IS proc_is;
273: PetscCall(MatPartitioningCreate(comm, &mpart));
274: PetscCall(MatPartitioningSetAdjacency(mpart, adj));
275: PetscCall(PCGetOptionsPrefix(pc, &pcpre));
276: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : ""));
277: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix));
278: PetscCall(MatPartitioningSetFromOptions(mpart));
279: PetscCall(MatPartitioningSetNParts(mpart, new_size));
280: PetscCall(MatPartitioningApply(mpart, &proc_is));
281: PetscCall(MatPartitioningDestroy(&mpart));
283: /* collect IS info */
284: PetscCall(PetscMalloc1(ncrs_eq, &newproc_idx));
285: PetscCall(ISGetIndices(proc_is, &is_idx));
286: for (kk = jj = 0; kk < nloc_old; kk++) {
287: for (ii = 0; ii < cr_bs; ii++, jj++) newproc_idx[jj] = is_idx[kk] * expand_factor; /* distribution */
288: }
289: PetscCall(ISRestoreIndices(proc_is, &is_idx));
290: PetscCall(ISDestroy(&proc_is));
291: }
292: PetscCall(MatDestroy(&adj));
294: PetscCall(ISCreateGeneral(comm, ncrs_eq, newproc_idx, PETSC_OWN_POINTER, &is_eq_newproc));
295: /*
296: Create an index set from the is_eq_newproc index set to indicate the mapping TO
297: */
298: PetscCall(ISPartitioningToNumbering(is_eq_newproc, &is_eq_num));
299: /*
300: Determine how many equations/vertices are assigned to each processor
301: */
302: PetscCall(PetscMalloc1(size, &counts));
303: PetscCall(ISPartitioningCount(is_eq_newproc, size, counts));
304: ncrs_eq_new = counts[rank];
305: PetscCall(ISDestroy(&is_eq_newproc));
306: PetscCall(PetscFree(counts));
307: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REPART], 0, 0, 0, 0));
308: } else { /* simple aggregation of parts -- 'is_eq_newproc' */
309: const PetscInt *ranges;
310: PetscInt newstart = 0;
311: PetscLayout ilay;
313: PetscCheck(new_size != nactive, PETSC_COMM_SELF, PETSC_ERR_PLIB, "new_size==nactive. Should not happen");
314: PetscCall(PetscInfo(pc, "%s: Number of equations (loc) %" PetscInt_FMT " with simple aggregation\n", ((PetscObject)pc)->prefix, ncrs_eq));
315: PetscCallMPI(MPI_Exscan(&ncrs_eq, &newstart, 1, MPIU_INT, MPI_SUM, comm));
316: PetscCall(ISCreateStride(comm, ncrs_eq, newstart, 1, &is_eq_num));
317: PetscCall(ISSetPermutation(is_eq_num));
318: PetscCall(ISGetLayout(is_eq_num, &ilay));
319: PetscCall(PetscLayoutGetRanges(ilay, &ranges));
320: ncrs_eq_new = 0;
321: for (PetscInt r = 0; r < size; r++)
322: if (rank == (r / rfactor) * expand_factor) ncrs_eq_new += ranges[r + 1] - ranges[r];
323: //targetPE = (rank / rfactor) * expand_factor;
324: //PetscCall(ISCreateStride(comm, ncrs_eq, targetPE, 0, &is_eq_newproc));
325: //PetscCall(ISPartitioningToNumbering(is_eq_newproc, &is_eq_num));
326: //PetscCall(PetscMalloc1(size, &counts));
327: //PetscCall(ISPartitioningCount(is_eq_newproc, size, counts));
328: //ncrs_eq_new = counts[rank];
329: //PetscCall(ISDestroy(&is_eq_newproc));
330: //PetscCall(PetscFree(counts));
331: } /* end simple 'is_eq_newproc' */
333: ncrs_new = ncrs_eq_new / cr_bs;
335: /* 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 */
336: {
337: Vec src_crd, dest_crd;
338: const PetscInt *idx, ndata_rows = pc_gamg->data_cell_rows, ndata_cols = pc_gamg->data_cell_cols, node_data_sz = ndata_rows * ndata_cols;
339: VecScatter vecscat;
340: PetscScalar *array;
341: IS isscat;
342: /* move data (for primal equations only) */
343: /* Create a vector to contain the newly ordered element information */
344: PetscCall(VecCreate(comm, &dest_crd));
345: PetscCall(VecSetSizes(dest_crd, node_data_sz * ncrs_new, PETSC_DECIDE));
346: PetscCall(VecSetType(dest_crd, VECSTANDARD)); /* this is needed! */
347: /*
348: There are 'ndata_rows*ndata_cols' data items per node, (one can think of the vectors of having
349: a block size of ...). Note, ISs are expanded into equation space by 'cr_bs'.
350: */
351: PetscCall(PetscMalloc1(ncrs * node_data_sz, &tidx));
352: PetscCall(ISGetIndices(is_eq_num, &idx));
353: for (ii = 0, jj = 0; ii < ncrs; ii++) {
354: PetscInt id = idx[ii * cr_bs] / cr_bs; /* get node back */
355: for (kk = 0; kk < node_data_sz; kk++, jj++) tidx[jj] = id * node_data_sz + kk;
356: }
357: PetscCall(ISRestoreIndices(is_eq_num, &idx));
358: PetscCall(ISCreateGeneral(comm, node_data_sz * ncrs, tidx, PETSC_COPY_VALUES, &isscat));
359: PetscCall(PetscFree(tidx));
360: /*
361: Create a vector to contain the original vertex information for each element
362: */
363: PetscCall(VecCreateSeq(PETSC_COMM_SELF, node_data_sz * ncrs, &src_crd));
364: for (jj = 0; jj < ndata_cols; jj++) {
365: const PetscInt stride0 = ncrs * pc_gamg->data_cell_rows;
366: for (ii = 0; ii < ncrs; ii++) {
367: for (kk = 0; kk < ndata_rows; kk++) {
368: PetscInt ix = ii * ndata_rows + kk + jj * stride0, jx = ii * node_data_sz + kk * ndata_cols + jj;
369: PetscScalar tt = pc_gamg->data[ix];
370: PetscCall(VecSetValues(src_crd, 1, &jx, &tt, INSERT_VALUES));
371: }
372: }
373: }
374: PetscCall(VecAssemblyBegin(src_crd));
375: PetscCall(VecAssemblyEnd(src_crd));
376: /*
377: Scatter the element vertex information (still in the original vertex ordering)
378: to the correct processor
379: */
380: PetscCall(VecScatterCreate(src_crd, NULL, dest_crd, isscat, &vecscat));
381: PetscCall(ISDestroy(&isscat));
382: PetscCall(VecScatterBegin(vecscat, src_crd, dest_crd, INSERT_VALUES, SCATTER_FORWARD));
383: PetscCall(VecScatterEnd(vecscat, src_crd, dest_crd, INSERT_VALUES, SCATTER_FORWARD));
384: PetscCall(VecScatterDestroy(&vecscat));
385: PetscCall(VecDestroy(&src_crd));
386: /*
387: Put the element vertex data into a new allocation of the gdata->ele
388: */
389: PetscCall(PetscFree(pc_gamg->data));
390: PetscCall(PetscMalloc1(node_data_sz * ncrs_new, &pc_gamg->data));
392: pc_gamg->data_sz = node_data_sz * ncrs_new;
393: strideNew = ncrs_new * ndata_rows;
395: PetscCall(VecGetArray(dest_crd, &array));
396: for (jj = 0; jj < ndata_cols; jj++) {
397: for (ii = 0; ii < ncrs_new; ii++) {
398: for (kk = 0; kk < ndata_rows; kk++) {
399: PetscInt ix = ii * ndata_rows + kk + jj * strideNew, jx = ii * node_data_sz + kk * ndata_cols + jj;
400: pc_gamg->data[ix] = PetscRealPart(array[jx]);
401: }
402: }
403: }
404: PetscCall(VecRestoreArray(dest_crd, &array));
405: PetscCall(VecDestroy(&dest_crd));
406: }
407: /* move A and P (columns) with new layout */
408: /*
409: Invert for MatCreateSubMatrix
410: */
411: PetscCall(ISInvertPermutation(is_eq_num, ncrs_eq_new, &new_eq_indices));
412: PetscCall(ISSort(new_eq_indices));
413: PetscCall(ISSetBlockSize(new_eq_indices, cr_bs));
414: if (Pcolumnperm) {
415: PetscCall(PetscObjectReference((PetscObject)new_eq_indices));
416: *Pcolumnperm = new_eq_indices;
417: }
418: PetscCall(ISDestroy(&is_eq_num));
420: /* 'a_Amat_crs' output */
421: if (Cmat) { /* repartitioning from Cmat adjacency case */
422: Mat mat;
423: PetscBool isset, isspd, isher;
424: #if !defined(PETSC_USE_COMPLEX)
425: PetscBool issym;
426: #endif
428: PetscCall(MatCreateSubMatrix(Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat));
429: PetscCall(MatIsSPDKnown(Cmat, &isset, &isspd)); // like MatPropagateSymmetryOptions, but should set MAT_STRUCTURALLY_SYMMETRIC ?
430: if (isset) PetscCall(MatSetOption(mat, MAT_SPD, isspd));
431: else {
432: PetscCall(MatIsHermitianKnown(Cmat, &isset, &isher));
433: if (isset) PetscCall(MatSetOption(mat, MAT_HERMITIAN, isher));
434: else {
435: #if !defined(PETSC_USE_COMPLEX)
436: PetscCall(MatIsSymmetricKnown(Cmat, &isset, &issym));
437: if (isset) PetscCall(MatSetOption(mat, MAT_SYMMETRIC, issym));
438: #endif
439: }
440: }
441: *a_Amat_crs = mat;
442: }
444: /* prolongator */
445: {
446: IS findices;
447: PetscInt Istart, Iend;
448: Mat Pnew;
450: PetscCall(MatGetOwnershipRange(Pold, &Istart, &Iend));
451: PetscCall(ISCreateStride(comm, Iend - Istart, Istart, 1, &findices));
452: PetscCall(ISSetBlockSize(findices, f_bs));
453: PetscCall(MatCreateSubMatrix(Pold, findices, new_eq_indices, MAT_INITIAL_MATRIX, &Pnew));
454: PetscCall(ISDestroy(&findices));
455: PetscCall(MatSetOption(Pnew, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
457: PetscCall(MatDestroy(a_P_inout));
459: /* output - repartitioned */
460: *a_P_inout = Pnew;
461: }
463: if (!Cmat) { /* simple repartitioning case */
464: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
465: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
466: PetscCall(MatPtAP(Amat_fine, *a_P_inout, MAT_INITIAL_MATRIX, 2.0, a_Amat_crs));
467: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][1], 0, 0, 0, 0));
468: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PTAP], 0, 0, 0, 0));
469: }
470: PetscCall(MatDestroy(&Cmat));
471: PetscCall(ISDestroy(&new_eq_indices));
473: *a_nactive_proc = new_size; /* output */
475: /* pinning on reduced grids, not a bad heuristic and optimization gets folded into process reduction optimization */
476: if (pc_gamg->cpu_pin_coarse_grids) {
477: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
478: static PetscInt llev = 2;
479: PetscCall(PetscInfo(pc, "%s: Pinning level %" PetscInt_FMT " to the CPU\n", ((PetscObject)pc)->prefix, llev++));
480: #endif
481: PetscCall(MatBindToCPU(*a_Amat_crs, PETSC_TRUE));
482: PetscCall(MatBindToCPU(*a_P_inout, PETSC_TRUE));
483: if (1) { /* HACK: move this to MatBindCPU_MPIAIJXXX; lvec is created, need to pin it, this is done in MatSetUpMultiply_MPIAIJ. Hack */
484: Mat A = *a_Amat_crs, P = *a_P_inout;
485: PetscMPIInt size;
486: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
487: if (size > 1) {
488: Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data, *p = (Mat_MPIAIJ *)P->data;
489: PetscCall(VecBindToCPU(a->lvec, PETSC_TRUE));
490: PetscCall(VecBindToCPU(p->lvec, PETSC_TRUE));
491: }
492: }
493: }
494: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_REDUCE], 0, 0, 0, 0));
495: } // processor reduce
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /* Computes the symbolic part of Gmat1^T * Gmat1 */
500: PetscErrorCode PCGAMGSquareGraph_GAMG(PC a_pc, Mat Gmat1, Mat *Gmat2)
501: {
502: const char *prefix;
503: char addp[32];
504: PC_MG *mg = (PC_MG *)a_pc->data;
505: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
507: PetscFunctionBegin;
508: PetscCall(PCGetOptionsPrefix(a_pc, &prefix));
509: PetscCall(PetscInfo(a_pc, "%s: Square Graph on level %" PetscInt_FMT "\n", ((PetscObject)a_pc)->prefix, pc_gamg->current_level + 1));
510: PetscCall(MatProductCreate(Gmat1, Gmat1, NULL, Gmat2));
511: PetscCall(MatSetOptionsPrefix(*Gmat2, prefix));
512: PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_square_%" PetscInt_FMT "_", pc_gamg->current_level));
513: PetscCall(MatAppendOptionsPrefix(*Gmat2, addp));
514: if ((*Gmat2)->structurally_symmetric == PETSC_BOOL3_TRUE) {
515: PetscCall(MatProductSetType(*Gmat2, MATPRODUCT_AB));
516: } else {
517: PetscCall(MatSetOption(Gmat1, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
518: PetscCall(MatProductSetType(*Gmat2, MATPRODUCT_AtB));
519: }
520: PetscCall(MatProductSetFromOptions(*Gmat2));
521: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0], 0, 0, 0, 0));
522: PetscCall(MatProductSymbolic(*Gmat2));
523: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][0], 0, 0, 0, 0));
524: PetscCall(MatProductClear(*Gmat2));
525: /* we only need the sparsity, cheat and tell PETSc the matrix has been assembled */
526: (*Gmat2)->assembled = PETSC_TRUE;
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /*
531: PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner
532: by setting data structures and options.
534: Input Parameter:
535: . pc - the preconditioner context
537: */
538: static PetscErrorCode PCSetUp_GAMG(PC pc)
539: {
540: PC_MG *mg = (PC_MG *)pc->data;
541: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
542: Mat Pmat = pc->pmat;
543: PetscInt fine_level, level, level1, bs, M, N, qq, lidx, nASMBlocksArr[PETSC_MG_MAXLEVELS], cr_bs;
544: MPI_Comm comm;
545: PetscMPIInt rank, size, nactivepe;
546: Mat Aarr[PETSC_MG_MAXLEVELS], Parr[PETSC_MG_MAXLEVELS];
547: IS *ASMLocalIDsArr[PETSC_MG_MAXLEVELS];
548: PetscBool is_last = PETSC_FALSE;
549: #if defined(PETSC_USE_INFO)
550: PetscLogDouble nnz0 = 0., nnztot = 0.;
551: MatInfo info;
552: #endif
554: PetscFunctionBegin;
555: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
556: PetscCallMPI(MPI_Comm_rank(comm, &rank));
557: PetscCallMPI(MPI_Comm_size(comm, &size));
558: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
559: if (pc->setupcalled) {
560: if (!pc_gamg->reuse_prol || pc->flag == DIFFERENT_NONZERO_PATTERN) {
561: /* reset everything */
562: PetscCall(PCReset_MG(pc));
563: pc->setupcalled = PETSC_FALSE;
564: } else {
565: PC_MG_Levels **mglevels = mg->levels;
566: /* just do Galerkin grids */
567: Mat B, dA, dB;
568: if (pc_gamg->Nlevels > 1) {
569: PetscInt gl;
570: /* currently only handle case where mat and pmat are the same on coarser levels */
571: PetscCall(KSPGetOperators(mglevels[pc_gamg->Nlevels - 1]->smoothd, &dA, &dB));
572: /* (re)set to get dirty flag */
573: PetscCall(KSPSetOperators(mglevels[pc_gamg->Nlevels - 1]->smoothd, dA, dB));
575: for (level = pc_gamg->Nlevels - 2, gl = 0; level >= 0; level--, gl++) {
576: MatReuse reuse = MAT_INITIAL_MATRIX;
577: #if defined(GAMG_STAGES)
578: PetscCall(PetscLogStagePush(gamg_stages[gl]));
579: #endif
580: /* matrix nonzero structure can change from repartitioning or process reduction but don't know if we have process reduction here. Should fix */
581: PetscCall(KSPGetOperators(mglevels[level]->smoothd, NULL, &B));
582: if (B->product) {
583: if (B->product->A == dB && B->product->B == mglevels[level + 1]->interpolate) reuse = MAT_REUSE_MATRIX;
584: }
585: if (reuse == MAT_INITIAL_MATRIX) PetscCall(MatDestroy(&mglevels[level]->A));
586: if (reuse == MAT_REUSE_MATRIX) {
587: PetscCall(PetscInfo(pc, "%s: RAP after initial setup, reuse matrix level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
588: } else {
589: PetscCall(PetscInfo(pc, "%s: RAP after initial setup, with repartitioning (new matrix) level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
590: }
591: PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[gl][1], 0, 0, 0, 0));
592: PetscCall(MatPtAP(dB, mglevels[level + 1]->interpolate, reuse, PETSC_DETERMINE, &B));
593: PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[gl][1], 0, 0, 0, 0));
594: if (reuse == MAT_INITIAL_MATRIX) mglevels[level]->A = B;
595: PetscCall(KSPSetOperators(mglevels[level]->smoothd, B, B));
596: // check for redoing eigen estimates
597: if (pc_gamg->recompute_esteig) {
598: PetscBool ischeb;
599: KSP smoother;
600: PetscCall(PCMGGetSmoother(pc, level + 1, &smoother));
601: PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb));
602: if (ischeb) {
603: KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data;
604: cheb->emin_provided = 0;
605: cheb->emax_provided = 0;
606: }
607: /* we could call PetscCall(KSPChebyshevSetEigenvalues(smoother, 0, 0)); but the logic does not work properly */
608: }
609: // inc
610: dB = B;
611: #if defined(GAMG_STAGES)
612: PetscCall(PetscLogStagePop());
613: #endif
614: }
615: }
617: PetscCall(PCSetUp_MG(pc));
618: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
621: }
623: if (!pc_gamg->data) {
624: if (pc_gamg->orig_data) {
625: PetscCall(MatGetBlockSize(Pmat, &bs));
626: PetscCall(MatGetLocalSize(Pmat, &qq, NULL));
628: pc_gamg->data_sz = (qq / bs) * pc_gamg->orig_data_cell_rows * pc_gamg->orig_data_cell_cols;
629: pc_gamg->data_cell_rows = pc_gamg->orig_data_cell_rows;
630: pc_gamg->data_cell_cols = pc_gamg->orig_data_cell_cols;
632: PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data));
633: for (qq = 0; qq < pc_gamg->data_sz; qq++) pc_gamg->data[qq] = pc_gamg->orig_data[qq];
634: } else {
635: PetscCheck(pc_gamg->ops->createdefaultdata, comm, PETSC_ERR_PLIB, "'createdefaultdata' not set(?) need to support NULL data");
636: PetscCall(pc_gamg->ops->createdefaultdata(pc, Pmat));
637: }
638: }
640: /* cache original data for reuse */
641: if (!pc_gamg->orig_data && (PetscBool)(!pc_gamg->reuse_prol)) {
642: PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->orig_data));
643: for (qq = 0; qq < pc_gamg->data_sz; qq++) pc_gamg->orig_data[qq] = pc_gamg->data[qq];
644: pc_gamg->orig_data_cell_rows = pc_gamg->data_cell_rows;
645: pc_gamg->orig_data_cell_cols = pc_gamg->data_cell_cols;
646: }
648: /* get basic dims */
649: PetscCall(MatGetBlockSize(Pmat, &bs));
650: PetscCall(MatGetSize(Pmat, &M, NULL));
652: #if defined(PETSC_USE_INFO)
653: PetscCall(MatGetInfo(Pmat, MAT_GLOBAL_SUM, &info)); /* global reduction */
654: nnz0 = info.nz_used;
655: nnztot = info.nz_used;
656: #endif
657: PetscCall(PetscInfo(pc, "%s: level %d) N=%" PetscInt_FMT ", n data rows=%" PetscInt_FMT ", n data cols=%" PetscInt_FMT ", nnz/row (ave)=%" PetscInt_FMT ", block size %" PetscInt_FMT ", np=%d\n", ((PetscObject)pc)->prefix, 0, M, pc_gamg->data_cell_rows,
658: pc_gamg->data_cell_cols, (PetscInt)(nnz0 / (PetscReal)M + 0.5), bs, size));
660: /* Get A_i and R_i */
661: for (level = 0, Aarr[0] = Pmat, nactivepe = size; level < (pc_gamg->Nlevels - 1) && (level == 0 || M > pc_gamg->coarse_eq_limit); level++) {
662: pc_gamg->current_level = level;
663: level1 = level + 1;
664: #if defined(GAMG_STAGES)
665: if (!gamg_stages[level]) {
666: char str[32];
667: PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %" PetscInt_FMT, level));
668: PetscCall(PetscLogStageRegister(str, &gamg_stages[level]));
669: }
670: PetscCall(PetscLogStagePush(gamg_stages[level]));
671: #endif
672: /* construct prolongator - Parr[level1] */
673: if (level == 0 && pc_gamg->injection_index_size > 0) {
674: Mat Prol;
675: MatType mtype;
676: PetscInt prol_m, prol_n, Prol_N = (M / bs) * pc_gamg->injection_index_size, Istart, Iend, nn, row;
677: PetscCall(PetscInfo(pc, "Create fine grid injection space prolongation %" PetscInt_FMT " x %" PetscInt_FMT ". %s\n", M, Prol_N, pc_gamg->data ? "delete null space data" : ""));
678: PetscCall(MatGetOwnershipRange(Pmat, &Istart, &Iend));
679: PetscCall(MatGetLocalSize(Pmat, &prol_m, NULL)); // rows m x n
680: prol_n = (prol_m / bs) * pc_gamg->injection_index_size;
681: PetscCheck(pc_gamg->injection_index_size < bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Injection size %" PetscInt_FMT " must be less that block size %" PetscInt_FMT, pc_gamg->injection_index_size, bs);
682: PetscCall(MatGetType(Pmat, &mtype));
683: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &Prol));
684: PetscCall(MatSetBlockSizes(Prol, bs, pc_gamg->injection_index_size));
685: PetscCall(MatSetSizes(Prol, prol_m, prol_n, M, Prol_N));
686: PetscCall(MatSetType(Prol, mtype));
687: #if PetscDefined(HAVE_DEVICE)
688: PetscBool flg;
689: PetscCall(MatBoundToCPU(Pmat, &flg));
690: PetscCall(MatBindToCPU(Prol, flg));
691: if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE));
692: #endif
693: PetscCall(MatSeqAIJSetPreallocation(Prol, 1, NULL));
694: PetscCall(MatMPIAIJSetPreallocation(Prol, 1, NULL, 0, NULL));
695: // set I \kron [1, 1, ... ]^T
696: for (PetscInt ii = Istart, col = (Istart / bs) * pc_gamg->injection_index_size; ii < Iend; ii += bs) {
697: const PetscScalar one = 1;
698: for (PetscInt jj = 0; jj < pc_gamg->injection_index_size; jj++, col++) {
699: PetscInt row = ii + pc_gamg->injection_index[jj];
700: PetscCall(MatSetValues(Prol, 1, &row, 1, &col, &one, INSERT_VALUES));
701: }
702: }
703: PetscCall(MatAssemblyBegin(Prol, MAT_FINAL_ASSEMBLY));
704: PetscCall(MatAssemblyEnd(Prol, MAT_FINAL_ASSEMBLY));
705: PetscCall(MatViewFromOptions(Prol, NULL, "-mat_view_injection"));
706: PetscCall(MatGetBlockSizes(Prol, NULL, &cr_bs)); // column size
707: Parr[level1] = Prol;
708: // can not deal with null space -- with array of 'injection cols' we could take 'injection rows and 'injection cols' to 'data'
709: if (pc_gamg->data) {
710: pc_gamg->data_cell_cols = pc_gamg->injection_index_size;
711: pc_gamg->data_cell_rows = pc_gamg->injection_index_size;
712: pc_gamg->orig_data_cell_cols = 0;
713: pc_gamg->orig_data_cell_rows = 0;
714: PetscCall(PetscFree(pc_gamg->data));
715: pc_gamg->data_sz = pc_gamg->injection_index_size * prol_n;
716: PetscCall(PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data));
717: for (row = nn = 0; row < prol_n; row += pc_gamg->injection_index_size) {
718: for (PetscInt jj = 0; jj < pc_gamg->injection_index_size; jj++) {
719: PetscInt idx = row * pc_gamg->injection_index_size + jj * pc_gamg->injection_index_size;
720: for (PetscInt kk = 0; kk < pc_gamg->injection_index_size; kk++, nn++) pc_gamg->data[idx + kk] = (jj == kk) ? 1 : 0;
721: }
722: }
723: PetscCheck(nn == pc_gamg->data_sz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nn != pc_gamg->data_sz %" PetscInt_FMT " %" PetscInt_FMT, pc_gamg->data_sz, nn);
724: }
725: } else {
726: Mat Gmat, mat;
727: PetscCoarsenData *agg_lists;
728: Mat Prol11;
730: PetscCall(PCGAMGCreateGraph(pc, Aarr[level], &Gmat));
731: PetscCall(pc_gamg->ops->coarsen(pc, &Gmat, &agg_lists)); // Gmat may have ghosts for QR aggregates not in matrix
732: PetscCall(PetscCDGetMat(agg_lists, &mat));
733: if (!mat) PetscCall(PetscCDSetMat(agg_lists, Gmat));
734: PetscCall(pc_gamg->ops->prolongator(pc, Aarr[level], agg_lists, &Prol11));
735: /* could have failed to create new level */
736: if (Prol11) {
737: const char *prefix;
738: char addp[32];
740: /* get new block size of coarse matrices */
741: PetscCall(MatGetBlockSizes(Prol11, NULL, &cr_bs)); // column size
743: if (pc_gamg->ops->optprolongator) {
744: /* smooth */
745: PetscCall(pc_gamg->ops->optprolongator(pc, Aarr[level], &Prol11));
746: }
748: if (pc_gamg->use_aggs_in_asm) {
749: PetscInt bs;
750: PetscCall(MatGetBlockSizes(Prol11, &bs, NULL)); // row block size
751: PetscCall(PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level]));
752: PetscCall(PetscInfo(pc, "%" PetscInt_FMT ": %" PetscInt_FMT " ASM local domains, bs = %" PetscInt_FMT "\n", level, nASMBlocksArr[level], bs));
753: } else if (pc_gamg->asm_hem_aggs) {
754: const char *prefix;
755: PetscInt bs;
757: /*
758: Do not use aggs created for defining coarser problems, instead create aggs specifically to use
759: to define PCASM blocks.
760: */
761: PetscCall(PetscCDGetMat(agg_lists, &mat));
762: if (mat == Gmat) PetscCall(PetscCDClearMat(agg_lists)); // take the Mat away from the list (yuck)
763: PetscCall(PetscCDDestroy(agg_lists));
764: PetscCall(PetscInfo(pc, "HEM ASM passes = %" PetscInt_FMT "\n", pc_gamg->asm_hem_aggs));
765: PetscCall(MatCoarsenDestroy(&pc_gamg->asm_crs));
766: PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg->asm_crs));
767: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix));
768: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg->asm_crs, prefix));
769: PetscCall(MatCoarsenSetFromOptions(pc_gamg->asm_crs)); // get strength args
770: PetscCall(MatCoarsenSetType(pc_gamg->asm_crs, MATCOARSENHEM));
771: PetscCall(MatCoarsenSetMaximumIterations(pc_gamg->asm_crs, pc_gamg->asm_hem_aggs));
772: PetscCall(MatCoarsenSetAdjacency(pc_gamg->asm_crs, Gmat));
773: PetscCall(MatCoarsenSetStrictAggs(pc_gamg->asm_crs, PETSC_TRUE));
774: PetscCall(MatCoarsenApply(pc_gamg->asm_crs));
775: PetscCall(MatCoarsenGetData(pc_gamg->asm_crs, &agg_lists)); /* output */
776: // create aggregates
777: PetscCall(MatGetBlockSizes(Aarr[level], &bs, NULL)); // row block size
778: PetscCall(PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level]));
779: }
780: PetscCall(PCGetOptionsPrefix(pc, &prefix));
781: PetscCall(MatSetOptionsPrefix(Prol11, prefix));
782: PetscCall(PetscSNPrintf(addp, sizeof(addp), "pc_gamg_prolongator_%" PetscInt_FMT "_", level));
783: PetscCall(MatAppendOptionsPrefix(Prol11, addp));
784: /* Always generate the transpose with CUDA
785: Such behaviour can be adapted with -pc_gamg_prolongator_ prefixed options */
786: PetscCall(MatSetOption(Prol11, MAT_FORM_EXPLICIT_TRANSPOSE, PETSC_TRUE));
787: PetscCall(MatSetFromOptions(Prol11));
788: Parr[level1] = Prol11;
789: } else Parr[level1] = NULL; /* failed to coarsen */
790: PetscCall(PetscCDGetMat(agg_lists, &mat));
791: if (mat == Gmat) PetscCall(PetscCDClearMat(agg_lists)); // take the Mat away from the list (yuck)
792: PetscCall(MatDestroy(&Gmat));
793: PetscCall(PetscCDDestroy(agg_lists));
794: } /* construct prolongator scope */
795: if (level == 0) Aarr[0] = Pmat; /* use Pmat for finest level setup */
796: if (!Parr[level1]) { /* failed to coarsen */
797: PetscCall(PetscInfo(pc, "%s: Stop gridding, level %" PetscInt_FMT "\n", ((PetscObject)pc)->prefix, level));
798: #if defined(GAMG_STAGES)
799: PetscCall(PetscLogStagePop());
800: #endif
801: break;
802: }
803: PetscCall(MatGetSize(Parr[level1], &M, &N)); /* N is next M, a loop test variables */
804: PetscCheck(!is_last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Is last ?");
805: if (N <= pc_gamg->coarse_eq_limit) is_last = PETSC_TRUE;
806: if (level1 == pc_gamg->Nlevels - 1) is_last = PETSC_TRUE;
807: if (level == PETSC_MG_MAXLEVELS - 2) is_last = PETSC_TRUE;
808: PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0));
809: PetscCall(pc_gamg->ops->createlevel(pc, Aarr[level], cr_bs, &Parr[level1], &Aarr[level1], &nactivepe, NULL, is_last));
810: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_LEVEL], 0, 0, 0, 0));
812: PetscCall(MatGetSize(Aarr[level1], &M, &N)); /* M is loop test variables */
813: #if defined(PETSC_USE_INFO)
814: PetscCall(MatGetInfo(Aarr[level1], MAT_GLOBAL_SUM, &info));
815: nnztot += info.nz_used;
816: #endif
817: PetscCall(PetscInfo(pc, "%s: %" PetscInt_FMT ") N=%" PetscInt_FMT ", n data cols=%" PetscInt_FMT ", nnz/row (ave)=%" PetscInt_FMT ", %d active pes\n", ((PetscObject)pc)->prefix, level1, M, pc_gamg->data_cell_cols, (PetscInt)(info.nz_used / (PetscReal)M), nactivepe));
819: #if defined(GAMG_STAGES)
820: PetscCall(PetscLogStagePop());
821: #endif
822: /* stop if one node or one proc -- could pull back for singular problems */
823: if ((pc_gamg->data_cell_cols && M / pc_gamg->data_cell_cols < 2) || (!pc_gamg->data_cell_cols && M / bs < 2)) {
824: 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));
825: level++;
826: break;
827: } else if (level == PETSC_MG_MAXLEVELS - 2) { /* stop if we are limited by PC_MG_MAXLEVELS */
828: PetscCall(PetscInfo(pc, "%s: HARD stop of coarsening on level %" PetscInt_FMT ". PC_MG_MAXLEVELS reached\n", ((PetscObject)pc)->prefix, level));
829: level++;
830: break;
831: }
832: } /* levels */
833: PetscCall(PetscFree(pc_gamg->data));
835: PetscCall(PetscInfo(pc, "%s: %" PetscInt_FMT " levels, operator complexity = %g\n", ((PetscObject)pc)->prefix, level + 1, nnztot / nnz0));
836: pc_gamg->Nlevels = level + 1;
837: fine_level = level;
838: PetscCall(PCMGSetLevels(pc, pc_gamg->Nlevels, NULL));
840: if (pc_gamg->Nlevels > 1) { /* don't setup MG if one level */
842: /* set default smoothers & set operators */
843: for (lidx = 1, level = pc_gamg->Nlevels - 2; lidx <= fine_level; lidx++, level--) {
844: KSP smoother;
845: PC subpc;
847: PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
848: PetscCall(KSPGetPC(smoother, &subpc));
850: PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE));
851: /* set ops */
852: PetscCall(KSPSetOperators(smoother, Aarr[level], Aarr[level]));
853: PetscCall(PCMGSetInterpolation(pc, lidx, Parr[level + 1]));
855: /* set defaults */
856: PetscCall(KSPSetType(smoother, KSPCHEBYSHEV));
858: /* set blocks for ASM smoother that uses the 'aggregates' */
859: if (pc_gamg->use_aggs_in_asm || pc_gamg->asm_hem_aggs) {
860: PetscInt sz;
861: IS *iss;
863: sz = nASMBlocksArr[level];
864: iss = ASMLocalIDsArr[level];
865: PetscCall(PCSetType(subpc, PCASM));
866: PetscCall(PCASMSetOverlap(subpc, 0));
867: PetscCall(PCASMSetType(subpc, PC_ASM_BASIC));
868: if (!sz) {
869: IS is;
870: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_COPY_VALUES, &is));
871: PetscCall(PCASMSetLocalSubdomains(subpc, 1, NULL, &is));
872: PetscCall(ISDestroy(&is));
873: } else {
874: PetscInt kk;
875: PetscCall(PCASMSetLocalSubdomains(subpc, sz, iss, NULL));
876: for (kk = 0; kk < sz; kk++) PetscCall(ISDestroy(&iss[kk]));
877: PetscCall(PetscFree(iss));
878: }
879: ASMLocalIDsArr[level] = NULL;
880: nASMBlocksArr[level] = 0;
881: } else {
882: PetscCall(PCSetType(subpc, PCJACOBI));
883: }
884: }
885: {
886: /* coarse grid */
887: KSP smoother, *k2;
888: PC subpc, pc2;
889: PetscInt ii, first;
890: Mat Lmat = Aarr[pc_gamg->Nlevels - 1];
891: lidx = 0;
893: PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
894: PetscCall(KSPSetOperators(smoother, Lmat, Lmat));
895: if (!pc_gamg->use_parallel_coarse_grid_solver) {
896: PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE));
897: PetscCall(KSPGetPC(smoother, &subpc));
898: PetscCall(PCSetType(subpc, PCBJACOBI));
899: PetscCall(PCSetUp(subpc));
900: PetscCall(PCBJacobiGetSubKSP(subpc, &ii, &first, &k2));
901: PetscCheck(ii == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " is not one", ii);
902: PetscCall(KSPGetPC(k2[0], &pc2));
903: PetscCall(PCSetType(pc2, PCLU));
904: PetscCall(PCFactorSetShiftType(pc2, MAT_SHIFT_INBLOCKS));
905: PetscCall(KSPSetTolerances(k2[0], PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 1));
906: PetscCall(KSPSetType(k2[0], KSPPREONLY));
907: }
908: }
910: /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */
911: PetscObjectOptionsBegin((PetscObject)pc);
912: PetscCall(PCSetFromOptions_MG(pc, PetscOptionsObject));
913: PetscOptionsEnd();
914: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
916: /* set cheby eigen estimates from SA to use in the solver */
917: if (pc_gamg->use_sa_esteig) {
918: for (lidx = 1, level = pc_gamg->Nlevels - 2; level >= 0; lidx++, level--) {
919: KSP smoother;
920: PetscBool ischeb;
922: PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
923: PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb));
924: if (ischeb) {
925: KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data;
927: // The command line will override these settings because KSPSetFromOptions is called in PCSetUp_MG
928: if (mg->max_eigen_DinvA[level] > 0) {
929: // SA uses Jacobi for P; we use SA estimates if the smoother is also Jacobi or if the user explicitly requested it.
930: // TODO: This should test whether it's the same Jacobi variant (DIAG, ROWSUM, etc.)
931: PetscReal emax, emin;
933: emin = mg->min_eigen_DinvA[level];
934: emax = mg->max_eigen_DinvA[level];
935: 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));
936: cheb->emin_provided = emin;
937: cheb->emax_provided = emax;
938: }
939: }
940: }
941: }
943: PetscCall(PCSetUp_MG(pc));
945: /* clean up */
946: for (level = 1; level < pc_gamg->Nlevels; level++) {
947: PetscCall(MatDestroy(&Parr[level]));
948: PetscCall(MatDestroy(&Aarr[level]));
949: }
950: } else {
951: KSP smoother;
953: PetscCall(PetscInfo(pc, "%s: One level solver used (system is seen as DD). Using default solver.\n", ((PetscObject)pc)->prefix));
954: PetscCall(PCMGGetSmoother(pc, 0, &smoother));
955: PetscCall(KSPSetOperators(smoother, Aarr[0], Aarr[0]));
956: PetscCall(KSPSetType(smoother, KSPPREONLY));
957: PetscCall(PCSetUp_MG(pc));
958: }
959: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
960: PetscFunctionReturn(PETSC_SUCCESS);
961: }
963: /*
964: PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner
965: that was created with PCCreate_GAMG().
967: Input Parameter:
968: . pc - the preconditioner context
970: Application Interface Routine: PCDestroy()
971: */
972: PetscErrorCode PCDestroy_GAMG(PC pc)
973: {
974: PC_MG *mg = (PC_MG *)pc->data;
975: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
977: PetscFunctionBegin;
978: PetscCall(PCReset_GAMG(pc));
979: if (pc_gamg->ops->destroy) PetscCall((*pc_gamg->ops->destroy)(pc));
980: PetscCall(PetscFree(pc_gamg->ops));
981: PetscCall(PetscFree(pc_gamg->gamg_type_name));
982: PetscCall(PetscFree(pc_gamg));
983: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", NULL));
984: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", NULL));
985: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", NULL));
986: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", NULL));
987: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", NULL));
988: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", NULL));
989: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", NULL));
990: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", NULL));
991: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", NULL));
992: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", NULL));
993: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", NULL));
994: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", NULL));
995: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", NULL));
996: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", NULL));
997: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", NULL));
998: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", NULL));
999: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", NULL));
1000: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetHEM_C", NULL));
1001: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndices_C", NULL));
1002: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndex_C", NULL));
1003: PetscCall(PCDestroy_MG(pc));
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: /*@
1008: PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction in `PCGAMG`
1010: Logically Collective
1012: Input Parameters:
1013: + pc - the preconditioner context
1014: - n - the number of equations
1016: Options Database Key:
1017: . -pc_gamg_process_eq_limit limit - set the limit
1019: Level: intermediate
1021: Note:
1022: `PCGAMG` will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process
1023: that has degrees of freedom
1025: Developer Note:
1026: Should be named `PCGAMGSetProcessEquationLimit()`.
1028: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()`
1029: @*/
1030: PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n)
1031: {
1032: PetscFunctionBegin;
1034: PetscTryMethod(pc, "PCGAMGSetProcEqLim_C", (PC, PetscInt), (pc, n));
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n)
1039: {
1040: PC_MG *mg = (PC_MG *)pc->data;
1041: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1043: PetscFunctionBegin;
1044: if (n > 0) pc_gamg->min_eq_proc = n;
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: /*@
1049: PCGAMGSetCoarseEqLim - Set maximum number of equations on the coarsest grid of `PCGAMG`
1051: Collective
1053: Input Parameters:
1054: + pc - the preconditioner context
1055: - n - maximum number of equations to aim for
1057: Options Database Key:
1058: . -pc_gamg_coarse_eq_limit limit - set the limit
1060: Level: intermediate
1062: Note:
1063: For example -pc_gamg_coarse_eq_limit 1000 will stop coarsening once the coarse grid
1064: has less than 1000 unknowns.
1066: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()`,
1067: `PCGAMGSetParallelCoarseGridSolve()`
1068: @*/
1069: PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n)
1070: {
1071: PetscFunctionBegin;
1073: PetscTryMethod(pc, "PCGAMGSetCoarseEqLim_C", (PC, PetscInt), (pc, n));
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n)
1078: {
1079: PC_MG *mg = (PC_MG *)pc->data;
1080: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1082: PetscFunctionBegin;
1083: if (n > 0) pc_gamg->coarse_eq_limit = n;
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }
1087: /*@
1088: PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids when reducing the number of MPI processes used
1090: Collective
1092: Input Parameters:
1093: + pc - the preconditioner context
1094: - n - `PETSC_TRUE` or `PETSC_FALSE`
1096: Options Database Key:
1097: . -pc_gamg_repartition (true|false) - turn on the repartitioning
1099: Level: intermediate
1101: Note:
1102: This will generally improve the load balance of the work on each level so the solves will be faster but it increases the
1103: preconditioner setup time.
1105: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`
1106: @*/
1107: PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n)
1108: {
1109: PetscFunctionBegin;
1111: PetscTryMethod(pc, "PCGAMGSetRepartition_C", (PC, PetscBool), (pc, n));
1112: PetscFunctionReturn(PETSC_SUCCESS);
1113: }
1115: static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n)
1116: {
1117: PC_MG *mg = (PC_MG *)pc->data;
1118: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1120: PetscFunctionBegin;
1121: pc_gamg->repart = n;
1122: PetscFunctionReturn(PETSC_SUCCESS);
1123: }
1125: /*@
1126: PCGAMGSetUseSAEstEig - Use the eigenvalue estimate from smoothed aggregation for the Chebyshev smoother during the solution process
1128: Collective
1130: Input Parameters:
1131: + pc - the preconditioner context
1132: - b - flag
1134: Options Database Key:
1135: . -pc_gamg_use_sa_esteig (true|false) - use the eigen estimate
1137: Level: advanced
1139: Notes:
1140: 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$.
1141: Eigenvalue estimates (based on a few `PCCG` or `PCGMRES` iterations) are computed to choose $\omega$ so that this is a stable smoothing operation.
1142: If `KSPCHEBYSHEV` with `PCJACOBI` (diagonal) preconditioning is used for smoothing on the finest level, then the eigenvalue estimates
1143: can be reused during the solution process.
1144: This option is only used when the smoother uses `PCJACOBI`, and should be turned off when a different `PCJacobiType` is used.
1146: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `PCGAMGSetRecomputeEstEig()`
1147: @*/
1148: PetscErrorCode PCGAMGSetUseSAEstEig(PC pc, PetscBool b)
1149: {
1150: PetscFunctionBegin;
1152: PetscTryMethod(pc, "PCGAMGSetUseSAEstEig_C", (PC, PetscBool), (pc, b));
1153: PetscFunctionReturn(PETSC_SUCCESS);
1154: }
1156: static PetscErrorCode PCGAMGSetUseSAEstEig_GAMG(PC pc, PetscBool b)
1157: {
1158: PC_MG *mg = (PC_MG *)pc->data;
1159: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1161: PetscFunctionBegin;
1162: pc_gamg->use_sa_esteig = b;
1163: PetscFunctionReturn(PETSC_SUCCESS);
1164: }
1166: /*@
1167: PCGAMGSetRecomputeEstEig - Set flag for Chebyshev smoothers to recompute the eigenvalue estimates when a new matrix is used
1169: Collective
1171: Input Parameters:
1172: + pc - the preconditioner context
1173: - b - flag, default is `PETSC_TRUE`
1175: Options Database Key:
1176: . -pc_gamg_recompute_esteig (true|false) - recompute the eigenvalue estimate
1178: Level: advanced
1180: Note:
1181: If the matrix changes only slightly in a new solve using `PETSC_FALSE` will save time in setting up the preconditioner
1182: and may not affect the solution time much.
1184: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`
1185: @*/
1186: PetscErrorCode PCGAMGSetRecomputeEstEig(PC pc, PetscBool b)
1187: {
1188: PetscFunctionBegin;
1190: PetscTryMethod(pc, "PCGAMGSetRecomputeEstEig_C", (PC, PetscBool), (pc, b));
1191: PetscFunctionReturn(PETSC_SUCCESS);
1192: }
1194: static PetscErrorCode PCGAMGSetRecomputeEstEig_GAMG(PC pc, PetscBool b)
1195: {
1196: PC_MG *mg = (PC_MG *)pc->data;
1197: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1199: PetscFunctionBegin;
1200: pc_gamg->recompute_esteig = b;
1201: PetscFunctionReturn(PETSC_SUCCESS);
1202: }
1204: /*@
1205: PCGAMGSetEigenvalues - Set WHAT eigenvalues WHY?
1207: Collective
1209: Input Parameters:
1210: + pc - the preconditioner context
1211: . emax - max eigenvalue
1212: - emin - min eigenvalue
1214: Options Database Key:
1215: . -pc_gamg_eigenvalues emin,emax - estimates of the eigenvalues
1217: Level: intermediate
1219: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetUseSAEstEig()`
1220: @*/
1221: PetscErrorCode PCGAMGSetEigenvalues(PC pc, PetscReal emax, PetscReal emin)
1222: {
1223: PetscFunctionBegin;
1225: PetscTryMethod(pc, "PCGAMGSetEigenvalues_C", (PC, PetscReal, PetscReal), (pc, emax, emin));
1226: PetscFunctionReturn(PETSC_SUCCESS);
1227: }
1229: static PetscErrorCode PCGAMGSetEigenvalues_GAMG(PC pc, PetscReal emax, PetscReal emin)
1230: {
1231: PC_MG *mg = (PC_MG *)pc->data;
1232: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1234: PetscFunctionBegin;
1235: 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);
1236: 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);
1237: pc_gamg->emax = emax;
1238: pc_gamg->emin = emin;
1239: PetscFunctionReturn(PETSC_SUCCESS);
1240: }
1242: /*@
1243: PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding a `PCGAMG` algebraic multigrid preconditioner when the matrix has changed
1245: Collective
1247: Input Parameters:
1248: + pc - the preconditioner context
1249: - n - `PETSC_TRUE` or `PETSC_FALSE`
1251: Options Database Key:
1252: . -pc_gamg_reuse_interpolation (true|false) - reuse the previous interpolation
1254: Level: intermediate
1256: Note:
1257: May negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows
1258: rebuilding the preconditioner quicker.
1260: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
1261: @*/
1262: PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n)
1263: {
1264: PetscFunctionBegin;
1266: PetscTryMethod(pc, "PCGAMGSetReuseInterpolation_C", (PC, PetscBool), (pc, n));
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n)
1271: {
1272: PC_MG *mg = (PC_MG *)pc->data;
1273: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1275: PetscFunctionBegin;
1276: pc_gamg->reuse_prol = n;
1277: PetscFunctionReturn(PETSC_SUCCESS);
1278: }
1280: /*@
1281: PCGAMGASMSetUseAggs - Have the `PCGAMG` smoother on each level use `PCASM` where the aggregates defined by the coarsening process are
1282: used as the subdomains for the additive Schwarz preconditioner smoother
1284: Collective
1286: Input Parameters:
1287: + pc - the preconditioner context
1288: - flg - `PETSC_TRUE` to use aggregates, `PETSC_FALSE` to not
1290: Options Database Key:
1291: . -pc_gamg_asm_use_agg (true|false) - use aggregates to define the additive Schwarz subdomains
1293: Level: intermediate
1295: Note:
1296: This option automatically sets the preconditioner on the levels to be `PCASM`.
1298: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCASM`, `PCSetType`
1299: @*/
1300: PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg)
1301: {
1302: PetscFunctionBegin;
1304: PetscTryMethod(pc, "PCGAMGASMSetUseAggs_C", (PC, PetscBool), (pc, flg));
1305: PetscFunctionReturn(PETSC_SUCCESS);
1306: }
1308: static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg)
1309: {
1310: PC_MG *mg = (PC_MG *)pc->data;
1311: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1313: PetscFunctionBegin;
1314: pc_gamg->use_aggs_in_asm = flg;
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: /*@
1319: PCGAMGSetParallelCoarseGridSolve - allow a parallel coarse grid solver
1321: Collective
1323: Input Parameters:
1324: + pc - the preconditioner context
1325: - flg - `PETSC_TRUE` to not force coarse grid onto one processor
1327: Options Database Key:
1328: . -pc_gamg_parallel_coarse_grid_solver (true|false) - use a parallel coarse grid solver
1330: Level: intermediate
1332: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetCpuPinCoarseGrids()`, `PCGAMGSetRankReductionFactors()`
1333: @*/
1334: PetscErrorCode PCGAMGSetParallelCoarseGridSolve(PC pc, PetscBool flg)
1335: {
1336: PetscFunctionBegin;
1338: PetscTryMethod(pc, "PCGAMGSetParallelCoarseGridSolve_C", (PC, PetscBool), (pc, flg));
1339: PetscFunctionReturn(PETSC_SUCCESS);
1340: }
1342: static PetscErrorCode PCGAMGSetParallelCoarseGridSolve_GAMG(PC pc, PetscBool flg)
1343: {
1344: PC_MG *mg = (PC_MG *)pc->data;
1345: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1347: PetscFunctionBegin;
1348: pc_gamg->use_parallel_coarse_grid_solver = flg;
1349: PetscFunctionReturn(PETSC_SUCCESS);
1350: }
1352: /*@
1353: 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
1355: Collective
1357: Input Parameters:
1358: + pc - the preconditioner context
1359: - flg - `PETSC_TRUE` to pin coarse grids to the CPU
1361: Options Database Key:
1362: . -pc_gamg_cpu_pin_coarse_grids (true|false) - pin the coarse grids to the CPU
1364: Level: intermediate
1366: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetParallelCoarseGridSolve()`
1367: @*/
1368: PetscErrorCode PCGAMGSetCpuPinCoarseGrids(PC pc, PetscBool flg)
1369: {
1370: PetscFunctionBegin;
1372: PetscTryMethod(pc, "PCGAMGSetCpuPinCoarseGrids_C", (PC, PetscBool), (pc, flg));
1373: PetscFunctionReturn(PETSC_SUCCESS);
1374: }
1376: static PetscErrorCode PCGAMGSetCpuPinCoarseGrids_GAMG(PC pc, PetscBool flg)
1377: {
1378: PC_MG *mg = (PC_MG *)pc->data;
1379: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1381: PetscFunctionBegin;
1382: pc_gamg->cpu_pin_coarse_grids = flg;
1383: PetscFunctionReturn(PETSC_SUCCESS);
1384: }
1386: /*@
1387: PCGAMGSetCoarseGridLayoutType - place coarse grids on processors with natural order (compact type)
1389: Collective
1391: Input Parameters:
1392: + pc - the preconditioner context
1393: - flg - `PCGAMGLayoutType` type, either `PCGAMG_LAYOUT_COMPACT` or `PCGAMG_LAYOUT_SPREAD`
1395: Options Database Key:
1396: . -pc_gamg_coarse_grid_layout_type (compact|spread) - place the coarse grids with natural ordering
1398: Level: intermediate
1400: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetCpuPinCoarseGrids()`, `PCGAMGLayoutType`, `PCGAMG_LAYOUT_COMPACT`, `PCGAMG_LAYOUT_SPREAD`
1401: @*/
1402: PetscErrorCode PCGAMGSetCoarseGridLayoutType(PC pc, PCGAMGLayoutType flg)
1403: {
1404: PetscFunctionBegin;
1406: PetscTryMethod(pc, "PCGAMGSetCoarseGridLayoutType_C", (PC, PCGAMGLayoutType), (pc, flg));
1407: PetscFunctionReturn(PETSC_SUCCESS);
1408: }
1410: static PetscErrorCode PCGAMGSetCoarseGridLayoutType_GAMG(PC pc, PCGAMGLayoutType flg)
1411: {
1412: PC_MG *mg = (PC_MG *)pc->data;
1413: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1415: PetscFunctionBegin;
1416: pc_gamg->layout_type = flg;
1417: PetscFunctionReturn(PETSC_SUCCESS);
1418: }
1420: /*@
1421: PCGAMGSetNlevels - Sets the maximum number of levels `PCGAMG` will use
1423: Collective
1425: Input Parameters:
1426: + pc - the preconditioner
1427: - n - the maximum number of levels to use
1429: Options Database Key:
1430: . -pc_mg_levels n - set the maximum number of levels to allow
1432: Level: intermediate
1434: Developer Notes:
1435: Should be called `PCGAMGSetMaximumNumberlevels()` and possible be shared with `PCMG`
1437: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
1438: @*/
1439: PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n)
1440: {
1441: PetscFunctionBegin;
1443: PetscTryMethod(pc, "PCGAMGSetNlevels_C", (PC, PetscInt), (pc, n));
1444: PetscFunctionReturn(PETSC_SUCCESS);
1445: }
1447: static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n)
1448: {
1449: PC_MG *mg = (PC_MG *)pc->data;
1450: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1452: PetscFunctionBegin;
1453: pc_gamg->Nlevels = n;
1454: PetscFunctionReturn(PETSC_SUCCESS);
1455: }
1457: /*@
1458: PCGAMGASMSetHEM - Sets the number of HEM matching passed
1460: Collective
1462: Input Parameters:
1463: + pc - the preconditioner
1464: - n - number of HEM matching passed to construct ASM subdomains
1466: Options Database Key:
1467: . -pc_gamg_asm_hem n - set the number of HEM matching passed
1469: Level: intermediate
1471: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
1472: @*/
1473: PetscErrorCode PCGAMGASMSetHEM(PC pc, PetscInt n)
1474: {
1475: PetscFunctionBegin;
1477: PetscTryMethod(pc, "PCGAMGASMSetHEM_C", (PC, PetscInt), (pc, n));
1478: PetscFunctionReturn(PETSC_SUCCESS);
1479: }
1481: static PetscErrorCode PCGAMGASMSetHEM_GAMG(PC pc, PetscInt n)
1482: {
1483: PC_MG *mg = (PC_MG *)pc->data;
1484: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1486: PetscFunctionBegin;
1487: pc_gamg->asm_hem_aggs = n;
1488: PetscFunctionReturn(PETSC_SUCCESS);
1489: }
1491: /*@
1492: PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph
1494: Not Collective
1496: Input Parameters:
1497: + pc - the preconditioner context
1498: . 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
1499: - n - number of threshold values provided in array
1501: Options Database Key:
1502: . -pc_gamg_threshold l0,l1,... - the thresholds to drop edges
1504: Level: intermediate
1506: Notes:
1507: Before coarsening or aggregating the graph, `PCGAMG` removes small values from the graph with this threshold, reducing the coupling
1508: in the graph and yielding a different (perhaps better) coarser set of points.
1510: Increasing the threshold decreases the rate of coarsening. Conversely reducing the threshold increases the rate of coarsening (aggressive coarsening)
1511: and thereby reduces the complexity of the coarse grids, and generally results in slower solver convergence rates.
1513: If `n` is less than the total number of coarsenings (see `PCGAMGSetNlevels()`), then threshold scaling (see `PCGAMGSetThresholdScale()`) is used for each successive coarsening.
1514: In this case, `PCGAMGSetThresholdScale()` must be called before `PCGAMGSetThreshold()`.
1516: If `n` is greater than the total number of levels, the excess entries in threshold are ignored.
1518: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`,
1519: `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetThresholdScale()`
1520: @*/
1521: PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal v[], PetscInt n)
1522: {
1523: PetscFunctionBegin;
1525: if (n) PetscAssertPointer(v, 2);
1526: PetscTryMethod(pc, "PCGAMGSetThreshold_C", (PC, PetscReal[], PetscInt), (pc, v, n));
1527: PetscFunctionReturn(PETSC_SUCCESS);
1528: }
1530: static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal v[], PetscInt n)
1531: {
1532: PC_MG *mg = (PC_MG *)pc->data;
1533: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1534: PetscInt i;
1536: PetscFunctionBegin;
1537: for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->threshold[i] = v[i];
1538: for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale;
1539: PetscFunctionReturn(PETSC_SUCCESS);
1540: }
1542: /*@
1543: PCGAMGSetRankReductionFactors - Set a manual schedule for MPI process reduction on coarse grids
1545: Collective
1547: Input Parameters:
1548: + pc - the preconditioner context
1549: . v - array of reduction factors. 0 for first value forces a reduction to one process/device on first level in CUDA
1550: - n - number of values provided in array
1552: Options Database Key:
1553: . -pc_gamg_rank_reduction_factors r0,r1,... - provide the schedule
1555: Level: intermediate
1557: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetParallelCoarseGridSolve()`
1558: @*/
1559: PetscErrorCode PCGAMGSetRankReductionFactors(PC pc, PetscInt v[], PetscInt n)
1560: {
1561: PetscFunctionBegin;
1563: if (n) PetscAssertPointer(v, 2);
1564: PetscTryMethod(pc, "PCGAMGSetRankReductionFactors_C", (PC, PetscInt[], PetscInt), (pc, v, n));
1565: PetscFunctionReturn(PETSC_SUCCESS);
1566: }
1568: static PetscErrorCode PCGAMGSetRankReductionFactors_GAMG(PC pc, PetscInt v[], PetscInt n)
1569: {
1570: PC_MG *mg = (PC_MG *)pc->data;
1571: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1572: PetscInt i;
1574: PetscFunctionBegin;
1575: for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->level_reduction_factors[i] = v[i];
1576: for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->level_reduction_factors[i] = -1; /* 0 stop putting one process/device on first level */
1577: PetscFunctionReturn(PETSC_SUCCESS);
1578: }
1580: /*@
1581: PCGAMGSetThresholdScale - Relative threshold reduction at each level
1583: Not Collective
1585: Input Parameters:
1586: + pc - the preconditioner context
1587: - v - the threshold value reduction, usually < 1.0
1589: Options Database Key:
1590: . -pc_gamg_threshold_scale v - set the relative threshold reduction on each level
1592: Level: advanced
1594: Note:
1595: The initial threshold (for an arbitrary number of levels starting from the finest) can be set with `PCGAMGSetThreshold()`.
1596: This scaling is used for each subsequent coarsening, but must be called before `PCGAMGSetThreshold()`.
1598: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`
1599: @*/
1600: PetscErrorCode PCGAMGSetThresholdScale(PC pc, PetscReal v)
1601: {
1602: PetscFunctionBegin;
1604: PetscTryMethod(pc, "PCGAMGSetThresholdScale_C", (PC, PetscReal), (pc, v));
1605: PetscFunctionReturn(PETSC_SUCCESS);
1606: }
1608: static PetscErrorCode PCGAMGSetThresholdScale_GAMG(PC pc, PetscReal v)
1609: {
1610: PC_MG *mg = (PC_MG *)pc->data;
1611: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1613: PetscFunctionBegin;
1614: pc_gamg->threshold_scale = v;
1615: PetscFunctionReturn(PETSC_SUCCESS);
1616: }
1618: /*@
1619: PCGAMGSetType - Set the type of algorithm `PCGAMG` should use
1621: Collective
1623: Input Parameters:
1624: + pc - the preconditioner context
1625: - type - `PCGAMGAGG`, `PCGAMGGEO`, or `PCGAMGCLASSICAL`
1627: Options Database Key:
1628: . -pc_gamg_type (agg|geo|classical) - type of algebraic multigrid to apply - only agg is supported
1630: Level: intermediate
1632: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMGGetType()`, `PCGAMG`, `PCGAMGType`
1633: @*/
1634: PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type)
1635: {
1636: PetscFunctionBegin;
1638: PetscTryMethod(pc, "PCGAMGSetType_C", (PC, PCGAMGType), (pc, type));
1639: PetscFunctionReturn(PETSC_SUCCESS);
1640: }
1642: /*@
1643: PCGAMGGetType - Get the type of algorithm `PCGAMG` will use
1645: Collective
1647: Input Parameter:
1648: . pc - the preconditioner context
1650: Output Parameter:
1651: . type - the type of algorithm used
1653: Level: intermediate
1655: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetType()`, `PCGAMGType`
1656: @*/
1657: PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type)
1658: {
1659: PetscFunctionBegin;
1661: PetscUseMethod(pc, "PCGAMGGetType_C", (PC, PCGAMGType *), (pc, type));
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type)
1666: {
1667: PC_MG *mg = (PC_MG *)pc->data;
1668: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1670: PetscFunctionBegin;
1671: *type = pc_gamg->type;
1672: PetscFunctionReturn(PETSC_SUCCESS);
1673: }
1675: static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type)
1676: {
1677: PC_MG *mg = (PC_MG *)pc->data;
1678: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1679: PetscErrorCode (*r)(PC);
1681: PetscFunctionBegin;
1682: pc_gamg->type = type;
1683: PetscCall(PetscFunctionListFind(GAMGList, type, &r));
1684: PetscCheck(r, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown GAMG type %s given", type);
1685: if (pc_gamg->ops->destroy) {
1686: PetscCall((*pc_gamg->ops->destroy)(pc));
1687: PetscCall(PetscMemzero(pc_gamg->ops, sizeof(struct _PCGAMGOps)));
1688: pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1689: /* cleaning up common data in pc_gamg - this should disappear someday */
1690: pc_gamg->data_cell_cols = 0;
1691: pc_gamg->data_cell_rows = 0;
1692: pc_gamg->orig_data_cell_cols = 0;
1693: pc_gamg->orig_data_cell_rows = 0;
1694: PetscCall(PetscFree(pc_gamg->data));
1695: pc_gamg->data_sz = 0;
1696: }
1697: PetscCall(PetscFree(pc_gamg->gamg_type_name));
1698: PetscCall(PetscStrallocpy(type, &pc_gamg->gamg_type_name));
1699: PetscCall((*r)(pc));
1700: PetscFunctionReturn(PETSC_SUCCESS);
1701: }
1703: static PetscErrorCode PCView_GAMG(PC pc, PetscViewer viewer)
1704: {
1705: PC_MG *mg = (PC_MG *)pc->data;
1706: PC_MG_Levels **mglevels = mg->levels;
1707: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1708: PetscReal gc, oc;
1710: PetscFunctionBegin;
1711: PetscCall(PetscViewerASCIIPrintf(viewer, " GAMG specific options\n"));
1712: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for dropping small values in graph on each level ="));
1713: for (PetscInt i = 0; i < mg->nlevels; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)pc_gamg->threshold[i]));
1714: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1715: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold scaling factor for each level not specified = %g\n", (double)pc_gamg->threshold_scale));
1716: if (pc_gamg->use_aggs_in_asm) PetscCall(PetscViewerASCIIPrintf(viewer, " Using aggregates from GAMG coarsening process to define subdomains for PCASM\n")); // this take precedence
1717: else if (pc_gamg->asm_hem_aggs) {
1718: PetscCall(PetscViewerASCIIPrintf(viewer, " Using aggregates made with %" PetscInt_FMT " applications of heavy edge matching (HEM) to define subdomains for PCASM\n", pc_gamg->asm_hem_aggs));
1719: PetscCall(PetscViewerASCIIPushTab(viewer));
1720: PetscCall(PetscViewerASCIIPushTab(viewer));
1721: PetscCall(PetscViewerASCIIPushTab(viewer));
1722: PetscCall(PetscViewerASCIIPushTab(viewer));
1723: PetscCall(MatCoarsenView(pc_gamg->asm_crs, viewer));
1724: PetscCall(PetscViewerASCIIPopTab(viewer));
1725: PetscCall(PetscViewerASCIIPopTab(viewer));
1726: PetscCall(PetscViewerASCIIPopTab(viewer));
1727: }
1728: 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"));
1729: if (pc_gamg->injection_index_size) {
1730: PetscCall(PetscViewerASCIIPrintf(viewer, " Using injection restriction/prolongation on first level, dofs:"));
1731: for (int i = 0; i < pc_gamg->injection_index_size; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, pc_gamg->injection_index[i]));
1732: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1733: }
1734: if (pc_gamg->ops->view) PetscCall((*pc_gamg->ops->view)(pc, viewer));
1735: gc = oc = 0;
1736: PetscCall(PCMGGetGridComplexity(pc, &gc, &oc));
1737: PetscCall(PetscViewerASCIIPrintf(viewer, " Complexity: grid = %g operator = %g\n", (double)gc, (double)oc));
1738: PetscCall(PetscViewerASCIIPrintf(viewer, " Per-level complexity: op = operator, int = interpolation\n"));
1739: PetscCall(PetscViewerASCIIPrintf(viewer, " #equations | #active PEs | avg nnz/row op | avg nnz/row int\n"));
1740: for (PetscInt i = 0; i < mg->nlevels; i++) {
1741: MatInfo info;
1742: Mat A;
1743: PetscReal rd[3];
1744: PetscInt rst, ren, N;
1746: PetscCall(KSPGetOperators(mglevels[i]->smoothd, NULL, &A));
1747: PetscCall(MatGetOwnershipRange(A, &rst, &ren));
1748: PetscCall(MatGetSize(A, &N, NULL));
1749: PetscCall(MatGetInfo(A, MAT_LOCAL, &info));
1750: rd[0] = (ren - rst > 0) ? 1 : 0;
1751: rd[1] = info.nz_used;
1752: rd[2] = 0;
1753: if (i) {
1754: Mat P;
1755: PetscCall(PCMGGetInterpolation(pc, i, &P));
1756: PetscCall(MatGetInfo(P, MAT_LOCAL, &info));
1757: rd[2] = info.nz_used;
1758: }
1759: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, rd, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)pc)));
1760: PetscCall(PetscViewerASCIIPrintf(viewer, " %12" PetscInt_FMT " %12" PetscInt_FMT " %12" PetscInt_FMT " %12" PetscInt_FMT "\n", N, (PetscInt)rd[0], (PetscInt)PetscCeilReal(rd[1] / N), (PetscInt)PetscCeilReal(rd[2] / N)));
1761: }
1762: PetscFunctionReturn(PETSC_SUCCESS);
1763: }
1765: /*@
1766: PCGAMGSetInjectionIndex - Array of subset of variables per vertex to inject into coarse grid space
1768: Logically Collective
1770: Input Parameters:
1771: + pc - the coarsen context
1772: . n - number of indices
1773: - idx - array of indices
1775: Options Database Key:
1776: . -pc_gamg_injection_index - array of subset of variables per vertex to use for injection coarse grid space
1778: Level: intermediate
1780: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`
1781: @*/
1782: PetscErrorCode PCGAMGSetInjectionIndex(PC pc, PetscInt n, PetscInt idx[])
1783: {
1784: PetscFunctionBegin;
1787: PetscTryMethod(pc, "PCGAMGSetInjectionIndex_C", (PC, PetscInt, PetscInt[]), (pc, n, idx));
1788: PetscFunctionReturn(PETSC_SUCCESS);
1789: }
1791: static PetscErrorCode PCGAMGSetInjectionIndex_GAMG(PC pc, PetscInt n, PetscInt idx[])
1792: {
1793: PC_MG *mg = (PC_MG *)pc->data;
1794: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1796: PetscFunctionBegin;
1797: pc_gamg->injection_index_size = n;
1798: PetscCheck(n < MAT_COARSEN_STRENGTH_INDEX_SIZE, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "array size %" PetscInt_FMT " larger than max %d", n, MAT_COARSEN_STRENGTH_INDEX_SIZE);
1799: for (PetscInt i = 0; i < n; i++) pc_gamg->injection_index[i] = idx[i];
1800: PetscFunctionReturn(PETSC_SUCCESS);
1801: }
1803: static PetscErrorCode PCSetFromOptions_GAMG(PC pc, PetscOptionItems PetscOptionsObject)
1804: {
1805: PC_MG *mg = (PC_MG *)pc->data;
1806: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1807: PetscBool flag;
1808: MPI_Comm comm;
1809: char prefix[256], tname[32];
1810: PetscInt i, n;
1811: const char *pcpre;
1812: static const char *LayoutTypes[] = {"compact", "spread", "PCGAMGLayoutType", "PC_GAMG_LAYOUT", NULL};
1814: PetscFunctionBegin;
1815: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1816: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG options");
1817: PetscCall(PetscOptionsInt("-pc_mg_levels", "Number of Levels", "PCMGSetLevels", -1, &n, &flag));
1818: PetscCheck(!flag, comm, PETSC_ERR_ARG_WRONG, "Invalid flag -pc_mg_levels. GAMG does not allow the number of levels to be set.");
1819: PetscCall(PetscOptionsFList("-pc_gamg_type", "Type of AMG method (only 'agg' supported and useful)", "PCGAMGSetType", GAMGList, pc_gamg->gamg_type_name, tname, sizeof(tname), &flag));
1820: if (flag) PetscCall(PCGAMGSetType(pc, tname));
1821: PetscCall(PetscOptionsBool("-pc_gamg_repartition", "Repartion coarse grids", "PCGAMGSetRepartition", pc_gamg->repart, &pc_gamg->repart, NULL));
1822: 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));
1823: PetscCall(PetscOptionsBool("-pc_gamg_recompute_esteig", "Set flag to recompute eigen estimates for Chebyshev when matrix changes", "PCGAMGSetRecomputeEstEig", pc_gamg->recompute_esteig, &pc_gamg->recompute_esteig, NULL));
1824: PetscCall(PetscOptionsBool("-pc_gamg_reuse_interpolation", "Reuse prolongation operator", "PCGAMGReuseInterpolation", pc_gamg->reuse_prol, &pc_gamg->reuse_prol, NULL));
1825: 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));
1826: PetscCall(PetscOptionsBool("-pc_gamg_parallel_coarse_grid_solver", "Use parallel coarse grid solver (otherwise put last grid on one process)", "PCGAMGSetParallelCoarseGridSolve", pc_gamg->use_parallel_coarse_grid_solver, &pc_gamg->use_parallel_coarse_grid_solver, NULL));
1827: 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));
1828: 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,
1829: (PetscEnum)pc_gamg->layout_type, (PetscEnum *)&pc_gamg->layout_type, NULL));
1830: 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));
1831: 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));
1832: PetscCall(PetscOptionsInt("-pc_gamg_asm_hem_aggs", "Number of HEM matching passed in aggregates for ASM smoother", "PCGAMGASMSetHEM", pc_gamg->asm_hem_aggs, &pc_gamg->asm_hem_aggs, NULL));
1833: PetscCall(PetscOptionsReal("-pc_gamg_threshold_scale", "Scaling of threshold for each level not specified", "PCGAMGSetThresholdScale", pc_gamg->threshold_scale, &pc_gamg->threshold_scale, NULL));
1834: n = PETSC_MG_MAXLEVELS;
1835: PetscCall(PetscOptionsRealArray("-pc_gamg_threshold", "Relative threshold to use for dropping edges in aggregation graph", "PCGAMGSetThreshold", pc_gamg->threshold, &n, &flag));
1836: if (!flag || n < PETSC_MG_MAXLEVELS) {
1837: if (!flag) n = 1;
1838: i = n;
1839: do {
1840: pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale;
1841: } while (++i < PETSC_MG_MAXLEVELS);
1842: }
1843: PetscCall(PetscOptionsInt("-pc_mg_levels", "Set number of MG levels (should get from base class)", "PCGAMGSetNlevels", pc_gamg->Nlevels, &pc_gamg->Nlevels, NULL));
1844: PetscCheck(pc_gamg->Nlevels <= PETSC_MG_MAXLEVELS, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_mg_levels (%" PetscInt_FMT ") >= PETSC_MG_MAXLEVELS (%d)", pc_gamg->Nlevels, PETSC_MG_MAXLEVELS);
1845: n = PETSC_MG_MAXLEVELS;
1846: 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));
1847: if (!flag) i = 0;
1848: else i = n;
1849: do {
1850: pc_gamg->level_reduction_factors[i] = -1;
1851: } while (++i < PETSC_MG_MAXLEVELS);
1852: {
1853: PetscReal eminmax[2] = {0., 0.};
1854: n = 2;
1855: PetscCall(PetscOptionsRealArray("-pc_gamg_eigenvalues", "extreme eigenvalues for smoothed aggregation", "PCGAMGSetEigenvalues", eminmax, &n, &flag));
1856: if (flag) {
1857: PetscCheck(n == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_gamg_eigenvalues: must specify 2 parameters, min and max eigenvalues");
1858: PetscCall(PCGAMGSetEigenvalues(pc, eminmax[1], eminmax[0]));
1859: }
1860: }
1861: pc_gamg->injection_index_size = MAT_COARSEN_STRENGTH_INDEX_SIZE;
1862: PetscCall(PetscOptionsIntArray("-pc_gamg_injection_index", "Array of indices to use to use injection coarse grid space", "PCGAMGSetInjectionIndex", pc_gamg->injection_index, &pc_gamg->injection_index_size, NULL));
1863: /* set options for subtype */
1864: PetscCall((*pc_gamg->ops->setfromoptions)(pc, PetscOptionsObject));
1866: PetscCall(PCGetOptionsPrefix(pc, &pcpre));
1867: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : ""));
1868: PetscOptionsHeadEnd();
1869: PetscFunctionReturn(PETSC_SUCCESS);
1870: }
1872: /*MC
1873: PCGAMG - Geometric algebraic multigrid (AMG) preconditioner
1875: Options Database Keys:
1876: + -pc_gamg_type type (agg|geo|classical) - see `PCGAMGType`
1877: . -pc_gamg_repartition (true|false) - repartition the degrees of freedom across the coarse grids as they are determined
1878: . -pc_gamg_asm_use_agg (true|false) - use the aggregates from the coarsening process to define the subdomains on each level for the `PCASM` smoother,
1879: this also forces using `-mg_levels_pc_type asm`
1880: . -pc_gamg_process_eq_limit limit - `PCGAMG` will reduce the number of MPI processes used directly on the coarse grids so that there are around limit
1881: equations on each process that has degrees of freedom
1882: . -pc_gamg_coarse_eq_limit limit - set maximum number of equations on coarsest grid to aim for
1883: . -pc_gamg_reuse_interpolation (true|false) - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations (should always be true)
1884: . -pc_gamg_threshold l0,l1,... - before aggregating the graph `PCGAMG` will remove small values from the graph on each level (< 0 does no filtering)
1885: - -pc_gamg_threshold_scale scale - scaling of threshold on each coarser grid if not specified
1887: Options Database Keys for Aggregation:
1888: + -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use with smooth aggregation to construct prolongation
1889: . -pc_gamg_aggressive_coarsening n - number of aggressive coarsening (MIS-2 or square graph) levels from finest.
1890: . -pc_gamg_aggressive_square_graph (true|false) - use square graph ($A^T A$) for coarsening. Otherwise, MIS-k (k=2) is used, see `PCGAMGMISkSetAggressive()`
1891: . -pc_gamg_square[_i]_mat_product_algorithm algorithm - `MatProductAlgorithm` to use when squaring the matrix for aggressive coarsening (on level i < `n`)
1892: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use minimum degree ordering in greedy MIS algorithm
1893: . -pc_gamg_asm_hem_aggs n - number of HEM aggregation steps for `PCASM` smoother
1894: - -pc_gamg_aggressive_mis_k n - number (k) distance in MIS coarsening (>2 is 'aggressive')
1896: Options Database Keys for Multigrid:
1897: + -pc_mg_cycle_type (v|w) - see `PCMGSetCycleType()`
1898: . -pc_mg_distinct_smoothup - configure the up and down (pre and post) smoothers separately, see `PCMGSetDistinctSmoothUp()`
1899: . -pc_mg_type (additive|multiplicative|full|kaskade) - see `PCMGType`
1900: - -pc_mg_levels levels - number of levels of multigrid to use; `PCGAMG` has a heuristic to determine the number of levels so
1901: this is not usually used with `PCGAMG`
1903: Level: intermediate
1905: Notes:
1906: To obtain good performance for `PCGAMG` for vector valued problems you must
1907: call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point
1908: call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator
1910: The many options for `PCMG` also work directly for `PCGAMG` such as controlling the smoothers on each level etc.
1912: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
1913: `MatSetBlockSize()`,
1914: `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`, `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`,
1915: `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`, `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`,
1916: `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`, `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
1917: M*/
1918: PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc)
1919: {
1920: PC_GAMG *pc_gamg;
1921: PC_MG *mg;
1923: PetscFunctionBegin;
1924: /* register AMG type */
1925: PetscCall(PCGAMGInitializePackage());
1927: /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */
1928: PetscCall(PCSetType(pc, PCMG));
1929: PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCGAMG));
1931: /* create a supporting struct and attach it to pc */
1932: PetscCall(PetscNew(&pc_gamg));
1933: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
1934: mg = (PC_MG *)pc->data;
1935: mg->innerctx = pc_gamg;
1937: PetscCall(PetscNew(&pc_gamg->ops));
1939: /* these should be in subctx but repartitioning needs simple arrays */
1940: pc_gamg->data_sz = 0;
1941: pc_gamg->data = NULL;
1943: /* overwrite the pointers of PCMG by the functions of base class PCGAMG */
1944: pc->ops->setfromoptions = PCSetFromOptions_GAMG;
1945: pc->ops->setup = PCSetUp_GAMG;
1946: pc->ops->reset = PCReset_GAMG;
1947: pc->ops->destroy = PCDestroy_GAMG;
1948: mg->view = PCView_GAMG;
1950: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", PCGAMGSetProcEqLim_GAMG));
1951: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", PCGAMGSetCoarseEqLim_GAMG));
1952: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", PCGAMGSetRepartition_GAMG));
1953: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", PCGAMGSetEigenvalues_GAMG));
1954: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", PCGAMGSetUseSAEstEig_GAMG));
1955: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", PCGAMGSetRecomputeEstEig_GAMG));
1956: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", PCGAMGSetReuseInterpolation_GAMG));
1957: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", PCGAMGASMSetUseAggs_GAMG));
1958: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", PCGAMGSetParallelCoarseGridSolve_GAMG));
1959: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", PCGAMGSetCpuPinCoarseGrids_GAMG));
1960: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", PCGAMGSetCoarseGridLayoutType_GAMG));
1961: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", PCGAMGSetThreshold_GAMG));
1962: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", PCGAMGSetRankReductionFactors_GAMG));
1963: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", PCGAMGSetThresholdScale_GAMG));
1964: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", PCGAMGSetType_GAMG));
1965: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", PCGAMGGetType_GAMG));
1966: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", PCGAMGSetNlevels_GAMG));
1967: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetHEM_C", PCGAMGASMSetHEM_GAMG));
1968: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndex_C", PCGAMGSetInjectionIndex_GAMG));
1969: pc_gamg->repart = PETSC_FALSE;
1970: pc_gamg->reuse_prol = PETSC_TRUE;
1971: pc_gamg->use_aggs_in_asm = PETSC_FALSE;
1972: pc_gamg->use_parallel_coarse_grid_solver = PETSC_FALSE;
1973: pc_gamg->cpu_pin_coarse_grids = PETSC_FALSE;
1974: pc_gamg->layout_type = PCGAMG_LAYOUT_SPREAD;
1975: pc_gamg->min_eq_proc = 50;
1976: pc_gamg->asm_hem_aggs = 0;
1977: pc_gamg->coarse_eq_limit = 50;
1978: for (int i = 0; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = -1;
1979: pc_gamg->threshold_scale = 1.;
1980: pc_gamg->Nlevels = PETSC_MG_MAXLEVELS;
1981: pc_gamg->current_level = 0; /* don't need to init really */
1982: pc_gamg->use_sa_esteig = PETSC_TRUE;
1983: pc_gamg->recompute_esteig = PETSC_TRUE;
1984: pc_gamg->emin = 0;
1985: pc_gamg->emax = 0;
1987: pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1989: /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */
1990: PetscCall(PCGAMGSetType(pc, PCGAMGAGG));
1991: PetscFunctionReturn(PETSC_SUCCESS);
1992: }
1994: /*@C
1995: PCGAMGInitializePackage - This function initializes everything in the `PCGAMG` package. It is called
1996: from `PCInitializePackage()`.
1998: Level: developer
2000: .seealso: [](ch_ksp), `PetscInitialize()`
2001: @*/
2002: PetscErrorCode PCGAMGInitializePackage(void)
2003: {
2004: PetscInt l;
2006: PetscFunctionBegin;
2007: if (PCGAMGPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
2008: PCGAMGPackageInitialized = PETSC_TRUE;
2009: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGGEO, PCCreateGAMG_GEO));
2010: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGAGG, PCCreateGAMG_AGG));
2011: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGCLASSICAL, PCCreateGAMG_Classical));
2012: PetscCall(PetscRegisterFinalize(PCGAMGFinalizePackage));
2014: /* general events */
2015: PetscCall(PetscLogEventRegister("PCSetUp_GAMG+", PC_CLASSID, &petsc_gamg_setup_events[GAMG_SETUP]));
2016: PetscCall(PetscLogEventRegister(" PCGAMGCreateG", PC_CLASSID, &petsc_gamg_setup_events[GAMG_GRAPH]));
2017: PetscCall(PetscLogEventRegister(" GAMG Coarsen", PC_CLASSID, &petsc_gamg_setup_events[GAMG_COARSEN]));
2018: PetscCall(PetscLogEventRegister(" GAMG MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[GAMG_MIS]));
2019: PetscCall(PetscLogEventRegister(" PCGAMGProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROL]));
2020: PetscCall(PetscLogEventRegister(" GAMG Prol-col", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLA]));
2021: PetscCall(PetscLogEventRegister(" GAMG Prol-lift", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLB]));
2022: PetscCall(PetscLogEventRegister(" PCGAMGOptProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPT]));
2023: PetscCall(PetscLogEventRegister(" GAMG smooth", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPTSM]));
2024: PetscCall(PetscLogEventRegister(" PCGAMGCreateL", PC_CLASSID, &petsc_gamg_setup_events[GAMG_LEVEL]));
2025: PetscCall(PetscLogEventRegister(" GAMG PtAP", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PTAP]));
2026: PetscCall(PetscLogEventRegister(" GAMG Reduce", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REDUCE]));
2027: PetscCall(PetscLogEventRegister(" GAMG Repart", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REPART]));
2028: /* PetscCall(PetscLogEventRegister(" GAMG Inv-Srt", PC_CLASSID, &petsc_gamg_setup_events[SET13])); */
2029: /* PetscCall(PetscLogEventRegister(" GAMG Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14])); */
2030: /* PetscCall(PetscLogEventRegister(" GAMG Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15])); */
2031: for (l = 0; l < PETSC_MG_MAXLEVELS; l++) {
2032: char ename[32];
2034: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Squ l%02" PetscInt_FMT, l));
2035: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][0]));
2036: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Gal l%02" PetscInt_FMT, l));
2037: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][1]));
2038: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Opt l%02" PetscInt_FMT, l));
2039: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][2]));
2040: }
2041: #if defined(GAMG_STAGES)
2042: { /* create timer stages */
2043: char str[32];
2044: PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", 0));
2045: PetscCall(PetscLogStageRegister(str, &gamg_stages[0]));
2046: }
2047: #endif
2048: PetscFunctionReturn(PETSC_SUCCESS);
2049: }
2051: /*@C
2052: PCGAMGFinalizePackage - This function frees everything from the `PCGAMG` package. It is
2053: called from `PetscFinalize()` automatically.
2055: Level: developer
2057: .seealso: [](ch_ksp), `PetscFinalize()`
2058: @*/
2059: PetscErrorCode PCGAMGFinalizePackage(void)
2060: {
2061: PetscFunctionBegin;
2062: PCGAMGPackageInitialized = PETSC_FALSE;
2063: PetscCall(PetscFunctionListDestroy(&GAMGList));
2064: PetscFunctionReturn(PETSC_SUCCESS);
2065: }
2067: /*@C
2068: PCGAMGRegister - Register a `PCGAMG` implementation.
2070: Input Parameters:
2071: + type - string that will be used as the name of the `PCGAMG` type.
2072: - create - function for creating the gamg context.
2074: Level: developer
2076: .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()`
2077: @*/
2078: PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC))
2079: {
2080: PetscFunctionBegin;
2081: PetscCall(PCGAMGInitializePackage());
2082: PetscCall(PetscFunctionListAdd(&GAMGList, type, create));
2083: PetscFunctionReturn(PETSC_SUCCESS);
2084: }
2086: /*@
2087: PCGAMGCreateGraph - Creates a graph that is used by the ``PCGAMGType`` in the coarsening process
2089: Input Parameters:
2090: + pc - the `PCGAMG`
2091: - A - the matrix, for any level
2093: Output Parameter:
2094: . G - the graph
2096: Level: advanced
2098: .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()`
2099: @*/
2100: PetscErrorCode PCGAMGCreateGraph(PC pc, Mat A, Mat *G)
2101: {
2102: PC_MG *mg = (PC_MG *)pc->data;
2103: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
2105: PetscFunctionBegin;
2106: PetscCall(pc_gamg->ops->creategraph(pc, A, G));
2107: PetscFunctionReturn(PETSC_SUCCESS);
2108: }