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: PetscCall(PCASMSetLocalSubdomains(subpc, sz, iss, NULL));
875: for (PetscInt kk = 0; kk < sz; kk++) PetscCall(ISDestroy(&iss[kk]));
876: PetscCall(PetscFree(iss));
877: }
878: ASMLocalIDsArr[level] = NULL;
879: nASMBlocksArr[level] = 0;
880: } else {
881: PetscCall(PCSetType(subpc, PCJACOBI));
882: }
883: }
884: {
885: /* coarse grid */
886: KSP smoother, *k2;
887: PC subpc, pc2;
888: PetscInt ii, first;
889: Mat Lmat = Aarr[pc_gamg->Nlevels - 1];
890: lidx = 0;
892: PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
893: PetscCall(KSPSetOperators(smoother, Lmat, Lmat));
894: if (!pc_gamg->use_parallel_coarse_grid_solver) {
895: PetscCall(KSPSetNormType(smoother, KSP_NORM_NONE));
896: PetscCall(KSPGetPC(smoother, &subpc));
897: PetscCall(PCSetType(subpc, PCBJACOBI));
898: PetscCall(PCSetUp(subpc));
899: PetscCall(PCBJacobiGetSubKSP(subpc, &ii, &first, &k2));
900: PetscCheck(ii == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " is not one", ii);
901: PetscCall(KSPGetPC(k2[0], &pc2));
902: PetscCall(PCSetType(pc2, PCLU));
903: PetscCall(PCFactorSetShiftType(pc2, MAT_SHIFT_INBLOCKS));
904: PetscCall(KSPSetTolerances(k2[0], PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 1));
905: PetscCall(KSPSetType(k2[0], KSPPREONLY));
906: }
907: }
909: /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */
910: PetscObjectOptionsBegin((PetscObject)pc);
911: PetscCall(PCSetFromOptions_MG(pc, PetscOptionsObject));
912: PetscOptionsEnd();
913: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
915: /* set cheby eigen estimates from SA to use in the solver */
916: if (pc_gamg->use_sa_esteig) {
917: for (lidx = 1, level = pc_gamg->Nlevels - 2; level >= 0; lidx++, level--) {
918: KSP smoother;
919: PetscBool ischeb;
921: PetscCall(PCMGGetSmoother(pc, lidx, &smoother));
922: PetscCall(PetscObjectTypeCompare((PetscObject)smoother, KSPCHEBYSHEV, &ischeb));
923: if (ischeb) {
924: KSP_Chebyshev *cheb = (KSP_Chebyshev *)smoother->data;
926: // The command line will override these settings because KSPSetFromOptions is called in PCSetUp_MG
927: if (mg->max_eigen_DinvA[level] > 0) {
928: // SA uses Jacobi for P; we use SA estimates if the smoother is also Jacobi or if the user explicitly requested it.
929: // TODO: This should test whether it's the same Jacobi variant (DIAG, ROWSUM, etc.)
930: PetscReal emax, emin;
932: emin = mg->min_eigen_DinvA[level];
933: emax = mg->max_eigen_DinvA[level];
934: 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));
935: cheb->emin_provided = emin;
936: cheb->emax_provided = emax;
937: }
938: }
939: }
940: }
942: PetscCall(PCSetUp_MG(pc));
944: /* clean up */
945: for (level = 1; level < pc_gamg->Nlevels; level++) {
946: PetscCall(MatDestroy(&Parr[level]));
947: PetscCall(MatDestroy(&Aarr[level]));
948: }
949: } else {
950: KSP smoother;
952: PetscCall(PetscInfo(pc, "%s: One level solver used (system is seen as DD). Using default solver.\n", ((PetscObject)pc)->prefix));
953: PetscCall(PCMGGetSmoother(pc, 0, &smoother));
954: PetscCall(KSPSetOperators(smoother, Aarr[0], Aarr[0]));
955: PetscCall(KSPSetType(smoother, KSPPREONLY));
956: PetscCall(PCSetUp_MG(pc));
957: }
958: PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_SETUP], 0, 0, 0, 0));
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: /*
963: PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner
964: that was created with PCCreate_GAMG().
966: Input Parameter:
967: . pc - the preconditioner context
969: Application Interface Routine: PCDestroy()
970: */
971: PetscErrorCode PCDestroy_GAMG(PC pc)
972: {
973: PC_MG *mg = (PC_MG *)pc->data;
974: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
976: PetscFunctionBegin;
977: PetscCall(PCReset_GAMG(pc));
978: if (pc_gamg->ops->destroy) PetscCall((*pc_gamg->ops->destroy)(pc));
979: PetscCall(PetscFree(pc_gamg->ops));
980: PetscCall(PetscFree(pc_gamg->gamg_type_name));
981: PetscCall(PetscFree(pc_gamg));
982: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetProcEqLim_C", NULL));
983: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseEqLim_C", NULL));
984: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRepartition_C", NULL));
985: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetEigenvalues_C", NULL));
986: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetUseSAEstEig_C", NULL));
987: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRecomputeEstEig_C", NULL));
988: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetReuseInterpolation_C", NULL));
989: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetUseAggs_C", NULL));
990: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetParallelCoarseGridSolve_C", NULL));
991: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCpuPinCoarseGrids_C", NULL));
992: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetCoarseGridLayoutType_C", NULL));
993: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThreshold_C", NULL));
994: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetRankReductionFactors_C", NULL));
995: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetThresholdScale_C", NULL));
996: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetType_C", NULL));
997: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGGetType_C", NULL));
998: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNlevels_C", NULL));
999: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGASMSetHEM_C", NULL));
1000: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndices_C", NULL));
1001: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetInjectionIndex_C", NULL));
1002: PetscCall(PCDestroy_MG(pc));
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: /*@
1007: PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction in `PCGAMG`
1009: Logically Collective
1011: Input Parameters:
1012: + pc - the preconditioner context
1013: - n - the number of equations
1015: Options Database Key:
1016: . -pc_gamg_process_eq_limit limit - set the limit
1018: Level: intermediate
1020: Note:
1021: `PCGAMG` will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process
1022: that has degrees of freedom
1024: Developer Note:
1025: Should be named `PCGAMGSetProcessEquationLimit()`.
1027: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()`
1028: @*/
1029: PetscErrorCode PCGAMGSetProcEqLim(PC pc, PetscInt n)
1030: {
1031: PetscFunctionBegin;
1033: PetscTryMethod(pc, "PCGAMGSetProcEqLim_C", (PC, PetscInt), (pc, n));
1034: PetscFunctionReturn(PETSC_SUCCESS);
1035: }
1037: static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n)
1038: {
1039: PC_MG *mg = (PC_MG *)pc->data;
1040: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1042: PetscFunctionBegin;
1043: if (n > 0) pc_gamg->min_eq_proc = n;
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: /*@
1048: PCGAMGSetCoarseEqLim - Set maximum number of equations on the coarsest grid of `PCGAMG`
1050: Collective
1052: Input Parameters:
1053: + pc - the preconditioner context
1054: - n - maximum number of equations to aim for
1056: Options Database Key:
1057: . -pc_gamg_coarse_eq_limit limit - set the limit
1059: Level: intermediate
1061: Note:
1062: For example -pc_gamg_coarse_eq_limit 1000 will stop coarsening once the coarse grid
1063: has less than 1000 unknowns.
1065: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`, `PCGAMGSetRepartition()`,
1066: `PCGAMGSetParallelCoarseGridSolve()`
1067: @*/
1068: PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n)
1069: {
1070: PetscFunctionBegin;
1072: PetscTryMethod(pc, "PCGAMGSetCoarseEqLim_C", (PC, PetscInt), (pc, n));
1073: PetscFunctionReturn(PETSC_SUCCESS);
1074: }
1076: static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n)
1077: {
1078: PC_MG *mg = (PC_MG *)pc->data;
1079: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1081: PetscFunctionBegin;
1082: if (n > 0) pc_gamg->coarse_eq_limit = n;
1083: PetscFunctionReturn(PETSC_SUCCESS);
1084: }
1086: /*@
1087: PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids when reducing the number of MPI processes used
1089: Collective
1091: Input Parameters:
1092: + pc - the preconditioner context
1093: - n - `PETSC_TRUE` or `PETSC_FALSE`
1095: Options Database Key:
1096: . -pc_gamg_repartition (true|false) - turn on the repartitioning
1098: Level: intermediate
1100: Note:
1101: This will generally improve the load balance of the work on each level so the solves will be faster but it increases the
1102: preconditioner setup time.
1104: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetRankReductionFactors()`
1105: @*/
1106: PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n)
1107: {
1108: PetscFunctionBegin;
1110: PetscTryMethod(pc, "PCGAMGSetRepartition_C", (PC, PetscBool), (pc, n));
1111: PetscFunctionReturn(PETSC_SUCCESS);
1112: }
1114: static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n)
1115: {
1116: PC_MG *mg = (PC_MG *)pc->data;
1117: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1119: PetscFunctionBegin;
1120: pc_gamg->repart = n;
1121: PetscFunctionReturn(PETSC_SUCCESS);
1122: }
1124: /*@
1125: PCGAMGSetUseSAEstEig - Use the eigenvalue estimate from smoothed aggregation for the Chebyshev smoother during the solution process
1127: Collective
1129: Input Parameters:
1130: + pc - the preconditioner context
1131: - b - flag
1133: Options Database Key:
1134: . -pc_gamg_use_sa_esteig (true|false) - use the eigen estimate
1136: Level: advanced
1138: Notes:
1139: 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$.
1140: Eigenvalue estimates (based on a few `PCCG` or `PCGMRES` iterations) are computed to choose $\omega$ so that this is a stable smoothing operation.
1141: If `KSPCHEBYSHEV` with `PCJACOBI` (diagonal) preconditioning is used for smoothing on the finest level, then the eigenvalue estimates
1142: can be reused during the solution process.
1143: This option is only used when the smoother uses `PCJACOBI`, and should be turned off when a different `PCJacobiType` is used.
1145: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `PCGAMGSetRecomputeEstEig()`
1146: @*/
1147: PetscErrorCode PCGAMGSetUseSAEstEig(PC pc, PetscBool b)
1148: {
1149: PetscFunctionBegin;
1151: PetscTryMethod(pc, "PCGAMGSetUseSAEstEig_C", (PC, PetscBool), (pc, b));
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: static PetscErrorCode PCGAMGSetUseSAEstEig_GAMG(PC pc, PetscBool b)
1156: {
1157: PC_MG *mg = (PC_MG *)pc->data;
1158: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1160: PetscFunctionBegin;
1161: pc_gamg->use_sa_esteig = b;
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: }
1165: /*@
1166: PCGAMGSetRecomputeEstEig - Set flag for Chebyshev smoothers to recompute the eigenvalue estimates when a new matrix is used
1168: Collective
1170: Input Parameters:
1171: + pc - the preconditioner context
1172: - b - flag, default is `PETSC_TRUE`
1174: Options Database Key:
1175: . -pc_gamg_recompute_esteig (true|false) - recompute the eigenvalue estimate
1177: Level: advanced
1179: Note:
1180: If the matrix changes only slightly in a new solve using `PETSC_FALSE` will save time in setting up the preconditioner
1181: and may not affect the solution time much.
1183: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`
1184: @*/
1185: PetscErrorCode PCGAMGSetRecomputeEstEig(PC pc, PetscBool b)
1186: {
1187: PetscFunctionBegin;
1189: PetscTryMethod(pc, "PCGAMGSetRecomputeEstEig_C", (PC, PetscBool), (pc, b));
1190: PetscFunctionReturn(PETSC_SUCCESS);
1191: }
1193: static PetscErrorCode PCGAMGSetRecomputeEstEig_GAMG(PC pc, PetscBool b)
1194: {
1195: PC_MG *mg = (PC_MG *)pc->data;
1196: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1198: PetscFunctionBegin;
1199: pc_gamg->recompute_esteig = b;
1200: PetscFunctionReturn(PETSC_SUCCESS);
1201: }
1203: /*@
1204: PCGAMGSetEigenvalues - Set WHAT eigenvalues WHY?
1206: Collective
1208: Input Parameters:
1209: + pc - the preconditioner context
1210: . emax - max eigenvalue
1211: - emin - min eigenvalue
1213: Options Database Key:
1214: . -pc_gamg_eigenvalues emin,emax - estimates of the eigenvalues
1216: Level: intermediate
1218: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetUseSAEstEig()`
1219: @*/
1220: PetscErrorCode PCGAMGSetEigenvalues(PC pc, PetscReal emax, PetscReal emin)
1221: {
1222: PetscFunctionBegin;
1224: PetscTryMethod(pc, "PCGAMGSetEigenvalues_C", (PC, PetscReal, PetscReal), (pc, emax, emin));
1225: PetscFunctionReturn(PETSC_SUCCESS);
1226: }
1228: static PetscErrorCode PCGAMGSetEigenvalues_GAMG(PC pc, PetscReal emax, PetscReal emin)
1229: {
1230: PC_MG *mg = (PC_MG *)pc->data;
1231: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1233: PetscFunctionBegin;
1234: 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);
1235: 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);
1236: pc_gamg->emax = emax;
1237: pc_gamg->emin = emin;
1238: PetscFunctionReturn(PETSC_SUCCESS);
1239: }
1241: /*@
1242: PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding a `PCGAMG` algebraic multigrid preconditioner when the matrix has changed
1244: Collective
1246: Input Parameters:
1247: + pc - the preconditioner context
1248: - n - `PETSC_TRUE` or `PETSC_FALSE`
1250: Options Database Key:
1251: . -pc_gamg_reuse_interpolation (true|false) - reuse the previous interpolation
1253: Level: intermediate
1255: Note:
1256: May negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows
1257: rebuilding the preconditioner quicker.
1259: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
1260: @*/
1261: PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n)
1262: {
1263: PetscFunctionBegin;
1265: PetscTryMethod(pc, "PCGAMGSetReuseInterpolation_C", (PC, PetscBool), (pc, n));
1266: PetscFunctionReturn(PETSC_SUCCESS);
1267: }
1269: static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n)
1270: {
1271: PC_MG *mg = (PC_MG *)pc->data;
1272: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1274: PetscFunctionBegin;
1275: pc_gamg->reuse_prol = n;
1276: PetscFunctionReturn(PETSC_SUCCESS);
1277: }
1279: /*@
1280: PCGAMGASMSetUseAggs - Have the `PCGAMG` smoother on each level use `PCASM` where the aggregates defined by the coarsening process are
1281: used as the subdomains for the additive Schwarz preconditioner smoother
1283: Collective
1285: Input Parameters:
1286: + pc - the preconditioner context
1287: - flg - `PETSC_TRUE` to use aggregates, `PETSC_FALSE` to not
1289: Options Database Key:
1290: . -pc_gamg_asm_use_agg (true|false) - use aggregates to define the additive Schwarz subdomains
1292: Level: intermediate
1294: Note:
1295: This option automatically sets the preconditioner on the levels to be `PCASM`.
1297: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCASM`, `PCSetType`
1298: @*/
1299: PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg)
1300: {
1301: PetscFunctionBegin;
1303: PetscTryMethod(pc, "PCGAMGASMSetUseAggs_C", (PC, PetscBool), (pc, flg));
1304: PetscFunctionReturn(PETSC_SUCCESS);
1305: }
1307: static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg)
1308: {
1309: PC_MG *mg = (PC_MG *)pc->data;
1310: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1312: PetscFunctionBegin;
1313: pc_gamg->use_aggs_in_asm = flg;
1314: PetscFunctionReturn(PETSC_SUCCESS);
1315: }
1317: /*@
1318: PCGAMGSetParallelCoarseGridSolve - allow a parallel coarse grid solver
1320: Collective
1322: Input Parameters:
1323: + pc - the preconditioner context
1324: - flg - `PETSC_TRUE` to not force coarse grid onto one processor
1326: Options Database Key:
1327: . -pc_gamg_parallel_coarse_grid_solver (true|false) - use a parallel coarse grid solver
1329: Level: intermediate
1331: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetCpuPinCoarseGrids()`, `PCGAMGSetRankReductionFactors()`
1332: @*/
1333: PetscErrorCode PCGAMGSetParallelCoarseGridSolve(PC pc, PetscBool flg)
1334: {
1335: PetscFunctionBegin;
1337: PetscTryMethod(pc, "PCGAMGSetParallelCoarseGridSolve_C", (PC, PetscBool), (pc, flg));
1338: PetscFunctionReturn(PETSC_SUCCESS);
1339: }
1341: static PetscErrorCode PCGAMGSetParallelCoarseGridSolve_GAMG(PC pc, PetscBool flg)
1342: {
1343: PC_MG *mg = (PC_MG *)pc->data;
1344: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1346: PetscFunctionBegin;
1347: pc_gamg->use_parallel_coarse_grid_solver = flg;
1348: PetscFunctionReturn(PETSC_SUCCESS);
1349: }
1351: /*@
1352: 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
1354: Collective
1356: Input Parameters:
1357: + pc - the preconditioner context
1358: - flg - `PETSC_TRUE` to pin coarse grids to the CPU
1360: Options Database Key:
1361: . -pc_gamg_cpu_pin_coarse_grids (true|false) - pin the coarse grids to the CPU
1363: Level: intermediate
1365: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetCoarseGridLayoutType()`, `PCGAMGSetParallelCoarseGridSolve()`
1366: @*/
1367: PetscErrorCode PCGAMGSetCpuPinCoarseGrids(PC pc, PetscBool flg)
1368: {
1369: PetscFunctionBegin;
1371: PetscTryMethod(pc, "PCGAMGSetCpuPinCoarseGrids_C", (PC, PetscBool), (pc, flg));
1372: PetscFunctionReturn(PETSC_SUCCESS);
1373: }
1375: static PetscErrorCode PCGAMGSetCpuPinCoarseGrids_GAMG(PC pc, PetscBool flg)
1376: {
1377: PC_MG *mg = (PC_MG *)pc->data;
1378: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1380: PetscFunctionBegin;
1381: pc_gamg->cpu_pin_coarse_grids = flg;
1382: PetscFunctionReturn(PETSC_SUCCESS);
1383: }
1385: /*@
1386: PCGAMGSetCoarseGridLayoutType - place coarse grids on processors with natural order (compact type)
1388: Collective
1390: Input Parameters:
1391: + pc - the preconditioner context
1392: - flg - `PCGAMGLayoutType` type, either `PCGAMG_LAYOUT_COMPACT` or `PCGAMG_LAYOUT_SPREAD`
1394: Options Database Key:
1395: . -pc_gamg_coarse_grid_layout_type (compact|spread) - place the coarse grids with natural ordering
1397: Level: intermediate
1399: .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`
1400: @*/
1401: PetscErrorCode PCGAMGSetCoarseGridLayoutType(PC pc, PCGAMGLayoutType flg)
1402: {
1403: PetscFunctionBegin;
1405: PetscTryMethod(pc, "PCGAMGSetCoarseGridLayoutType_C", (PC, PCGAMGLayoutType), (pc, flg));
1406: PetscFunctionReturn(PETSC_SUCCESS);
1407: }
1409: static PetscErrorCode PCGAMGSetCoarseGridLayoutType_GAMG(PC pc, PCGAMGLayoutType flg)
1410: {
1411: PC_MG *mg = (PC_MG *)pc->data;
1412: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1414: PetscFunctionBegin;
1415: pc_gamg->layout_type = flg;
1416: PetscFunctionReturn(PETSC_SUCCESS);
1417: }
1419: /*@
1420: PCGAMGSetNlevels - Sets the maximum number of levels `PCGAMG` will use
1422: Collective
1424: Input Parameters:
1425: + pc - the preconditioner
1426: - n - the maximum number of levels to use
1428: Options Database Key:
1429: . -pc_mg_levels n - set the maximum number of levels to allow
1431: Level: intermediate
1433: Developer Notes:
1434: Should be called `PCGAMGSetMaximumNumberlevels()` and possible be shared with `PCMG`
1436: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
1437: @*/
1438: PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n)
1439: {
1440: PetscFunctionBegin;
1442: PetscTryMethod(pc, "PCGAMGSetNlevels_C", (PC, PetscInt), (pc, n));
1443: PetscFunctionReturn(PETSC_SUCCESS);
1444: }
1446: static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n)
1447: {
1448: PC_MG *mg = (PC_MG *)pc->data;
1449: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1451: PetscFunctionBegin;
1452: pc_gamg->Nlevels = n;
1453: PetscFunctionReturn(PETSC_SUCCESS);
1454: }
1456: /*@
1457: PCGAMGASMSetHEM - Sets the number of HEM matching passed
1459: Collective
1461: Input Parameters:
1462: + pc - the preconditioner
1463: - n - number of HEM matching passed to construct ASM subdomains
1465: Options Database Key:
1466: . -pc_gamg_asm_hem n - set the number of HEM matching passed
1468: Level: intermediate
1470: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`
1471: @*/
1472: PetscErrorCode PCGAMGASMSetHEM(PC pc, PetscInt n)
1473: {
1474: PetscFunctionBegin;
1476: PetscTryMethod(pc, "PCGAMGASMSetHEM_C", (PC, PetscInt), (pc, n));
1477: PetscFunctionReturn(PETSC_SUCCESS);
1478: }
1480: static PetscErrorCode PCGAMGASMSetHEM_GAMG(PC pc, PetscInt n)
1481: {
1482: PC_MG *mg = (PC_MG *)pc->data;
1483: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1485: PetscFunctionBegin;
1486: pc_gamg->asm_hem_aggs = n;
1487: PetscFunctionReturn(PETSC_SUCCESS);
1488: }
1490: /*@
1491: PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph
1493: Not Collective
1495: Input Parameters:
1496: + pc - the preconditioner context
1497: . 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
1498: - n - number of threshold values provided in array
1500: Options Database Key:
1501: . -pc_gamg_threshold l0,l1,... - the thresholds to drop edges
1503: Level: intermediate
1505: Notes:
1506: Before coarsening or aggregating the graph, `PCGAMG` removes small values from the graph with this threshold, reducing the coupling
1507: in the graph and yielding a different (perhaps better) coarser set of points.
1509: Increasing the threshold decreases the rate of coarsening. Conversely reducing the threshold increases the rate of coarsening (aggressive coarsening)
1510: and thereby reduces the complexity of the coarse grids, and generally results in slower solver convergence rates.
1512: If `n` is less than the total number of coarsenings (see `PCGAMGSetNlevels()`), then threshold scaling (see `PCGAMGSetThresholdScale()`) is used for each successive coarsening.
1513: In this case, `PCGAMGSetThresholdScale()` must be called before `PCGAMGSetThreshold()`.
1515: If `n` is greater than the total number of levels, the excess entries in threshold are ignored.
1517: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`,
1518: `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetThresholdScale()`
1519: @*/
1520: PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal v[], PetscInt n)
1521: {
1522: PetscFunctionBegin;
1524: if (n) PetscAssertPointer(v, 2);
1525: PetscTryMethod(pc, "PCGAMGSetThreshold_C", (PC, PetscReal[], PetscInt), (pc, v, n));
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal v[], PetscInt n)
1530: {
1531: PC_MG *mg = (PC_MG *)pc->data;
1532: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1533: PetscInt i;
1535: PetscFunctionBegin;
1536: for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->threshold[i] = v[i];
1537: for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale;
1538: PetscFunctionReturn(PETSC_SUCCESS);
1539: }
1541: /*@
1542: PCGAMGSetRankReductionFactors - Set a manual schedule for MPI process reduction on coarse grids
1544: Collective
1546: Input Parameters:
1547: + pc - the preconditioner context
1548: . v - array of reduction factors. 0 for first value forces a reduction to one process/device on first level in CUDA
1549: - n - number of values provided in array
1551: Options Database Key:
1552: . -pc_gamg_rank_reduction_factors r0,r1,... - provide the schedule
1554: Level: intermediate
1556: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetParallelCoarseGridSolve()`
1557: @*/
1558: PetscErrorCode PCGAMGSetRankReductionFactors(PC pc, PetscInt v[], PetscInt n)
1559: {
1560: PetscFunctionBegin;
1562: if (n) PetscAssertPointer(v, 2);
1563: PetscTryMethod(pc, "PCGAMGSetRankReductionFactors_C", (PC, PetscInt[], PetscInt), (pc, v, n));
1564: PetscFunctionReturn(PETSC_SUCCESS);
1565: }
1567: static PetscErrorCode PCGAMGSetRankReductionFactors_GAMG(PC pc, PetscInt v[], PetscInt n)
1568: {
1569: PC_MG *mg = (PC_MG *)pc->data;
1570: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1571: PetscInt i;
1573: PetscFunctionBegin;
1574: for (i = 0; i < PetscMin(n, PETSC_MG_MAXLEVELS); i++) pc_gamg->level_reduction_factors[i] = v[i];
1575: for (; i < PETSC_MG_MAXLEVELS; i++) pc_gamg->level_reduction_factors[i] = -1; /* 0 stop putting one process/device on first level */
1576: PetscFunctionReturn(PETSC_SUCCESS);
1577: }
1579: /*@
1580: PCGAMGSetThresholdScale - Relative threshold reduction at each level
1582: Not Collective
1584: Input Parameters:
1585: + pc - the preconditioner context
1586: - v - the threshold value reduction, usually < 1.0
1588: Options Database Key:
1589: . -pc_gamg_threshold_scale v - set the relative threshold reduction on each level
1591: Level: advanced
1593: Note:
1594: The initial threshold (for an arbitrary number of levels starting from the finest) can be set with `PCGAMGSetThreshold()`.
1595: This scaling is used for each subsequent coarsening, but must be called before `PCGAMGSetThreshold()`.
1597: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`
1598: @*/
1599: PetscErrorCode PCGAMGSetThresholdScale(PC pc, PetscReal v)
1600: {
1601: PetscFunctionBegin;
1603: PetscTryMethod(pc, "PCGAMGSetThresholdScale_C", (PC, PetscReal), (pc, v));
1604: PetscFunctionReturn(PETSC_SUCCESS);
1605: }
1607: static PetscErrorCode PCGAMGSetThresholdScale_GAMG(PC pc, PetscReal v)
1608: {
1609: PC_MG *mg = (PC_MG *)pc->data;
1610: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1612: PetscFunctionBegin;
1613: pc_gamg->threshold_scale = v;
1614: PetscFunctionReturn(PETSC_SUCCESS);
1615: }
1617: /*@
1618: PCGAMGSetType - Set the type of algorithm `PCGAMG` should use
1620: Collective
1622: Input Parameters:
1623: + pc - the preconditioner context
1624: - type - `PCGAMGAGG`, `PCGAMGGEO`, or `PCGAMGCLASSICAL`
1626: Options Database Key:
1627: . -pc_gamg_type (agg|geo|classical) - type of algebraic multigrid to apply - only agg is supported
1629: Level: intermediate
1631: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMGGetType()`, `PCGAMG`, `PCGAMGType`
1632: @*/
1633: PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type)
1634: {
1635: PetscFunctionBegin;
1637: PetscTryMethod(pc, "PCGAMGSetType_C", (PC, PCGAMGType), (pc, type));
1638: PetscFunctionReturn(PETSC_SUCCESS);
1639: }
1641: /*@
1642: PCGAMGGetType - Get the type of algorithm `PCGAMG` will use
1644: Collective
1646: Input Parameter:
1647: . pc - the preconditioner context
1649: Output Parameter:
1650: . type - the type of algorithm used
1652: Level: intermediate
1654: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetType()`, `PCGAMGType`
1655: @*/
1656: PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type)
1657: {
1658: PetscFunctionBegin;
1660: PetscUseMethod(pc, "PCGAMGGetType_C", (PC, PCGAMGType *), (pc, type));
1661: PetscFunctionReturn(PETSC_SUCCESS);
1662: }
1664: static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type)
1665: {
1666: PC_MG *mg = (PC_MG *)pc->data;
1667: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1669: PetscFunctionBegin;
1670: *type = pc_gamg->type;
1671: PetscFunctionReturn(PETSC_SUCCESS);
1672: }
1674: static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type)
1675: {
1676: PC_MG *mg = (PC_MG *)pc->data;
1677: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1678: PetscErrorCode (*r)(PC);
1680: PetscFunctionBegin;
1681: pc_gamg->type = type;
1682: PetscCall(PetscFunctionListFind(GAMGList, type, &r));
1683: PetscCheck(r, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown GAMG type %s given", type);
1684: if (pc_gamg->ops->destroy) {
1685: PetscCall((*pc_gamg->ops->destroy)(pc));
1686: PetscCall(PetscMemzero(pc_gamg->ops, sizeof(struct _PCGAMGOps)));
1687: pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1688: /* cleaning up common data in pc_gamg - this should disappear someday */
1689: pc_gamg->data_cell_cols = 0;
1690: pc_gamg->data_cell_rows = 0;
1691: pc_gamg->orig_data_cell_cols = 0;
1692: pc_gamg->orig_data_cell_rows = 0;
1693: PetscCall(PetscFree(pc_gamg->data));
1694: pc_gamg->data_sz = 0;
1695: }
1696: PetscCall(PetscFree(pc_gamg->gamg_type_name));
1697: PetscCall(PetscStrallocpy(type, &pc_gamg->gamg_type_name));
1698: PetscCall((*r)(pc));
1699: PetscFunctionReturn(PETSC_SUCCESS);
1700: }
1702: static PetscErrorCode PCView_GAMG(PC pc, PetscViewer viewer)
1703: {
1704: PC_MG *mg = (PC_MG *)pc->data;
1705: PC_MG_Levels **mglevels = mg->levels;
1706: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1707: PetscReal gc, oc;
1709: PetscFunctionBegin;
1710: PetscCall(PetscViewerASCIIPrintf(viewer, " GAMG specific options\n"));
1711: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for dropping small values in graph on each level ="));
1712: for (PetscInt i = 0; i < mg->nlevels; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)pc_gamg->threshold[i]));
1713: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1714: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold scaling factor for each level not specified = %g\n", (double)pc_gamg->threshold_scale));
1715: 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
1716: else if (pc_gamg->asm_hem_aggs) {
1717: 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));
1718: PetscCall(PetscViewerASCIIPushTab(viewer));
1719: PetscCall(PetscViewerASCIIPushTab(viewer));
1720: PetscCall(PetscViewerASCIIPushTab(viewer));
1721: PetscCall(PetscViewerASCIIPushTab(viewer));
1722: PetscCall(MatCoarsenView(pc_gamg->asm_crs, viewer));
1723: PetscCall(PetscViewerASCIIPopTab(viewer));
1724: PetscCall(PetscViewerASCIIPopTab(viewer));
1725: PetscCall(PetscViewerASCIIPopTab(viewer));
1726: }
1727: 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"));
1728: if (pc_gamg->injection_index_size) {
1729: PetscCall(PetscViewerASCIIPrintf(viewer, " Using injection restriction/prolongation on first level, dofs:"));
1730: for (int i = 0; i < pc_gamg->injection_index_size; i++) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, pc_gamg->injection_index[i]));
1731: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1732: }
1733: if (pc_gamg->ops->view) PetscCall((*pc_gamg->ops->view)(pc, viewer));
1734: gc = oc = 0;
1735: PetscCall(PCMGGetGridComplexity(pc, &gc, &oc));
1736: PetscCall(PetscViewerASCIIPrintf(viewer, " Complexity: grid = %g operator = %g\n", (double)gc, (double)oc));
1737: PetscCall(PetscViewerASCIIPrintf(viewer, " Per-level complexity: op = operator, int = interpolation\n"));
1738: PetscCall(PetscViewerASCIIPrintf(viewer, " #equations | #active PEs | avg nnz/row op | avg nnz/row int\n"));
1739: for (PetscInt i = 0; i < mg->nlevels; i++) {
1740: MatInfo info;
1741: Mat A;
1742: PetscReal rd[3];
1743: PetscInt rst, ren, N;
1745: PetscCall(KSPGetOperators(mglevels[i]->smoothd, NULL, &A));
1746: PetscCall(MatGetOwnershipRange(A, &rst, &ren));
1747: PetscCall(MatGetSize(A, &N, NULL));
1748: PetscCall(MatGetInfo(A, MAT_LOCAL, &info));
1749: rd[0] = (ren - rst > 0) ? 1 : 0;
1750: rd[1] = info.nz_used;
1751: rd[2] = 0;
1752: if (i) {
1753: Mat P;
1754: PetscCall(PCMGGetInterpolation(pc, i, &P));
1755: PetscCall(MatGetInfo(P, MAT_LOCAL, &info));
1756: rd[2] = info.nz_used;
1757: }
1758: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, rd, 3, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)pc)));
1759: 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)));
1760: }
1761: PetscFunctionReturn(PETSC_SUCCESS);
1762: }
1764: /*@
1765: PCGAMGSetInjectionIndex - Array of subset of variables per vertex to inject into coarse grid space
1767: Logically Collective
1769: Input Parameters:
1770: + pc - the coarsen context
1771: . n - number of indices
1772: - idx - array of indices
1774: Options Database Key:
1775: . -pc_gamg_injection_index - array of subset of variables per vertex to use for injection coarse grid space
1777: Level: intermediate
1779: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`
1780: @*/
1781: PetscErrorCode PCGAMGSetInjectionIndex(PC pc, PetscInt n, PetscInt idx[])
1782: {
1783: PetscFunctionBegin;
1786: PetscTryMethod(pc, "PCGAMGSetInjectionIndex_C", (PC, PetscInt, PetscInt[]), (pc, n, idx));
1787: PetscFunctionReturn(PETSC_SUCCESS);
1788: }
1790: static PetscErrorCode PCGAMGSetInjectionIndex_GAMG(PC pc, PetscInt n, PetscInt idx[])
1791: {
1792: PC_MG *mg = (PC_MG *)pc->data;
1793: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1795: PetscFunctionBegin;
1796: pc_gamg->injection_index_size = n;
1797: 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);
1798: for (PetscInt i = 0; i < n; i++) pc_gamg->injection_index[i] = idx[i];
1799: PetscFunctionReturn(PETSC_SUCCESS);
1800: }
1802: static PetscErrorCode PCSetFromOptions_GAMG(PC pc, PetscOptionItems PetscOptionsObject)
1803: {
1804: PC_MG *mg = (PC_MG *)pc->data;
1805: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
1806: PetscBool flag;
1807: MPI_Comm comm;
1808: char prefix[256], tname[32];
1809: PetscInt i, n;
1810: const char *pcpre;
1811: static const char *LayoutTypes[] = {"compact", "spread", "PCGAMGLayoutType", "PC_GAMG_LAYOUT", NULL};
1813: PetscFunctionBegin;
1814: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1815: PetscOptionsHeadBegin(PetscOptionsObject, "GAMG options");
1816: PetscCall(PetscOptionsInt("-pc_mg_levels", "Number of Levels", "PCMGSetLevels", -1, &n, &flag));
1817: PetscCheck(!flag, comm, PETSC_ERR_ARG_WRONG, "Invalid flag -pc_mg_levels. GAMG does not allow the number of levels to be set.");
1818: 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));
1819: if (flag) PetscCall(PCGAMGSetType(pc, tname));
1820: PetscCall(PetscOptionsBool("-pc_gamg_repartition", "Repartion coarse grids", "PCGAMGSetRepartition", pc_gamg->repart, &pc_gamg->repart, NULL));
1821: 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));
1822: 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));
1823: PetscCall(PetscOptionsBool("-pc_gamg_reuse_interpolation", "Reuse prolongation operator", "PCGAMGReuseInterpolation", pc_gamg->reuse_prol, &pc_gamg->reuse_prol, NULL));
1824: 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));
1825: 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));
1826: 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));
1827: 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,
1828: (PetscEnum)pc_gamg->layout_type, (PetscEnum *)&pc_gamg->layout_type, NULL));
1829: 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));
1830: 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));
1831: 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));
1832: PetscCall(PetscOptionsReal("-pc_gamg_threshold_scale", "Scaling of threshold for each level not specified", "PCGAMGSetThresholdScale", pc_gamg->threshold_scale, &pc_gamg->threshold_scale, NULL));
1833: n = PETSC_MG_MAXLEVELS;
1834: PetscCall(PetscOptionsRealArray("-pc_gamg_threshold", "Relative threshold to use for dropping edges in aggregation graph", "PCGAMGSetThreshold", pc_gamg->threshold, &n, &flag));
1835: if (!flag || n < PETSC_MG_MAXLEVELS) {
1836: if (!flag) n = 1;
1837: i = n;
1838: do {
1839: pc_gamg->threshold[i] = pc_gamg->threshold[i - 1] * pc_gamg->threshold_scale;
1840: } while (++i < PETSC_MG_MAXLEVELS);
1841: }
1842: PetscCall(PetscOptionsInt("-pc_mg_levels", "Set number of MG levels (should get from base class)", "PCGAMGSetNlevels", pc_gamg->Nlevels, &pc_gamg->Nlevels, NULL));
1843: 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);
1844: n = PETSC_MG_MAXLEVELS;
1845: 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));
1846: if (!flag) i = 0;
1847: else i = n;
1848: do {
1849: pc_gamg->level_reduction_factors[i] = -1;
1850: } while (++i < PETSC_MG_MAXLEVELS);
1851: {
1852: PetscReal eminmax[2] = {0., 0.};
1853: n = 2;
1854: PetscCall(PetscOptionsRealArray("-pc_gamg_eigenvalues", "extreme eigenvalues for smoothed aggregation", "PCGAMGSetEigenvalues", eminmax, &n, &flag));
1855: if (flag) {
1856: PetscCheck(n == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_gamg_eigenvalues: must specify 2 parameters, min and max eigenvalues");
1857: PetscCall(PCGAMGSetEigenvalues(pc, eminmax[1], eminmax[0]));
1858: }
1859: }
1860: pc_gamg->injection_index_size = MAT_COARSEN_STRENGTH_INDEX_SIZE;
1861: 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));
1862: /* set options for subtype */
1863: PetscCall((*pc_gamg->ops->setfromoptions)(pc, PetscOptionsObject));
1865: PetscCall(PCGetOptionsPrefix(pc, &pcpre));
1866: PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%spc_gamg_", pcpre ? pcpre : ""));
1867: PetscOptionsHeadEnd();
1868: PetscFunctionReturn(PETSC_SUCCESS);
1869: }
1871: /*MC
1872: PCGAMG - Geometric algebraic multigrid (AMG) preconditioner
1874: Options Database Keys:
1875: + -pc_gamg_type type (agg|geo|classical) - see `PCGAMGType`
1876: . -pc_gamg_repartition (true|false) - repartition the degrees of freedom across the coarse grids as they are determined
1877: . -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,
1878: this also forces using `-mg_levels_pc_type asm`
1879: . -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
1880: equations on each process that has degrees of freedom
1881: . -pc_gamg_coarse_eq_limit limit - set maximum number of equations on coarsest grid to aim for
1882: . -pc_gamg_reuse_interpolation (true|false) - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations (should always be true)
1883: . -pc_gamg_threshold l0,l1,... - before aggregating the graph `PCGAMG` will remove small values from the graph on each level (< 0 does no filtering)
1884: - -pc_gamg_threshold_scale scale - scaling of threshold on each coarser grid if not specified
1886: Options Database Keys for Aggregation:
1887: + -pc_gamg_agg_nsmooths nsmooth - number of smoothing steps to use with smooth aggregation to construct prolongation
1888: . -pc_gamg_aggressive_coarsening n - number of aggressive coarsening (MIS-2 or square graph) levels from finest.
1889: . -pc_gamg_aggressive_square_graph (true|false) - use square graph ($A^T A$) for coarsening. Otherwise, MIS-k (k=2) is used, see `PCGAMGMISkSetAggressive()`
1890: . -pc_gamg_square[_i]_mat_product_algorithm algorithm - `MatProductAlgorithm` to use when squaring the matrix for aggressive coarsening (on level i < `n`)
1891: . -pc_gamg_mis_k_minimum_degree_ordering (true|false) - use minimum degree ordering in greedy MIS algorithm
1892: . -pc_gamg_asm_hem_aggs n - number of HEM aggregation steps for `PCASM` smoother
1893: - -pc_gamg_aggressive_mis_k n - number (k) distance in MIS coarsening (>2 is 'aggressive')
1895: Options Database Keys for Multigrid:
1896: + -pc_mg_cycle_type (v|w) - see `PCMGSetCycleType()`
1897: . -pc_mg_distinct_smoothup - configure the up and down (pre and post) smoothers separately, see `PCMGSetDistinctSmoothUp()`
1898: . -pc_mg_type (additive|multiplicative|full|kaskade) - see `PCMGType`
1899: - -pc_mg_levels levels - number of levels of multigrid to use; `PCGAMG` has a heuristic to determine the number of levels so
1900: this is not usually used with `PCGAMG`
1902: Level: intermediate
1904: Notes:
1905: To obtain good performance for `PCGAMG` for vector valued problems you must
1906: call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point
1907: call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator
1909: The many options for `PCMG` also work directly for `PCGAMG` such as controlling the smoothers on each level etc.
1911: .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`,
1912: `MatSetBlockSize()`,
1913: `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`, `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`,
1914: `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`, `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`,
1915: `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`, `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()`
1916: M*/
1917: PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc)
1918: {
1919: PC_GAMG *pc_gamg;
1920: PC_MG *mg;
1922: PetscFunctionBegin;
1923: /* register AMG type */
1924: PetscCall(PCGAMGInitializePackage());
1926: /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */
1927: PetscCall(PCSetType(pc, PCMG));
1928: PetscCall(PetscObjectChangeTypeName((PetscObject)pc, PCGAMG));
1930: /* create a supporting struct and attach it to pc */
1931: PetscCall(PetscNew(&pc_gamg));
1932: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_EXTERNAL));
1933: mg = (PC_MG *)pc->data;
1934: mg->innerctx = pc_gamg;
1936: PetscCall(PetscNew(&pc_gamg->ops));
1938: /* these should be in subctx but repartitioning needs simple arrays */
1939: pc_gamg->prolongator_filter = 0.0;
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: PetscFunctionBegin;
2005: if (PCGAMGPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
2006: PCGAMGPackageInitialized = PETSC_TRUE;
2007: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGGEO, PCCreateGAMG_GEO));
2008: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGAGG, PCCreateGAMG_AGG));
2009: PetscCall(PetscFunctionListAdd(&GAMGList, PCGAMGCLASSICAL, PCCreateGAMG_Classical));
2010: PetscCall(PetscRegisterFinalize(PCGAMGFinalizePackage));
2012: /* general events */
2013: PetscCall(PetscLogEventRegister("PCSetUp_GAMG+", PC_CLASSID, &petsc_gamg_setup_events[GAMG_SETUP]));
2014: PetscCall(PetscLogEventRegister(" PCGAMGCreateG", PC_CLASSID, &petsc_gamg_setup_events[GAMG_GRAPH]));
2015: PetscCall(PetscLogEventRegister(" GAMG Coarsen", PC_CLASSID, &petsc_gamg_setup_events[GAMG_COARSEN]));
2016: PetscCall(PetscLogEventRegister(" GAMG MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[GAMG_MIS]));
2017: PetscCall(PetscLogEventRegister(" PCGAMGProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROL]));
2018: PetscCall(PetscLogEventRegister(" GAMG Prol-col", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLA]));
2019: PetscCall(PetscLogEventRegister(" GAMG Prol-lift", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PROLB]));
2020: PetscCall(PetscLogEventRegister(" PCGAMGOptProl", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPT]));
2021: PetscCall(PetscLogEventRegister(" GAMG smooth", PC_CLASSID, &petsc_gamg_setup_events[GAMG_OPTSM]));
2022: PetscCall(PetscLogEventRegister(" PCGAMGCreateL", PC_CLASSID, &petsc_gamg_setup_events[GAMG_LEVEL]));
2023: PetscCall(PetscLogEventRegister(" GAMG PtAP", PC_CLASSID, &petsc_gamg_setup_events[GAMG_PTAP]));
2024: PetscCall(PetscLogEventRegister(" GAMG Reduce", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REDUCE]));
2025: PetscCall(PetscLogEventRegister(" GAMG Repart", PC_CLASSID, &petsc_gamg_setup_events[GAMG_REPART]));
2026: /* PetscCall(PetscLogEventRegister(" GAMG Inv-Srt", PC_CLASSID, &petsc_gamg_setup_events[SET13])); */
2027: /* PetscCall(PetscLogEventRegister(" GAMG Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14])); */
2028: /* PetscCall(PetscLogEventRegister(" GAMG Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15])); */
2029: for (PetscInt l = 0; l < PETSC_MG_MAXLEVELS; l++) {
2030: char ename[32];
2032: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Squ l%02" PetscInt_FMT, l));
2033: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][0]));
2034: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Gal l%02" PetscInt_FMT, l));
2035: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][1]));
2036: PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCGAMG Opt l%02" PetscInt_FMT, l));
2037: PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &petsc_gamg_setup_matmat_events[l][2]));
2038: }
2039: #if defined(GAMG_STAGES)
2040: { /* create timer stages */
2041: char str[32];
2042: PetscCall(PetscSNPrintf(str, PETSC_STATIC_ARRAY_LENGTH(str), "GAMG Level %d", 0));
2043: PetscCall(PetscLogStageRegister(str, &gamg_stages[0]));
2044: }
2045: #endif
2046: PetscFunctionReturn(PETSC_SUCCESS);
2047: }
2049: /*@C
2050: PCGAMGFinalizePackage - This function frees everything from the `PCGAMG` package. It is
2051: called from `PetscFinalize()` automatically.
2053: Level: developer
2055: .seealso: [](ch_ksp), `PetscFinalize()`
2056: @*/
2057: PetscErrorCode PCGAMGFinalizePackage(void)
2058: {
2059: PetscFunctionBegin;
2060: PCGAMGPackageInitialized = PETSC_FALSE;
2061: PetscCall(PetscFunctionListDestroy(&GAMGList));
2062: PetscFunctionReturn(PETSC_SUCCESS);
2063: }
2065: /*@C
2066: PCGAMGRegister - Register a `PCGAMG` implementation.
2068: Input Parameters:
2069: + type - string that will be used as the name of the `PCGAMG` type.
2070: - create - function for creating the gamg context.
2072: Level: developer
2074: .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()`
2075: @*/
2076: PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC))
2077: {
2078: PetscFunctionBegin;
2079: PetscCall(PCGAMGInitializePackage());
2080: PetscCall(PetscFunctionListAdd(&GAMGList, type, create));
2081: PetscFunctionReturn(PETSC_SUCCESS);
2082: }
2084: /*@
2085: PCGAMGCreateGraph - Creates a graph that is used by the ``PCGAMGType`` in the coarsening process
2087: Input Parameters:
2088: + pc - the `PCGAMG`
2089: - A - the matrix, for any level
2091: Output Parameter:
2092: . G - the graph
2094: Level: advanced
2096: .seealso: [](ch_ksp), `PCGAMGType`, `PCGAMG`, `PCGAMGSetType()`
2097: @*/
2098: PetscErrorCode PCGAMGCreateGraph(PC pc, Mat A, Mat *G)
2099: {
2100: PC_MG *mg = (PC_MG *)pc->data;
2101: PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx;
2103: PetscFunctionBegin;
2104: PetscCall(pc_gamg->ops->creategraph(pc, A, G));
2105: PetscFunctionReturn(PETSC_SUCCESS);
2106: }