1 /* 2 GAMG geometric-algebric multigrid PC - Mark Adams 2011 3 */ 4 5 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 6 #include <petscblaslapack.h> 7 #include <petscdm.h> 8 #include <petsc/private/kspimpl.h> 9 10 typedef struct { 11 PetscInt nsmooths; // number of smoothing steps to construct prolongation 12 PetscInt aggressive_coarsening_levels; // number of aggressive coarsening levels (square or MISk) 13 PetscInt aggressive_mis_k; // the k in MIS-k 14 PetscBool use_aggressive_square_graph; 15 PetscBool use_minimum_degree_ordering; 16 PetscBool use_low_mem_filter; 17 PetscBool graph_symmetrize; 18 MatCoarsen crs; 19 } PC_GAMG_AGG; 20 21 /*@ 22 PCGAMGSetNSmooths - Set number of smoothing steps (1 is typical) used to construct the prolongation operator 23 24 Logically Collective 25 26 Input Parameters: 27 + pc - the preconditioner context 28 - n - the number of smooths 29 30 Options Database Key: 31 . -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use 32 33 Level: intermediate 34 35 Note: 36 This is a different concept from the number smoothing steps used during the linear solution process which 37 can be set with `-mg_levels_ksp_max_it` 38 39 Developer Note: 40 This should be named `PCGAMGAGGSetNSmooths()`. 41 42 .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCMG`, `PCGAMG` 43 @*/ 44 PetscErrorCode PCGAMGSetNSmooths(PC pc, PetscInt n) 45 { 46 PetscFunctionBegin; 47 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 48 PetscValidLogicalCollectiveInt(pc, n, 2); 49 PetscTryMethod(pc, "PCGAMGSetNSmooths_C", (PC, PetscInt), (pc, n)); 50 PetscFunctionReturn(PETSC_SUCCESS); 51 } 52 53 static PetscErrorCode PCGAMGSetNSmooths_AGG(PC pc, PetscInt n) 54 { 55 PC_MG *mg = (PC_MG *)pc->data; 56 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 57 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 58 59 PetscFunctionBegin; 60 pc_gamg_agg->nsmooths = n; 61 PetscFunctionReturn(PETSC_SUCCESS); 62 } 63 64 /*@ 65 PCGAMGSetAggressiveLevels - Use aggressive coarsening on first n levels 66 67 Logically Collective 68 69 Input Parameters: 70 + pc - the preconditioner context 71 - n - 0, 1 or more 72 73 Options Database Key: 74 . -pc_gamg_aggressive_coarsening <n,default = 1> - Number of levels on which to square the graph on before aggregating it 75 76 Level: intermediate 77 78 .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()` 79 @*/ 80 PetscErrorCode PCGAMGSetAggressiveLevels(PC pc, PetscInt n) 81 { 82 PetscFunctionBegin; 83 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 84 PetscValidLogicalCollectiveInt(pc, n, 2); 85 PetscTryMethod(pc, "PCGAMGSetAggressiveLevels_C", (PC, PetscInt), (pc, n)); 86 PetscFunctionReturn(PETSC_SUCCESS); 87 } 88 89 /*@ 90 PCGAMGMISkSetAggressive - Number (k) distance in MIS coarsening (>2 is 'aggressive') 91 92 Logically Collective 93 94 Input Parameters: 95 + pc - the preconditioner context 96 - n - 1 or more (default = 2) 97 98 Options Database Key: 99 . -pc_gamg_aggressive_mis_k <n,default=2> - Number (k) distance in MIS coarsening (>2 is 'aggressive') 100 101 Level: intermediate 102 103 .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()` 104 @*/ 105 PetscErrorCode PCGAMGMISkSetAggressive(PC pc, PetscInt n) 106 { 107 PetscFunctionBegin; 108 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 109 PetscValidLogicalCollectiveInt(pc, n, 2); 110 PetscTryMethod(pc, "PCGAMGMISkSetAggressive_C", (PC, PetscInt), (pc, n)); 111 PetscFunctionReturn(PETSC_SUCCESS); 112 } 113 114 /*@ 115 PCGAMGSetAggressiveSquareGraph - Use graph square A'A for aggressive coarsening, old method 116 117 Logically Collective 118 119 Input Parameters: 120 + pc - the preconditioner context 121 - b - default false - MIS-k is faster 122 123 Options Database Key: 124 . -pc_gamg_aggressive_square_graph <bool,default=false> - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening 125 126 Level: intermediate 127 128 .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGMISkSetMinDegreeOrdering()`, `PCGAMGSetLowMemoryFilter()` 129 @*/ 130 PetscErrorCode PCGAMGSetAggressiveSquareGraph(PC pc, PetscBool b) 131 { 132 PetscFunctionBegin; 133 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 134 PetscValidLogicalCollectiveBool(pc, b, 2); 135 PetscTryMethod(pc, "PCGAMGSetAggressiveSquareGraph_C", (PC, PetscBool), (pc, b)); 136 PetscFunctionReturn(PETSC_SUCCESS); 137 } 138 139 /*@ 140 PCGAMGMISkSetMinDegreeOrdering - Use minimum degree ordering in greedy MIS algorithm 141 142 Logically Collective 143 144 Input Parameters: 145 + pc - the preconditioner context 146 - b - default true 147 148 Options Database Key: 149 . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm 150 151 Level: intermediate 152 153 .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGSetLowMemoryFilter()` 154 @*/ 155 PetscErrorCode PCGAMGMISkSetMinDegreeOrdering(PC pc, PetscBool b) 156 { 157 PetscFunctionBegin; 158 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 159 PetscValidLogicalCollectiveBool(pc, b, 2); 160 PetscTryMethod(pc, "PCGAMGMISkSetMinDegreeOrdering_C", (PC, PetscBool), (pc, b)); 161 PetscFunctionReturn(PETSC_SUCCESS); 162 } 163 164 /*@ 165 PCGAMGSetLowMemoryFilter - Use low memory graph/matrix filter 166 167 Logically Collective 168 169 Input Parameters: 170 + pc - the preconditioner context 171 - b - default false 172 173 Options Database Key: 174 . -pc_gamg_low_memory_threshold_filter <bool,default=false> - Use low memory graph/matrix filter 175 176 Level: intermediate 177 178 .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, 179 `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()` 180 @*/ 181 PetscErrorCode PCGAMGSetLowMemoryFilter(PC pc, PetscBool b) 182 { 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 185 PetscValidLogicalCollectiveBool(pc, b, 2); 186 PetscTryMethod(pc, "PCGAMGSetLowMemoryFilter_C", (PC, PetscBool), (pc, b)); 187 PetscFunctionReturn(PETSC_SUCCESS); 188 } 189 190 /*@ 191 PCGAMGSetGraphSymmetrize - Set the flag to symmetrize the graph used in coarsening 192 193 Logically Collective 194 195 Input Parameters: 196 + pc - the preconditioner context 197 - b - default false 198 199 Options Database Key: 200 . -pc_gamg_graph_symmetrize <bool,default=false> - Symmetrize the graph 201 202 Level: intermediate 203 204 .seealso: [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), `PCGAMG`, `PCGAMGSetThreshold()`, `PCGAMGSetAggressiveLevels()`, 205 `PCGAMGMISkSetAggressive()`, `PCGAMGSetAggressiveSquareGraph()`, `PCGAMGMISkSetMinDegreeOrdering()` 206 @*/ 207 PetscErrorCode PCGAMGSetGraphSymmetrize(PC pc, PetscBool b) 208 { 209 PetscFunctionBegin; 210 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 211 PetscValidLogicalCollectiveBool(pc, b, 2); 212 PetscTryMethod(pc, "PCGAMGSetGraphSymmetrize_C", (PC, PetscBool), (pc, b)); 213 PetscFunctionReturn(PETSC_SUCCESS); 214 } 215 216 static PetscErrorCode PCGAMGSetAggressiveLevels_AGG(PC pc, PetscInt n) 217 { 218 PC_MG *mg = (PC_MG *)pc->data; 219 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 220 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 221 222 PetscFunctionBegin; 223 pc_gamg_agg->aggressive_coarsening_levels = n; 224 PetscFunctionReturn(PETSC_SUCCESS); 225 } 226 227 static PetscErrorCode PCGAMGMISkSetAggressive_AGG(PC pc, PetscInt n) 228 { 229 PC_MG *mg = (PC_MG *)pc->data; 230 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 231 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 232 233 PetscFunctionBegin; 234 pc_gamg_agg->aggressive_mis_k = n; 235 PetscFunctionReturn(PETSC_SUCCESS); 236 } 237 238 static PetscErrorCode PCGAMGSetAggressiveSquareGraph_AGG(PC pc, PetscBool b) 239 { 240 PC_MG *mg = (PC_MG *)pc->data; 241 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 242 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 243 244 PetscFunctionBegin; 245 pc_gamg_agg->use_aggressive_square_graph = b; 246 PetscFunctionReturn(PETSC_SUCCESS); 247 } 248 249 static PetscErrorCode PCGAMGSetLowMemoryFilter_AGG(PC pc, PetscBool b) 250 { 251 PC_MG *mg = (PC_MG *)pc->data; 252 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 253 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 254 255 PetscFunctionBegin; 256 pc_gamg_agg->use_low_mem_filter = b; 257 PetscFunctionReturn(PETSC_SUCCESS); 258 } 259 260 static PetscErrorCode PCGAMGSetGraphSymmetrize_AGG(PC pc, PetscBool b) 261 { 262 PC_MG *mg = (PC_MG *)pc->data; 263 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 264 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 265 266 PetscFunctionBegin; 267 pc_gamg_agg->graph_symmetrize = b; 268 PetscFunctionReturn(PETSC_SUCCESS); 269 } 270 271 static PetscErrorCode PCGAMGMISkSetMinDegreeOrdering_AGG(PC pc, PetscBool b) 272 { 273 PC_MG *mg = (PC_MG *)pc->data; 274 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 275 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 276 277 PetscFunctionBegin; 278 pc_gamg_agg->use_minimum_degree_ordering = b; 279 PetscFunctionReturn(PETSC_SUCCESS); 280 } 281 282 static PetscErrorCode PCSetFromOptions_GAMG_AGG(PC pc, PetscOptionItems PetscOptionsObject) 283 { 284 PC_MG *mg = (PC_MG *)pc->data; 285 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 286 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 287 PetscBool n_aggressive_flg, old_sq_provided = PETSC_FALSE, new_sq_provided = PETSC_FALSE, new_sqr_graph = pc_gamg_agg->use_aggressive_square_graph; 288 PetscInt nsq_graph_old = 0; 289 290 PetscFunctionBegin; 291 PetscOptionsHeadBegin(PetscOptionsObject, "GAMG-AGG options"); 292 PetscCall(PetscOptionsInt("-pc_gamg_agg_nsmooths", "number of smoothing steps to construct prolongation, usually 1", "PCGAMGSetNSmooths", pc_gamg_agg->nsmooths, &pc_gamg_agg->nsmooths, NULL)); 293 // aggressive coarsening logic with deprecated -pc_gamg_square_graph 294 PetscCall(PetscOptionsInt("-pc_gamg_aggressive_coarsening", "Number of aggressive coarsening (MIS-2) levels from finest", "PCGAMGSetAggressiveLevels", pc_gamg_agg->aggressive_coarsening_levels, &pc_gamg_agg->aggressive_coarsening_levels, &n_aggressive_flg)); 295 if (!n_aggressive_flg) 296 PetscCall(PetscOptionsInt("-pc_gamg_square_graph", "Number of aggressive coarsening (MIS-2) levels from finest (deprecated alias for -pc_gamg_aggressive_coarsening)", "PCGAMGSetAggressiveLevels", nsq_graph_old, &nsq_graph_old, &old_sq_provided)); 297 PetscCall(PetscOptionsBool("-pc_gamg_aggressive_square_graph", "Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening", "PCGAMGSetAggressiveSquareGraph", new_sqr_graph, &pc_gamg_agg->use_aggressive_square_graph, &new_sq_provided)); 298 if (!new_sq_provided && old_sq_provided) { 299 pc_gamg_agg->aggressive_coarsening_levels = nsq_graph_old; // could be zero 300 pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE; 301 } 302 if (new_sq_provided && old_sq_provided) 303 PetscCall(PetscInfo(pc, "Warning: both -pc_gamg_square_graph and -pc_gamg_aggressive_coarsening are used. -pc_gamg_square_graph is deprecated, Number of aggressive levels is %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels)); 304 PetscCall(PetscOptionsBool("-pc_gamg_mis_k_minimum_degree_ordering", "Use minimum degree ordering for greedy MIS", "PCGAMGMISkSetMinDegreeOrdering", pc_gamg_agg->use_minimum_degree_ordering, &pc_gamg_agg->use_minimum_degree_ordering, NULL)); 305 PetscCall(PetscOptionsBool("-pc_gamg_low_memory_threshold_filter", "Use the (built-in) low memory graph/matrix filter", "PCGAMGSetLowMemoryFilter", pc_gamg_agg->use_low_mem_filter, &pc_gamg_agg->use_low_mem_filter, NULL)); 306 PetscCall(PetscOptionsInt("-pc_gamg_aggressive_mis_k", "Number of levels of multigrid to use.", "PCGAMGMISkSetAggressive", pc_gamg_agg->aggressive_mis_k, &pc_gamg_agg->aggressive_mis_k, NULL)); 307 PetscCall(PetscOptionsBool("-pc_gamg_graph_symmetrize", "Symmetrize graph for coarsening", "PCGAMGSetGraphSymmetrize", pc_gamg_agg->graph_symmetrize, &pc_gamg_agg->graph_symmetrize, NULL)); 308 PetscOptionsHeadEnd(); 309 PetscFunctionReturn(PETSC_SUCCESS); 310 } 311 312 static PetscErrorCode PCDestroy_GAMG_AGG(PC pc) 313 { 314 PC_MG *mg = (PC_MG *)pc->data; 315 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 316 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 317 318 PetscFunctionBegin; 319 PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs)); 320 PetscCall(PetscFree(pc_gamg->subctx)); 321 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", NULL)); 322 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", NULL)); 323 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", NULL)); 324 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", NULL)); 325 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", NULL)); 326 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", NULL)); 327 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", NULL)); 328 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL)); 329 PetscFunctionReturn(PETSC_SUCCESS); 330 } 331 332 /* 333 PCSetCoordinates_AGG 334 335 Collective 336 337 Input Parameter: 338 . pc - the preconditioner context 339 . ndm - dimension of data (used for dof/vertex for Stokes) 340 . a_nloc - number of vertices local 341 . coords - [a_nloc][ndm] - interleaved coordinate data: {x_0, y_0, z_0, x_1, y_1, ...} 342 */ 343 344 static PetscErrorCode PCSetCoordinates_AGG(PC pc, PetscInt ndm, PetscInt a_nloc, PetscReal *coords) 345 { 346 PC_MG *mg = (PC_MG *)pc->data; 347 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 348 PetscInt arrsz, kk, ii, jj, nloc, ndatarows, ndf; 349 Mat mat = pc->pmat; 350 351 PetscFunctionBegin; 352 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 353 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 354 nloc = a_nloc; 355 356 /* SA: null space vectors */ 357 PetscCall(MatGetBlockSize(mat, &ndf)); /* this does not work for Stokes */ 358 if (coords && ndf == 1) pc_gamg->data_cell_cols = 1; /* scalar w/ coords and SA (not needed) */ 359 else if (coords) { 360 PetscCheck(ndm <= ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "degrees of motion %" PetscInt_FMT " > block size %" PetscInt_FMT, ndm, ndf); 361 pc_gamg->data_cell_cols = (ndm == 2 ? 3 : 6); /* displacement elasticity */ 362 if (ndm != ndf) PetscCheck(pc_gamg->data_cell_cols == ndf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Don't know how to create null space for ndm=%" PetscInt_FMT ", ndf=%" PetscInt_FMT ". Use MatSetNearNullSpace().", ndm, ndf); 363 } else pc_gamg->data_cell_cols = ndf; /* no data, force SA with constant null space vectors */ 364 pc_gamg->data_cell_rows = ndatarows = ndf; 365 PetscCheck(pc_gamg->data_cell_cols > 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "pc_gamg->data_cell_cols %" PetscInt_FMT " <= 0", pc_gamg->data_cell_cols); 366 arrsz = nloc * pc_gamg->data_cell_rows * pc_gamg->data_cell_cols; 367 368 if (!pc_gamg->data || (pc_gamg->data_sz != arrsz)) { 369 PetscCall(PetscFree(pc_gamg->data)); 370 PetscCall(PetscMalloc1(arrsz + 1, &pc_gamg->data)); 371 } 372 /* copy data in - column-oriented */ 373 for (kk = 0; kk < nloc; kk++) { 374 const PetscInt M = nloc * pc_gamg->data_cell_rows; /* stride into data */ 375 PetscReal *data = &pc_gamg->data[kk * ndatarows]; /* start of cell */ 376 377 if (pc_gamg->data_cell_cols == 1) *data = 1.0; 378 else { 379 /* translational modes */ 380 for (ii = 0; ii < ndatarows; ii++) { 381 for (jj = 0; jj < ndatarows; jj++) { 382 if (ii == jj) data[ii * M + jj] = 1.0; 383 else data[ii * M + jj] = 0.0; 384 } 385 } 386 387 /* rotational modes */ 388 if (coords) { 389 if (ndm == 2) { 390 data += 2 * M; 391 data[0] = -coords[2 * kk + 1]; 392 data[1] = coords[2 * kk]; 393 } else { 394 data += 3 * M; 395 data[0] = 0.0; 396 data[M + 0] = coords[3 * kk + 2]; 397 data[2 * M + 0] = -coords[3 * kk + 1]; 398 data[1] = -coords[3 * kk + 2]; 399 data[M + 1] = 0.0; 400 data[2 * M + 1] = coords[3 * kk]; 401 data[2] = coords[3 * kk + 1]; 402 data[M + 2] = -coords[3 * kk]; 403 data[2 * M + 2] = 0.0; 404 } 405 } 406 } 407 } 408 pc_gamg->data_sz = arrsz; 409 PetscFunctionReturn(PETSC_SUCCESS); 410 } 411 412 /* 413 PCSetData_AGG - called if data is not set with PCSetCoordinates. 414 Looks in Mat for near null space. 415 Does not work for Stokes 416 417 Input Parameter: 418 . pc - 419 . a_A - matrix to get (near) null space out of. 420 */ 421 static PetscErrorCode PCSetData_AGG(PC pc, Mat a_A) 422 { 423 PC_MG *mg = (PC_MG *)pc->data; 424 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 425 MatNullSpace mnull; 426 427 PetscFunctionBegin; 428 PetscCall(MatGetNearNullSpace(a_A, &mnull)); 429 if (!mnull) { 430 DM dm; 431 432 PetscCall(PCGetDM(pc, &dm)); 433 if (!dm) PetscCall(MatGetDM(a_A, &dm)); 434 if (dm) { 435 PetscObject deformation; 436 PetscInt Nf; 437 438 PetscCall(DMGetNumFields(dm, &Nf)); 439 if (Nf) { 440 PetscCall(DMGetField(dm, 0, NULL, &deformation)); 441 if (deformation) { 442 PetscCall(PetscObjectQuery(deformation, "nearnullspace", (PetscObject *)&mnull)); 443 if (!mnull) PetscCall(PetscObjectQuery(deformation, "nullspace", (PetscObject *)&mnull)); 444 } 445 } 446 } 447 } 448 449 if (!mnull) { 450 PetscInt bs, NN, MM; 451 452 PetscCall(MatGetBlockSize(a_A, &bs)); 453 PetscCall(MatGetLocalSize(a_A, &MM, &NN)); 454 PetscCheck(MM % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MM %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, MM, bs); 455 PetscCall(PCSetCoordinates_AGG(pc, bs, MM / bs, NULL)); 456 } else { 457 PetscReal *nullvec; 458 PetscBool has_const; 459 PetscInt i, j, mlocal, nvec, bs; 460 const Vec *vecs; 461 const PetscScalar *v; 462 463 PetscCall(MatGetLocalSize(a_A, &mlocal, NULL)); 464 PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs)); 465 for (i = 0; i < nvec; i++) { 466 PetscCall(VecGetLocalSize(vecs[i], &j)); 467 PetscCheck(j == mlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Attached null space vector size %" PetscInt_FMT " != matrix size %" PetscInt_FMT, j, mlocal); 468 } 469 pc_gamg->data_sz = (nvec + !!has_const) * mlocal; 470 PetscCall(PetscMalloc1((nvec + !!has_const) * mlocal, &nullvec)); 471 if (has_const) 472 for (i = 0; i < mlocal; i++) nullvec[i] = 1.0; 473 for (i = 0; i < nvec; i++) { 474 PetscCall(VecGetArrayRead(vecs[i], &v)); 475 for (j = 0; j < mlocal; j++) nullvec[(i + !!has_const) * mlocal + j] = PetscRealPart(v[j]); 476 PetscCall(VecRestoreArrayRead(vecs[i], &v)); 477 } 478 pc_gamg->data = nullvec; 479 pc_gamg->data_cell_cols = (nvec + !!has_const); 480 PetscCall(MatGetBlockSize(a_A, &bs)); 481 pc_gamg->data_cell_rows = bs; 482 } 483 PetscFunctionReturn(PETSC_SUCCESS); 484 } 485 486 /* 487 formProl0 - collect null space data for each aggregate, do QR, put R in coarse grid data and Q in P_0 488 489 Input Parameter: 490 . agg_llists - list of arrays with aggregates -- list from selected vertices of aggregate unselected vertices 491 . bs - row block size 492 . nSAvec - column bs of new P 493 . my0crs - global index of start of locals 494 . data_stride - bs*(nloc nodes + ghost nodes) [data_stride][nSAvec] 495 . data_in[data_stride*nSAvec] - local data on fine grid 496 . flid_fgid[data_stride/bs] - make local to global IDs, includes ghosts in 'locals_llist' 497 498 Output Parameter: 499 . a_data_out - in with fine grid data (w/ghosts), out with coarse grid data 500 . a_Prol - prolongation operator 501 */ 502 static PetscErrorCode formProl0(PetscCoarsenData *agg_llists, PetscInt bs, PetscInt nSAvec, PetscInt my0crs, PetscInt data_stride, PetscReal data_in[], const PetscInt flid_fgid[], PetscReal **a_data_out, Mat a_Prol) 503 { 504 PetscInt Istart, my0, Iend, nloc, clid, flid = 0, aggID, kk, jj, ii, mm, nSelected, minsz, nghosts, out_data_stride; 505 MPI_Comm comm; 506 PetscReal *out_data; 507 PetscCDIntNd *pos; 508 PCGAMGHashTable fgid_flid; 509 510 PetscFunctionBegin; 511 PetscCall(PetscObjectGetComm((PetscObject)a_Prol, &comm)); 512 PetscCall(MatGetOwnershipRange(a_Prol, &Istart, &Iend)); 513 nloc = (Iend - Istart) / bs; 514 my0 = Istart / bs; 515 PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT " must be divisible by bs %" PetscInt_FMT, Iend, Istart, bs); 516 Iend /= bs; 517 nghosts = data_stride / bs - nloc; 518 519 PetscCall(PCGAMGHashTableCreate(2 * nghosts + 1, &fgid_flid)); 520 for (kk = 0; kk < nghosts; kk++) PetscCall(PCGAMGHashTableAdd(&fgid_flid, flid_fgid[nloc + kk], nloc + kk)); 521 522 /* count selected -- same as number of cols of P */ 523 for (nSelected = mm = 0; mm < nloc; mm++) { 524 PetscBool ise; 525 526 PetscCall(PetscCDIsEmptyAt(agg_llists, mm, &ise)); 527 if (!ise) nSelected++; 528 } 529 PetscCall(MatGetOwnershipRangeColumn(a_Prol, &ii, &jj)); 530 PetscCheck((ii / nSAvec) == my0crs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ii %" PetscInt_FMT " /nSAvec %" PetscInt_FMT " != my0crs %" PetscInt_FMT, ii, nSAvec, my0crs); 531 PetscCheck(nSelected == (jj - ii) / nSAvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "nSelected %" PetscInt_FMT " != (jj %" PetscInt_FMT " - ii %" PetscInt_FMT ")/nSAvec %" PetscInt_FMT, nSelected, jj, ii, nSAvec); 532 533 /* aloc space for coarse point data (output) */ 534 out_data_stride = nSelected * nSAvec; 535 536 PetscCall(PetscMalloc1(out_data_stride * nSAvec, &out_data)); 537 for (ii = 0; ii < out_data_stride * nSAvec; ii++) out_data[ii] = PETSC_MAX_REAL; 538 *a_data_out = out_data; /* output - stride nSelected*nSAvec */ 539 540 /* find points and set prolongation */ 541 minsz = 100; 542 for (mm = clid = 0; mm < nloc; mm++) { 543 PetscCall(PetscCDCountAt(agg_llists, mm, &jj)); 544 if (jj > 0) { 545 const PetscInt lid = mm, cgid = my0crs + clid; 546 PetscInt cids[100]; /* max bs */ 547 PetscBLASInt asz, M, N, INFO; 548 PetscBLASInt Mdata, LDA, LWORK; 549 PetscScalar *qqc, *qqr, *TAU, *WORK; 550 PetscInt *fids; 551 PetscReal *data; 552 553 PetscCall(PetscBLASIntCast(jj, &asz)); 554 PetscCall(PetscBLASIntCast(asz * bs, &M)); 555 PetscCall(PetscBLASIntCast(nSAvec, &N)); 556 PetscCall(PetscBLASIntCast(M + ((N - M > 0) ? N - M : 0), &Mdata)); 557 PetscCall(PetscBLASIntCast(Mdata, &LDA)); 558 PetscCall(PetscBLASIntCast(N * bs, &LWORK)); 559 /* count agg */ 560 if (asz < minsz) minsz = asz; 561 562 /* get block */ 563 PetscCall(PetscMalloc5(Mdata * N, &qqc, M * N, &qqr, N, &TAU, LWORK, &WORK, M, &fids)); 564 565 aggID = 0; 566 PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos)); 567 while (pos) { 568 PetscInt gid1; 569 570 PetscCall(PetscCDIntNdGetID(pos, &gid1)); 571 PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos)); 572 573 if (gid1 >= my0 && gid1 < Iend) flid = gid1 - my0; 574 else { 575 PetscCall(PCGAMGHashTableFind(&fgid_flid, gid1, &flid)); 576 PetscCheck(flid >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find gid1 in table"); 577 } 578 /* copy in B_i matrix - column-oriented */ 579 data = &data_in[flid * bs]; 580 for (ii = 0; ii < bs; ii++) { 581 for (jj = 0; jj < N; jj++) { 582 PetscReal d = data[jj * data_stride + ii]; 583 584 qqc[jj * Mdata + aggID * bs + ii] = d; 585 } 586 } 587 /* set fine IDs */ 588 for (kk = 0; kk < bs; kk++) fids[aggID * bs + kk] = flid_fgid[flid] * bs + kk; 589 aggID++; 590 } 591 592 /* pad with zeros */ 593 for (ii = asz * bs; ii < Mdata; ii++) { 594 for (jj = 0; jj < N; jj++, kk++) qqc[jj * Mdata + ii] = .0; 595 } 596 597 /* QR */ 598 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 599 PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&Mdata, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO)); 600 PetscCall(PetscFPTrapPop()); 601 PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xGEQRF error"); 602 /* get R - column-oriented - output B_{i+1} */ 603 { 604 PetscReal *data = &out_data[clid * nSAvec]; 605 606 for (jj = 0; jj < nSAvec; jj++) { 607 for (ii = 0; ii < nSAvec; ii++) { 608 PetscCheck(data[jj * out_data_stride + ii] == PETSC_MAX_REAL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "data[jj*out_data_stride + ii] != %e", (double)PETSC_MAX_REAL); 609 if (ii <= jj) data[jj * out_data_stride + ii] = PetscRealPart(qqc[jj * Mdata + ii]); 610 else data[jj * out_data_stride + ii] = 0.; 611 } 612 } 613 } 614 615 /* get Q - row-oriented */ 616 PetscCallBLAS("LAPACKorgqr", LAPACKorgqr_(&Mdata, &N, &N, qqc, &LDA, TAU, WORK, &LWORK, &INFO)); 617 PetscCheck(INFO == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xORGQR error arg %" PetscBLASInt_FMT, -INFO); 618 619 for (ii = 0; ii < M; ii++) { 620 for (jj = 0; jj < N; jj++) qqr[N * ii + jj] = qqc[jj * Mdata + ii]; 621 } 622 623 /* add diagonal block of P0 */ 624 for (kk = 0; kk < N; kk++) { cids[kk] = N * cgid + kk; /* global col IDs in P0 */ } 625 PetscCall(MatSetValues(a_Prol, M, fids, N, cids, qqr, INSERT_VALUES)); 626 PetscCall(PetscFree5(qqc, qqr, TAU, WORK, fids)); 627 clid++; 628 } /* coarse agg */ 629 } /* for all fine nodes */ 630 PetscCall(MatAssemblyBegin(a_Prol, MAT_FINAL_ASSEMBLY)); 631 PetscCall(MatAssemblyEnd(a_Prol, MAT_FINAL_ASSEMBLY)); 632 PetscCall(PCGAMGHashTableDestroy(&fgid_flid)); 633 PetscFunctionReturn(PETSC_SUCCESS); 634 } 635 636 static PetscErrorCode PCView_GAMG_AGG(PC pc, PetscViewer viewer) 637 { 638 PC_MG *mg = (PC_MG *)pc->data; 639 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 640 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 641 642 PetscFunctionBegin; 643 PetscCall(PetscViewerASCIIPrintf(viewer, " AGG specific options\n")); 644 PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %" PetscInt_FMT "\n", pc_gamg_agg->aggressive_coarsening_levels)); 645 if (pc_gamg_agg->aggressive_coarsening_levels > 0) { 646 PetscCall(PetscViewerASCIIPrintf(viewer, " %s aggressive coarsening\n", !pc_gamg_agg->use_aggressive_square_graph ? "MIS-k" : "Square graph")); 647 if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(PetscViewerASCIIPrintf(viewer, " MIS-%" PetscInt_FMT " coarsening on aggressive levels\n", pc_gamg_agg->aggressive_mis_k)); 648 } 649 PetscCall(PetscViewerASCIIPushTab(viewer)); 650 PetscCall(PetscViewerASCIIPushTab(viewer)); 651 PetscCall(PetscViewerASCIIPushTab(viewer)); 652 PetscCall(PetscViewerASCIIPushTab(viewer)); 653 if (pc_gamg_agg->crs) PetscCall(MatCoarsenView(pc_gamg_agg->crs, viewer)); 654 else PetscCall(PetscViewerASCIIPrintf(viewer, "Coarsening algorithm not yet selected\n")); 655 PetscCall(PetscViewerASCIIPopTab(viewer)); 656 PetscCall(PetscViewerASCIIPopTab(viewer)); 657 PetscCall(PetscViewerASCIIPopTab(viewer)); 658 PetscCall(PetscViewerASCIIPopTab(viewer)); 659 PetscCall(PetscViewerASCIIPrintf(viewer, " Number smoothing steps to construct prolongation %" PetscInt_FMT "\n", pc_gamg_agg->nsmooths)); 660 PetscFunctionReturn(PETSC_SUCCESS); 661 } 662 663 static PetscErrorCode PCGAMGCreateGraph_AGG(PC pc, Mat Amat, Mat *a_Gmat) 664 { 665 PC_MG *mg = (PC_MG *)pc->data; 666 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 667 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 668 const PetscReal vfilter = pc_gamg->threshold[pc_gamg->current_level]; 669 PetscBool ishem, ismis; 670 const char *prefix; 671 MatInfo info0, info1; 672 PetscInt bs; 673 674 PetscFunctionBegin; 675 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0)); 676 /* Note: depending on the algorithm that will be used for computing the coarse grid points this should pass PETSC_TRUE or PETSC_FALSE as the first argument */ 677 /* MATCOARSENHEM requires numerical weights for edges so ensure they are computed */ 678 PetscCall(MatCoarsenDestroy(&pc_gamg_agg->crs)); 679 PetscCall(MatCoarsenCreate(PetscObjectComm((PetscObject)pc), &pc_gamg_agg->crs)); 680 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc, &prefix)); 681 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix)); 682 PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)pc_gamg_agg->crs, "pc_gamg_")); 683 PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); 684 PetscCall(MatGetBlockSize(Amat, &bs)); 685 // check for valid indices wrt bs 686 for (int ii = 0; ii < pc_gamg_agg->crs->strength_index_size; ii++) { 687 PetscCheck(pc_gamg_agg->crs->strength_index[ii] >= 0 && pc_gamg_agg->crs->strength_index[ii] < bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Indices (%" PetscInt_FMT ") must be non-negative and < block size (%" PetscInt_FMT "), NB, can not use -mat_coarsen_strength_index with -mat_coarsen_strength_index", 688 pc_gamg_agg->crs->strength_index[ii], bs); 689 } 690 PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENHEM, &ishem)); 691 if (ishem) { 692 if (pc_gamg_agg->aggressive_coarsening_levels) PetscCall(PetscInfo(pc, "HEM and aggressive coarsening ignored: HEM using %" PetscInt_FMT " iterations\n", pc_gamg_agg->crs->max_it)); 693 pc_gamg_agg->aggressive_coarsening_levels = 0; // aggressive and HEM does not make sense 694 PetscCall(MatCoarsenSetMaximumIterations(pc_gamg_agg->crs, pc_gamg_agg->crs->max_it)); // for code coverage 695 PetscCall(MatCoarsenSetThreshold(pc_gamg_agg->crs, vfilter)); // for code coverage 696 } else { 697 PetscCall(PetscObjectTypeCompare((PetscObject)pc_gamg_agg->crs, MATCOARSENMIS, &ismis)); 698 if (ismis && pc_gamg_agg->aggressive_coarsening_levels && !pc_gamg_agg->use_aggressive_square_graph) { 699 PetscCall(PetscInfo(pc, "MIS and aggressive coarsening and no square graph: force square graph\n")); 700 pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE; 701 } 702 } 703 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0)); 704 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0)); 705 PetscCall(MatGetInfo(Amat, MAT_LOCAL, &info0)); /* global reduction */ 706 707 if (ishem || pc_gamg_agg->use_low_mem_filter) { 708 PetscCall(MatCreateGraph(Amat, pc_gamg_agg->graph_symmetrize, (vfilter >= 0 || ishem) ? PETSC_TRUE : PETSC_FALSE, vfilter, pc_gamg_agg->crs->strength_index_size, pc_gamg_agg->crs->strength_index, a_Gmat)); 709 } else { 710 // make scalar graph, symmetrize if not known to be symmetric, scale, but do not filter (expensive) 711 PetscCall(MatCreateGraph(Amat, pc_gamg_agg->graph_symmetrize, PETSC_TRUE, -1, pc_gamg_agg->crs->strength_index_size, pc_gamg_agg->crs->strength_index, a_Gmat)); 712 if (vfilter >= 0) { 713 PetscInt Istart, Iend, ncols, nnz0, nnz1, NN, MM, nloc; 714 Mat tGmat, Gmat = *a_Gmat; 715 MPI_Comm comm; 716 const PetscScalar *vals; 717 const PetscInt *idx; 718 PetscInt *d_nnz, *o_nnz, kk, *garray = NULL, *AJ, maxcols = 0; 719 MatScalar *AA; // this is checked in graph 720 PetscBool isseqaij; 721 Mat a, b, c; 722 MatType jtype; 723 724 PetscCall(PetscObjectGetComm((PetscObject)Gmat, &comm)); 725 PetscCall(PetscObjectBaseTypeCompare((PetscObject)Gmat, MATSEQAIJ, &isseqaij)); 726 PetscCall(MatGetType(Gmat, &jtype)); 727 PetscCall(MatCreate(comm, &tGmat)); 728 PetscCall(MatSetType(tGmat, jtype)); 729 730 /* TODO GPU: this can be called when filter = 0 -> Probably provide MatAIJThresholdCompress that compresses the entries below a threshold? 731 Also, if the matrix is symmetric, can we skip this 732 operation? It can be very expensive on large matrices. */ 733 734 // global sizes 735 PetscCall(MatGetSize(Gmat, &MM, &NN)); 736 PetscCall(MatGetOwnershipRange(Gmat, &Istart, &Iend)); 737 nloc = Iend - Istart; 738 PetscCall(PetscMalloc2(nloc, &d_nnz, nloc, &o_nnz)); 739 if (isseqaij) { 740 a = Gmat; 741 b = NULL; 742 } else { 743 Mat_MPIAIJ *d = (Mat_MPIAIJ *)Gmat->data; 744 745 a = d->A; 746 b = d->B; 747 garray = d->garray; 748 } 749 /* Determine upper bound on non-zeros needed in new filtered matrix */ 750 for (PetscInt row = 0; row < nloc; row++) { 751 PetscCall(MatGetRow(a, row, &ncols, NULL, NULL)); 752 d_nnz[row] = ncols; 753 if (ncols > maxcols) maxcols = ncols; 754 PetscCall(MatRestoreRow(a, row, &ncols, NULL, NULL)); 755 } 756 if (b) { 757 for (PetscInt row = 0; row < nloc; row++) { 758 PetscCall(MatGetRow(b, row, &ncols, NULL, NULL)); 759 o_nnz[row] = ncols; 760 if (ncols > maxcols) maxcols = ncols; 761 PetscCall(MatRestoreRow(b, row, &ncols, NULL, NULL)); 762 } 763 } 764 PetscCall(MatSetSizes(tGmat, nloc, nloc, MM, MM)); 765 PetscCall(MatSetBlockSizes(tGmat, 1, 1)); 766 PetscCall(MatSeqAIJSetPreallocation(tGmat, 0, d_nnz)); 767 PetscCall(MatMPIAIJSetPreallocation(tGmat, 0, d_nnz, 0, o_nnz)); 768 PetscCall(MatSetOption(tGmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 769 PetscCall(PetscFree2(d_nnz, o_nnz)); 770 PetscCall(PetscMalloc2(maxcols, &AA, maxcols, &AJ)); 771 nnz0 = nnz1 = 0; 772 for (c = a, kk = 0; c && kk < 2; c = b, kk++) { 773 for (PetscInt row = 0, grow = Istart, ncol_row, jj; row < nloc; row++, grow++) { 774 PetscCall(MatGetRow(c, row, &ncols, &idx, &vals)); 775 for (ncol_row = jj = 0; jj < ncols; jj++, nnz0++) { 776 PetscScalar sv = PetscAbs(PetscRealPart(vals[jj])); 777 if (PetscRealPart(sv) > vfilter) { 778 PetscInt cid = idx[jj] + Istart; //diag 779 780 nnz1++; 781 if (c != a) cid = garray[idx[jj]]; 782 AA[ncol_row] = vals[jj]; 783 AJ[ncol_row] = cid; 784 ncol_row++; 785 } 786 } 787 PetscCall(MatRestoreRow(c, row, &ncols, &idx, &vals)); 788 PetscCall(MatSetValues(tGmat, 1, &grow, ncol_row, AJ, AA, INSERT_VALUES)); 789 } 790 } 791 PetscCall(PetscFree2(AA, AJ)); 792 PetscCall(MatAssemblyBegin(tGmat, MAT_FINAL_ASSEMBLY)); 793 PetscCall(MatAssemblyEnd(tGmat, MAT_FINAL_ASSEMBLY)); 794 PetscCall(MatPropagateSymmetryOptions(Gmat, tGmat)); /* Normal Mat options are not relevant ? */ 795 PetscCall(PetscInfo(pc, "\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%" PetscInt_FMT ", max row size %" PetscInt_FMT "\n", (!nnz0) ? 1. : 100. * (double)nnz1 / (double)nnz0, (double)vfilter, (!nloc) ? 1. : (double)nnz0 / (double)nloc, MM, maxcols)); 796 PetscCall(MatViewFromOptions(tGmat, NULL, "-mat_filter_graph_view")); 797 PetscCall(MatDestroy(&Gmat)); 798 *a_Gmat = tGmat; 799 } 800 } 801 802 PetscCall(MatGetInfo(*a_Gmat, MAT_LOCAL, &info1)); /* global reduction */ 803 if (info0.nz_used > 0) PetscCall(PetscInfo(pc, "Filtering left %g %% edges in graph (%e %e)\n", 100.0 * info1.nz_used * (double)(bs * bs) / info0.nz_used, info0.nz_used, info1.nz_used)); 804 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_GRAPH], 0, 0, 0, 0)); 805 PetscFunctionReturn(PETSC_SUCCESS); 806 } 807 808 typedef PetscInt NState; 809 static const NState NOT_DONE = -2; 810 static const NState DELETED = -1; 811 static const NState REMOVED = -3; 812 #define IS_SELECTED(s) (s != DELETED && s != NOT_DONE && s != REMOVED) 813 814 /* 815 fixAggregatesWithSquare - greedy grab of with G1 (unsquared graph) -- AIJ specific -- change to fixAggregatesWithSquare -- TODD 816 - AGG-MG specific: clears singletons out of 'selected_2' 817 818 Input Parameter: 819 . Gmat_2 - global matrix of squared graph (data not defined) 820 . Gmat_1 - base graph to grab with base graph 821 Input/Output Parameter: 822 . aggs_2 - linked list of aggs with gids) 823 */ 824 static PetscErrorCode fixAggregatesWithSquare(PC pc, Mat Gmat_2, Mat Gmat_1, PetscCoarsenData *aggs_2) 825 { 826 PetscBool isMPI; 827 Mat_SeqAIJ *matA_1, *matB_1 = NULL; 828 MPI_Comm comm; 829 PetscInt lid, *ii, *idx, ix, Iend, my0, kk, n, j; 830 Mat_MPIAIJ *mpimat_2 = NULL, *mpimat_1 = NULL; 831 const PetscInt nloc = Gmat_2->rmap->n; 832 PetscScalar *cpcol_1_state, *cpcol_2_state, *cpcol_2_par_orig, *lid_parent_gid; 833 PetscInt *lid_cprowID_1 = NULL; 834 NState *lid_state; 835 Vec ghost_par_orig2; 836 PetscMPIInt rank; 837 838 PetscFunctionBegin; 839 PetscCall(PetscObjectGetComm((PetscObject)Gmat_2, &comm)); 840 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 841 PetscCall(MatGetOwnershipRange(Gmat_1, &my0, &Iend)); 842 843 /* get submatrices */ 844 PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATMPIAIJ, &isMPI)); 845 PetscCall(PetscInfo(pc, "isMPI = %s\n", isMPI ? "yes" : "no")); 846 PetscCall(PetscMalloc3(nloc, &lid_state, nloc, &lid_parent_gid, nloc, &lid_cprowID_1)); 847 for (lid = 0; lid < nloc; lid++) lid_cprowID_1[lid] = -1; 848 if (isMPI) { 849 /* grab matrix objects */ 850 mpimat_2 = (Mat_MPIAIJ *)Gmat_2->data; 851 mpimat_1 = (Mat_MPIAIJ *)Gmat_1->data; 852 matA_1 = (Mat_SeqAIJ *)mpimat_1->A->data; 853 matB_1 = (Mat_SeqAIJ *)mpimat_1->B->data; 854 855 /* force compressed row storage for B matrix in AuxMat */ 856 PetscCall(MatCheckCompressedRow(mpimat_1->B, matB_1->nonzerorowcnt, &matB_1->compressedrow, matB_1->i, Gmat_1->rmap->n, -1.0)); 857 for (ix = 0; ix < matB_1->compressedrow.nrows; ix++) { 858 PetscInt lid = matB_1->compressedrow.rindex[ix]; 859 860 PetscCheck(lid <= nloc && lid >= -1, PETSC_COMM_SELF, PETSC_ERR_USER, "lid %" PetscInt_FMT " out of range. nloc = %" PetscInt_FMT, lid, nloc); 861 if (lid != -1) lid_cprowID_1[lid] = ix; 862 } 863 } else { 864 PetscBool isAIJ; 865 866 PetscCall(PetscStrbeginswith(((PetscObject)Gmat_1)->type_name, MATSEQAIJ, &isAIJ)); 867 PetscCheck(isAIJ, PETSC_COMM_SELF, PETSC_ERR_USER, "Require AIJ matrix."); 868 matA_1 = (Mat_SeqAIJ *)Gmat_1->data; 869 } 870 if (nloc > 0) { PetscCheck(!matB_1 || matB_1->compressedrow.use, PETSC_COMM_SELF, PETSC_ERR_PLIB, "matB_1 && !matB_1->compressedrow.use: PETSc bug???"); } 871 /* get state of locals and selected gid for deleted */ 872 for (lid = 0; lid < nloc; lid++) { 873 lid_parent_gid[lid] = -1.0; 874 lid_state[lid] = DELETED; 875 } 876 877 /* set lid_state */ 878 for (lid = 0; lid < nloc; lid++) { 879 PetscCDIntNd *pos; 880 881 PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos)); 882 if (pos) { 883 PetscInt gid1; 884 885 PetscCall(PetscCDIntNdGetID(pos, &gid1)); 886 PetscCheck(gid1 == lid + my0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "gid1 %" PetscInt_FMT " != lid %" PetscInt_FMT " + my0 %" PetscInt_FMT, gid1, lid, my0); 887 lid_state[lid] = gid1; 888 } 889 } 890 891 /* map local to selected local, DELETED means a ghost owns it */ 892 for (lid = 0; lid < nloc; lid++) { 893 NState state = lid_state[lid]; 894 895 if (IS_SELECTED(state)) { 896 PetscCDIntNd *pos; 897 898 PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos)); 899 while (pos) { 900 PetscInt gid1; 901 902 PetscCall(PetscCDIntNdGetID(pos, &gid1)); 903 PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos)); 904 if (gid1 >= my0 && gid1 < Iend) lid_parent_gid[gid1 - my0] = (PetscScalar)(lid + my0); 905 } 906 } 907 } 908 /* get 'cpcol_1/2_state' & cpcol_2_par_orig - uses mpimat_1/2->lvec for temp space */ 909 if (isMPI) { 910 Vec tempVec; 911 912 /* get 'cpcol_1_state' */ 913 PetscCall(MatCreateVecs(Gmat_1, &tempVec, NULL)); 914 for (kk = 0, j = my0; kk < nloc; kk++, j++) { 915 PetscScalar v = (PetscScalar)lid_state[kk]; 916 917 PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES)); 918 } 919 PetscCall(VecAssemblyBegin(tempVec)); 920 PetscCall(VecAssemblyEnd(tempVec)); 921 PetscCall(VecScatterBegin(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD)); 922 PetscCall(VecScatterEnd(mpimat_1->Mvctx, tempVec, mpimat_1->lvec, INSERT_VALUES, SCATTER_FORWARD)); 923 PetscCall(VecGetArray(mpimat_1->lvec, &cpcol_1_state)); 924 /* get 'cpcol_2_state' */ 925 PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD)); 926 PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, mpimat_2->lvec, INSERT_VALUES, SCATTER_FORWARD)); 927 PetscCall(VecGetArray(mpimat_2->lvec, &cpcol_2_state)); 928 /* get 'cpcol_2_par_orig' */ 929 for (kk = 0, j = my0; kk < nloc; kk++, j++) { 930 PetscScalar v = lid_parent_gid[kk]; 931 932 PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES)); 933 } 934 PetscCall(VecAssemblyBegin(tempVec)); 935 PetscCall(VecAssemblyEnd(tempVec)); 936 PetscCall(VecDuplicate(mpimat_2->lvec, &ghost_par_orig2)); 937 PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD)); 938 PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghost_par_orig2, INSERT_VALUES, SCATTER_FORWARD)); 939 PetscCall(VecGetArray(ghost_par_orig2, &cpcol_2_par_orig)); 940 941 PetscCall(VecDestroy(&tempVec)); 942 } /* ismpi */ 943 for (lid = 0; lid < nloc; lid++) { 944 NState state = lid_state[lid]; 945 946 if (IS_SELECTED(state)) { 947 /* steal locals */ 948 ii = matA_1->i; 949 n = ii[lid + 1] - ii[lid]; 950 idx = matA_1->j + ii[lid]; 951 for (j = 0; j < n; j++) { 952 PetscInt lidj = idx[j], sgid; 953 NState statej = lid_state[lidj]; 954 955 if (statej == DELETED && (sgid = (PetscInt)PetscRealPart(lid_parent_gid[lidj])) != lid + my0) { /* steal local */ 956 lid_parent_gid[lidj] = (PetscScalar)(lid + my0); /* send this if sgid is not local */ 957 if (sgid >= my0 && sgid < Iend) { /* I'm stealing this local from a local sgid */ 958 PetscInt hav = 0, slid = sgid - my0, gidj = lidj + my0; 959 PetscCDIntNd *pos, *last = NULL; 960 961 /* looking for local from local so id_llist_2 works */ 962 PetscCall(PetscCDGetHeadPos(aggs_2, slid, &pos)); 963 while (pos) { 964 PetscInt gid; 965 966 PetscCall(PetscCDIntNdGetID(pos, &gid)); 967 if (gid == gidj) { 968 PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null"); 969 PetscCall(PetscCDRemoveNextNode(aggs_2, slid, last)); 970 PetscCall(PetscCDAppendNode(aggs_2, lid, pos)); 971 hav = 1; 972 break; 973 } else last = pos; 974 PetscCall(PetscCDGetNextPos(aggs_2, slid, &pos)); 975 } 976 if (hav != 1) { 977 PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find adj in 'selected' lists - structurally unsymmetric matrix"); 978 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav); 979 } 980 } else { /* I'm stealing this local, owned by a ghost */ 981 PetscCheck(sgid == -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mat has an un-symmetric graph. Use '-%spc_gamg_sym_graph true' to symmetrize the graph or '-%spc_gamg_threshold -1' if the matrix is structurally symmetric.", 982 ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : ""); 983 PetscCall(PetscCDAppendID(aggs_2, lid, lidj + my0)); 984 } 985 } 986 } /* local neighbors */ 987 } else if (state == DELETED /* && lid_cprowID_1 */) { 988 PetscInt sgidold = (PetscInt)PetscRealPart(lid_parent_gid[lid]); 989 990 /* see if I have a selected ghost neighbor that will steal me */ 991 if ((ix = lid_cprowID_1[lid]) != -1) { 992 ii = matB_1->compressedrow.i; 993 n = ii[ix + 1] - ii[ix]; 994 idx = matB_1->j + ii[ix]; 995 for (j = 0; j < n; j++) { 996 PetscInt cpid = idx[j]; 997 NState statej = (NState)PetscRealPart(cpcol_1_state[cpid]); 998 999 if (IS_SELECTED(statej) && sgidold != statej) { /* ghost will steal this, remove from my list */ 1000 lid_parent_gid[lid] = (PetscScalar)statej; /* send who selected */ 1001 if (sgidold >= my0 && sgidold < Iend) { /* this was mine */ 1002 PetscInt hav = 0, oldslidj = sgidold - my0; 1003 PetscCDIntNd *pos, *last = NULL; 1004 1005 /* remove from 'oldslidj' list */ 1006 PetscCall(PetscCDGetHeadPos(aggs_2, oldslidj, &pos)); 1007 while (pos) { 1008 PetscInt gid; 1009 1010 PetscCall(PetscCDIntNdGetID(pos, &gid)); 1011 if (lid + my0 == gid) { 1012 /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove lid from oldslidj list *\/ */ 1013 PetscCheck(last, PETSC_COMM_SELF, PETSC_ERR_PLIB, "last cannot be null"); 1014 PetscCall(PetscCDRemoveNextNode(aggs_2, oldslidj, last)); 1015 /* ghost (PetscScalar)statej will add this later */ 1016 hav = 1; 1017 break; 1018 } else last = pos; 1019 PetscCall(PetscCDGetNextPos(aggs_2, oldslidj, &pos)); 1020 } 1021 if (hav != 1) { 1022 PetscCheck(hav, PETSC_COMM_SELF, PETSC_ERR_PLIB, "failed to find (hav=%" PetscInt_FMT ") adj in 'selected' lists - structurally unsymmetric matrix", hav); 1023 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "found node %" PetscInt_FMT " times???", hav); 1024 } 1025 } else { 1026 /* TODO: ghosts remove this later */ 1027 } 1028 } 1029 } 1030 } 1031 } /* selected/deleted */ 1032 } /* node loop */ 1033 1034 if (isMPI) { 1035 PetscScalar *cpcol_2_parent, *cpcol_2_gid; 1036 Vec tempVec, ghostgids2, ghostparents2; 1037 PetscInt cpid, nghost_2; 1038 PCGAMGHashTable gid_cpid; 1039 1040 PetscCall(VecGetSize(mpimat_2->lvec, &nghost_2)); 1041 PetscCall(MatCreateVecs(Gmat_2, &tempVec, NULL)); 1042 1043 /* get 'cpcol_2_parent' */ 1044 for (kk = 0, j = my0; kk < nloc; kk++, j++) { PetscCall(VecSetValues(tempVec, 1, &j, &lid_parent_gid[kk], INSERT_VALUES)); } 1045 PetscCall(VecAssemblyBegin(tempVec)); 1046 PetscCall(VecAssemblyEnd(tempVec)); 1047 PetscCall(VecDuplicate(mpimat_2->lvec, &ghostparents2)); 1048 PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD)); 1049 PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostparents2, INSERT_VALUES, SCATTER_FORWARD)); 1050 PetscCall(VecGetArray(ghostparents2, &cpcol_2_parent)); 1051 1052 /* get 'cpcol_2_gid' */ 1053 for (kk = 0, j = my0; kk < nloc; kk++, j++) { 1054 PetscScalar v = (PetscScalar)j; 1055 1056 PetscCall(VecSetValues(tempVec, 1, &j, &v, INSERT_VALUES)); 1057 } 1058 PetscCall(VecAssemblyBegin(tempVec)); 1059 PetscCall(VecAssemblyEnd(tempVec)); 1060 PetscCall(VecDuplicate(mpimat_2->lvec, &ghostgids2)); 1061 PetscCall(VecScatterBegin(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD)); 1062 PetscCall(VecScatterEnd(mpimat_2->Mvctx, tempVec, ghostgids2, INSERT_VALUES, SCATTER_FORWARD)); 1063 PetscCall(VecGetArray(ghostgids2, &cpcol_2_gid)); 1064 PetscCall(VecDestroy(&tempVec)); 1065 1066 /* look for deleted ghosts and add to table */ 1067 PetscCall(PCGAMGHashTableCreate(2 * nghost_2 + 1, &gid_cpid)); 1068 for (cpid = 0; cpid < nghost_2; cpid++) { 1069 NState state = (NState)PetscRealPart(cpcol_2_state[cpid]); 1070 1071 if (state == DELETED) { 1072 PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]); 1073 PetscInt sgid_old = (PetscInt)PetscRealPart(cpcol_2_par_orig[cpid]); 1074 1075 if (sgid_old == -1 && sgid_new != -1) { 1076 PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]); 1077 1078 PetscCall(PCGAMGHashTableAdd(&gid_cpid, gid, cpid)); 1079 } 1080 } 1081 } 1082 1083 /* look for deleted ghosts and see if they moved - remove it */ 1084 for (lid = 0; lid < nloc; lid++) { 1085 NState state = lid_state[lid]; 1086 1087 if (IS_SELECTED(state)) { 1088 PetscCDIntNd *pos, *last = NULL; 1089 1090 /* look for deleted ghosts and see if they moved */ 1091 PetscCall(PetscCDGetHeadPos(aggs_2, lid, &pos)); 1092 while (pos) { 1093 PetscInt gid; 1094 1095 PetscCall(PetscCDIntNdGetID(pos, &gid)); 1096 if (gid < my0 || gid >= Iend) { 1097 PetscCall(PCGAMGHashTableFind(&gid_cpid, gid, &cpid)); 1098 if (cpid != -1) { 1099 /* a moved ghost - */ 1100 /* id_llist_2[lastid] = id_llist_2[flid]; /\* remove 'flid' from list *\/ */ 1101 PetscCall(PetscCDRemoveNextNode(aggs_2, lid, last)); 1102 } else last = pos; 1103 } else last = pos; 1104 1105 PetscCall(PetscCDGetNextPos(aggs_2, lid, &pos)); 1106 } /* loop over list of deleted */ 1107 } /* selected */ 1108 } 1109 PetscCall(PCGAMGHashTableDestroy(&gid_cpid)); 1110 1111 /* look at ghosts, see if they changed - and it */ 1112 for (cpid = 0; cpid < nghost_2; cpid++) { 1113 PetscInt sgid_new = (PetscInt)PetscRealPart(cpcol_2_parent[cpid]); 1114 1115 if (sgid_new >= my0 && sgid_new < Iend) { /* this is mine */ 1116 PetscInt gid = (PetscInt)PetscRealPart(cpcol_2_gid[cpid]); 1117 PetscInt slid_new = sgid_new - my0, hav = 0; 1118 PetscCDIntNd *pos; 1119 1120 /* search for this gid to see if I have it */ 1121 PetscCall(PetscCDGetHeadPos(aggs_2, slid_new, &pos)); 1122 while (pos) { 1123 PetscInt gidj; 1124 1125 PetscCall(PetscCDIntNdGetID(pos, &gidj)); 1126 PetscCall(PetscCDGetNextPos(aggs_2, slid_new, &pos)); 1127 1128 if (gidj == gid) { 1129 hav = 1; 1130 break; 1131 } 1132 } 1133 if (hav != 1) { 1134 /* insert 'flidj' into head of llist */ 1135 PetscCall(PetscCDAppendID(aggs_2, slid_new, gid)); 1136 } 1137 } 1138 } 1139 PetscCall(VecRestoreArray(mpimat_1->lvec, &cpcol_1_state)); 1140 PetscCall(VecRestoreArray(mpimat_2->lvec, &cpcol_2_state)); 1141 PetscCall(VecRestoreArray(ghostparents2, &cpcol_2_parent)); 1142 PetscCall(VecRestoreArray(ghostgids2, &cpcol_2_gid)); 1143 PetscCall(VecDestroy(&ghostgids2)); 1144 PetscCall(VecDestroy(&ghostparents2)); 1145 PetscCall(VecDestroy(&ghost_par_orig2)); 1146 } 1147 PetscCall(PetscFree3(lid_state, lid_parent_gid, lid_cprowID_1)); 1148 PetscFunctionReturn(PETSC_SUCCESS); 1149 } 1150 1151 /* 1152 PCGAMGCoarsen_AGG - supports squaring the graph (deprecated) and new graph for 1153 communication of QR data used with HEM and MISk coarsening 1154 1155 Input Parameter: 1156 . a_pc - this 1157 1158 Input/Output Parameter: 1159 . a_Gmat1 - graph to coarsen (in), graph off processor edges for QR gather scatter (out) 1160 1161 Output Parameter: 1162 . agg_lists - list of aggregates 1163 1164 */ 1165 static PetscErrorCode PCGAMGCoarsen_AGG(PC a_pc, Mat *a_Gmat1, PetscCoarsenData **agg_lists) 1166 { 1167 PC_MG *mg = (PC_MG *)a_pc->data; 1168 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1169 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 1170 Mat Gmat2, Gmat1 = *a_Gmat1; /* aggressive graph */ 1171 IS perm; 1172 PetscInt Istart, Iend, Ii, nloc, bs, nn; 1173 PetscInt *permute, *degree; 1174 PetscBool *bIndexSet; 1175 PetscReal hashfact; 1176 PetscInt iSwapIndex; 1177 PetscRandom random; 1178 MPI_Comm comm; 1179 1180 PetscFunctionBegin; 1181 PetscCall(PetscObjectGetComm((PetscObject)Gmat1, &comm)); 1182 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0)); 1183 PetscCall(MatGetLocalSize(Gmat1, &nn, NULL)); 1184 PetscCall(MatGetBlockSize(Gmat1, &bs)); 1185 PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bs %" PetscInt_FMT " must be 1", bs); 1186 nloc = nn / bs; 1187 /* get MIS aggs - randomize */ 1188 PetscCall(PetscMalloc2(nloc, &permute, nloc, °ree)); 1189 PetscCall(PetscCalloc1(nloc, &bIndexSet)); 1190 for (Ii = 0; Ii < nloc; Ii++) permute[Ii] = Ii; 1191 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &random)); 1192 PetscCall(MatGetOwnershipRange(Gmat1, &Istart, &Iend)); 1193 for (Ii = 0; Ii < nloc; Ii++) { 1194 PetscInt nc; 1195 1196 PetscCall(MatGetRow(Gmat1, Istart + Ii, &nc, NULL, NULL)); 1197 degree[Ii] = nc; 1198 PetscCall(MatRestoreRow(Gmat1, Istart + Ii, &nc, NULL, NULL)); 1199 } 1200 for (Ii = 0; Ii < nloc; Ii++) { 1201 PetscCall(PetscRandomGetValueReal(random, &hashfact)); 1202 iSwapIndex = (PetscInt)(hashfact * nloc) % nloc; 1203 if (!bIndexSet[iSwapIndex] && iSwapIndex != Ii) { 1204 PetscInt iTemp = permute[iSwapIndex]; 1205 1206 permute[iSwapIndex] = permute[Ii]; 1207 permute[Ii] = iTemp; 1208 iTemp = degree[iSwapIndex]; 1209 degree[iSwapIndex] = degree[Ii]; 1210 degree[Ii] = iTemp; 1211 bIndexSet[iSwapIndex] = PETSC_TRUE; 1212 } 1213 } 1214 // apply minimum degree ordering -- NEW 1215 if (pc_gamg_agg->use_minimum_degree_ordering) { PetscCall(PetscSortIntWithArray(nloc, degree, permute)); } 1216 PetscCall(PetscFree(bIndexSet)); 1217 PetscCall(PetscRandomDestroy(&random)); 1218 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nloc, permute, PETSC_USE_POINTER, &perm)); 1219 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0)); 1220 // square graph 1221 if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels && pc_gamg_agg->use_aggressive_square_graph) { 1222 PetscCall(PCGAMGSquareGraph_GAMG(a_pc, Gmat1, &Gmat2)); 1223 } else Gmat2 = Gmat1; 1224 // switch to old MIS-1 for square graph 1225 if (pc_gamg->current_level < pc_gamg_agg->aggressive_coarsening_levels) { 1226 if (!pc_gamg_agg->use_aggressive_square_graph) PetscCall(MatCoarsenMISKSetDistance(pc_gamg_agg->crs, pc_gamg_agg->aggressive_mis_k)); // hardwire to MIS-2 1227 else PetscCall(MatCoarsenSetType(pc_gamg_agg->crs, MATCOARSENMIS)); // old MIS -- side effect 1228 } else if (pc_gamg_agg->use_aggressive_square_graph && pc_gamg_agg->aggressive_coarsening_levels > 0) { // we reset the MIS 1229 const char *prefix; 1230 1231 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)a_pc, &prefix)); 1232 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)pc_gamg_agg->crs, prefix)); 1233 PetscCall(MatCoarsenSetFromOptions(pc_gamg_agg->crs)); // get the default back on non-aggressive levels when square graph switched to old MIS 1234 } 1235 PetscCall(MatCoarsenSetAdjacency(pc_gamg_agg->crs, Gmat2)); 1236 PetscCall(MatCoarsenSetStrictAggs(pc_gamg_agg->crs, PETSC_TRUE)); 1237 PetscCall(MatCoarsenSetGreedyOrdering(pc_gamg_agg->crs, perm)); 1238 PetscCall(MatCoarsenApply(pc_gamg_agg->crs)); 1239 PetscCall(MatCoarsenGetData(pc_gamg_agg->crs, agg_lists)); /* output */ 1240 1241 PetscCall(ISDestroy(&perm)); 1242 PetscCall(PetscFree2(permute, degree)); 1243 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_MIS], 0, 0, 0, 0)); 1244 1245 if (Gmat2 != Gmat1) { // square graph, we need ghosts for selected 1246 PetscCoarsenData *llist = *agg_lists; 1247 1248 PetscCall(fixAggregatesWithSquare(a_pc, Gmat2, Gmat1, *agg_lists)); 1249 PetscCall(MatDestroy(&Gmat1)); 1250 *a_Gmat1 = Gmat2; /* output */ 1251 PetscCall(PetscCDSetMat(llist, *a_Gmat1)); /* Need a graph with ghosts here */ 1252 } 1253 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_COARSEN], 0, 0, 0, 0)); 1254 PetscFunctionReturn(PETSC_SUCCESS); 1255 } 1256 1257 /* 1258 PCGAMGConstructProlongator_AGG 1259 1260 Input Parameter: 1261 . pc - this 1262 . Amat - matrix on this fine level 1263 . Graph - used to get ghost data for nodes in 1264 . agg_lists - list of aggregates 1265 Output Parameter: 1266 . a_P_out - prolongation operator to the next level 1267 */ 1268 static PetscErrorCode PCGAMGConstructProlongator_AGG(PC pc, Mat Amat, PetscCoarsenData *agg_lists, Mat *a_P_out) 1269 { 1270 PC_MG *mg = (PC_MG *)pc->data; 1271 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1272 const PetscInt col_bs = pc_gamg->data_cell_cols; 1273 PetscInt Istart, Iend, nloc, ii, jj, kk, my0, nLocalSelected, bs; 1274 Mat Gmat, Prol; 1275 PetscMPIInt size; 1276 MPI_Comm comm; 1277 PetscReal *data_w_ghost; 1278 PetscInt myCrs0, nbnodes = 0, *flid_fgid; 1279 MatType mtype; 1280 1281 PetscFunctionBegin; 1282 PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm)); 1283 PetscCheck(col_bs >= 1, comm, PETSC_ERR_PLIB, "Column bs cannot be less than 1"); 1284 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0)); 1285 PetscCallMPI(MPI_Comm_size(comm, &size)); 1286 PetscCall(MatGetOwnershipRange(Amat, &Istart, &Iend)); 1287 PetscCall(MatGetBlockSize(Amat, &bs)); 1288 nloc = (Iend - Istart) / bs; 1289 my0 = Istart / bs; 1290 PetscCheck((Iend - Istart) % bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(Iend %" PetscInt_FMT " - Istart %" PetscInt_FMT ") not divisible by bs %" PetscInt_FMT, Iend, Istart, bs); 1291 PetscCall(PetscCDGetMat(agg_lists, &Gmat)); // get auxiliary matrix for ghost edges for size > 1 1292 1293 /* get 'nLocalSelected' */ 1294 for (ii = 0, nLocalSelected = 0; ii < nloc; ii++) { 1295 PetscBool ise; 1296 1297 /* filter out singletons 0 or 1? */ 1298 PetscCall(PetscCDIsEmptyAt(agg_lists, ii, &ise)); 1299 if (!ise) nLocalSelected++; 1300 } 1301 1302 /* create prolongator, create P matrix */ 1303 PetscCall(MatGetType(Amat, &mtype)); 1304 PetscCall(MatCreate(comm, &Prol)); 1305 PetscCall(MatSetSizes(Prol, nloc * bs, nLocalSelected * col_bs, PETSC_DETERMINE, PETSC_DETERMINE)); 1306 PetscCall(MatSetBlockSizes(Prol, bs, col_bs)); // should this be before MatSetSizes? 1307 PetscCall(MatSetType(Prol, mtype)); 1308 #if PetscDefined(HAVE_DEVICE) 1309 PetscBool flg; 1310 PetscCall(MatBoundToCPU(Amat, &flg)); 1311 PetscCall(MatBindToCPU(Prol, flg)); 1312 if (flg) PetscCall(MatSetBindingPropagates(Prol, PETSC_TRUE)); 1313 #endif 1314 PetscCall(MatSeqAIJSetPreallocation(Prol, col_bs, NULL)); 1315 PetscCall(MatMPIAIJSetPreallocation(Prol, col_bs, NULL, col_bs, NULL)); 1316 1317 /* can get all points "removed" */ 1318 PetscCall(MatGetSize(Prol, &kk, &ii)); 1319 if (!ii) { 1320 PetscCall(PetscInfo(pc, "%s: No selected points on coarse grid\n", ((PetscObject)pc)->prefix)); 1321 PetscCall(MatDestroy(&Prol)); 1322 *a_P_out = NULL; /* out */ 1323 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0)); 1324 PetscFunctionReturn(PETSC_SUCCESS); 1325 } 1326 PetscCall(PetscInfo(pc, "%s: New grid %" PetscInt_FMT " nodes\n", ((PetscObject)pc)->prefix, ii / col_bs)); 1327 PetscCall(MatGetOwnershipRangeColumn(Prol, &myCrs0, &kk)); 1328 1329 PetscCheck((kk - myCrs0) % col_bs == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(kk %" PetscInt_FMT " -myCrs0 %" PetscInt_FMT ") not divisible by col_bs %" PetscInt_FMT, kk, myCrs0, col_bs); 1330 myCrs0 = myCrs0 / col_bs; 1331 PetscCheck((kk / col_bs - myCrs0) == nLocalSelected, PETSC_COMM_SELF, PETSC_ERR_PLIB, "(kk %" PetscInt_FMT "/col_bs %" PetscInt_FMT " - myCrs0 %" PetscInt_FMT ") != nLocalSelected %" PetscInt_FMT ")", kk, col_bs, myCrs0, nLocalSelected); 1332 1333 /* create global vector of data in 'data_w_ghost' */ 1334 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0)); 1335 if (size > 1) { /* get ghost null space data */ 1336 PetscReal *tmp_gdata, *tmp_ldata, *tp2; 1337 1338 PetscCall(PetscMalloc1(nloc, &tmp_ldata)); 1339 for (jj = 0; jj < col_bs; jj++) { 1340 for (kk = 0; kk < bs; kk++) { 1341 PetscInt ii, stride; 1342 const PetscReal *tp = PetscSafePointerPlusOffset(pc_gamg->data, jj * bs * nloc + kk); 1343 1344 for (ii = 0; ii < nloc; ii++, tp += bs) tmp_ldata[ii] = *tp; 1345 1346 PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, tmp_ldata, &stride, &tmp_gdata)); 1347 1348 if (!jj && !kk) { /* now I know how many total nodes - allocate TODO: move below and do in one 'col_bs' call */ 1349 PetscCall(PetscMalloc1(stride * bs * col_bs, &data_w_ghost)); 1350 nbnodes = bs * stride; 1351 } 1352 tp2 = PetscSafePointerPlusOffset(data_w_ghost, jj * bs * stride + kk); 1353 for (ii = 0; ii < stride; ii++, tp2 += bs) *tp2 = tmp_gdata[ii]; 1354 PetscCall(PetscFree(tmp_gdata)); 1355 } 1356 } 1357 PetscCall(PetscFree(tmp_ldata)); 1358 } else { 1359 nbnodes = bs * nloc; 1360 data_w_ghost = pc_gamg->data; 1361 } 1362 1363 /* get 'flid_fgid' TODO - move up to get 'stride' and do get null space data above in one step (jj loop) */ 1364 if (size > 1) { 1365 PetscReal *fid_glid_loc, *fiddata; 1366 PetscInt stride; 1367 1368 PetscCall(PetscMalloc1(nloc, &fid_glid_loc)); 1369 for (kk = 0; kk < nloc; kk++) fid_glid_loc[kk] = (PetscReal)(my0 + kk); 1370 PetscCall(PCGAMGGetDataWithGhosts(Gmat, 1, fid_glid_loc, &stride, &fiddata)); 1371 PetscCall(PetscMalloc1(stride, &flid_fgid)); /* copy real data to in */ 1372 for (kk = 0; kk < stride; kk++) flid_fgid[kk] = (PetscInt)fiddata[kk]; 1373 PetscCall(PetscFree(fiddata)); 1374 1375 PetscCheck(stride == nbnodes / bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "stride %" PetscInt_FMT " != nbnodes %" PetscInt_FMT "/bs %" PetscInt_FMT, stride, nbnodes, bs); 1376 PetscCall(PetscFree(fid_glid_loc)); 1377 } else { 1378 PetscCall(PetscMalloc1(nloc, &flid_fgid)); 1379 for (kk = 0; kk < nloc; kk++) flid_fgid[kk] = my0 + kk; 1380 } 1381 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLA], 0, 0, 0, 0)); 1382 /* get P0 */ 1383 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0)); 1384 { 1385 PetscReal *data_out = NULL; 1386 1387 PetscCall(formProl0(agg_lists, bs, col_bs, myCrs0, nbnodes, data_w_ghost, flid_fgid, &data_out, Prol)); 1388 PetscCall(PetscFree(pc_gamg->data)); 1389 1390 pc_gamg->data = data_out; 1391 pc_gamg->data_cell_rows = col_bs; 1392 pc_gamg->data_sz = col_bs * col_bs * nLocalSelected; 1393 } 1394 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROLB], 0, 0, 0, 0)); 1395 if (size > 1) PetscCall(PetscFree(data_w_ghost)); 1396 PetscCall(PetscFree(flid_fgid)); 1397 1398 *a_P_out = Prol; /* out */ 1399 PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_initial_prolongation")); 1400 1401 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_PROL], 0, 0, 0, 0)); 1402 PetscFunctionReturn(PETSC_SUCCESS); 1403 } 1404 1405 /* 1406 PCGAMGOptimizeProlongator_AGG - given the initial prolongator optimizes it by smoothed aggregation pc_gamg_agg->nsmooths times 1407 1408 Input Parameter: 1409 . pc - this 1410 . Amat - matrix on this fine level 1411 In/Output Parameter: 1412 . a_P - prolongation operator to the next level 1413 */ 1414 static PetscErrorCode PCGAMGOptimizeProlongator_AGG(PC pc, Mat Amat, Mat *a_P) 1415 { 1416 PC_MG *mg = (PC_MG *)pc->data; 1417 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1418 PC_GAMG_AGG *pc_gamg_agg = (PC_GAMG_AGG *)pc_gamg->subctx; 1419 PetscInt jj; 1420 Mat Prol = *a_P; 1421 MPI_Comm comm; 1422 KSP eksp; 1423 Vec bb, xx; 1424 PC epc; 1425 PetscReal alpha, emax, emin; 1426 1427 PetscFunctionBegin; 1428 PetscCall(PetscObjectGetComm((PetscObject)Amat, &comm)); 1429 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0)); 1430 1431 /* compute maximum singular value of operator to be used in smoother */ 1432 if (0 < pc_gamg_agg->nsmooths) { 1433 /* get eigen estimates */ 1434 if (pc_gamg->emax > 0) { 1435 emin = pc_gamg->emin; 1436 emax = pc_gamg->emax; 1437 } else { 1438 const char *prefix; 1439 1440 PetscCall(MatCreateVecs(Amat, &bb, NULL)); 1441 PetscCall(MatCreateVecs(Amat, &xx, NULL)); 1442 PetscCall(KSPSetNoisy_Private(Amat, bb)); 1443 1444 PetscCall(KSPCreate(comm, &eksp)); 1445 PetscCall(KSPSetNestLevel(eksp, pc->kspnestlevel)); 1446 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 1447 PetscCall(KSPSetOptionsPrefix(eksp, prefix)); 1448 PetscCall(KSPAppendOptionsPrefix(eksp, "pc_gamg_esteig_")); 1449 { 1450 PetscBool isset, sflg; 1451 1452 PetscCall(MatIsSPDKnown(Amat, &isset, &sflg)); 1453 if (isset && sflg) PetscCall(KSPSetType(eksp, KSPCG)); 1454 } 1455 PetscCall(KSPSetErrorIfNotConverged(eksp, pc->erroriffailure)); 1456 PetscCall(KSPSetNormType(eksp, KSP_NORM_NONE)); 1457 1458 PetscCall(KSPSetInitialGuessNonzero(eksp, PETSC_FALSE)); 1459 PetscCall(KSPSetOperators(eksp, Amat, Amat)); 1460 1461 PetscCall(KSPGetPC(eksp, &epc)); 1462 PetscCall(PCSetType(epc, PCJACOBI)); /* smoother in smoothed agg. */ 1463 1464 PetscCall(KSPSetTolerances(eksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, 10)); // 10 is safer, but 5 is often fine, can override with -pc_gamg_esteig_ksp_max_it -mg_levels_ksp_chebyshev_esteig 0,0.25,0,1.2 1465 1466 PetscCall(KSPSetFromOptions(eksp)); 1467 PetscCall(KSPSetComputeSingularValues(eksp, PETSC_TRUE)); 1468 PetscCall(KSPSolve(eksp, bb, xx)); 1469 PetscCall(KSPCheckSolve(eksp, pc, xx)); 1470 1471 PetscCall(KSPComputeExtremeSingularValues(eksp, &emax, &emin)); 1472 PetscCall(PetscInfo(pc, "%s: Smooth P0: max eigen=%e min=%e PC=%s\n", ((PetscObject)pc)->prefix, (double)emax, (double)emin, PCJACOBI)); 1473 PetscCall(VecDestroy(&xx)); 1474 PetscCall(VecDestroy(&bb)); 1475 PetscCall(KSPDestroy(&eksp)); 1476 } 1477 if (pc_gamg->use_sa_esteig) { 1478 mg->min_eigen_DinvA[pc_gamg->current_level] = emin; 1479 mg->max_eigen_DinvA[pc_gamg->current_level] = emax; 1480 PetscCall(PetscInfo(pc, "%s: Smooth P0: level %" PetscInt_FMT ", cache spectra %g %g\n", ((PetscObject)pc)->prefix, pc_gamg->current_level, (double)emin, (double)emax)); 1481 } else { 1482 mg->min_eigen_DinvA[pc_gamg->current_level] = 0; 1483 mg->max_eigen_DinvA[pc_gamg->current_level] = 0; 1484 } 1485 } else { 1486 mg->min_eigen_DinvA[pc_gamg->current_level] = 0; 1487 mg->max_eigen_DinvA[pc_gamg->current_level] = 0; 1488 } 1489 1490 /* smooth P0 */ 1491 if (pc_gamg_agg->nsmooths > 0) { 1492 Vec diag; 1493 1494 /* TODO: Set a PCFailedReason and exit the building of the AMG preconditioner */ 1495 PetscCheck(emax != 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Computed maximum singular value as zero"); 1496 1497 PetscCall(MatCreateVecs(Amat, &diag, NULL)); 1498 PetscCall(MatGetDiagonal(Amat, diag)); /* effectively PCJACOBI */ 1499 PetscCall(VecReciprocal(diag)); 1500 1501 for (jj = 0; jj < pc_gamg_agg->nsmooths; jj++) { 1502 Mat tMat; 1503 1504 PetscCall(PetscLogEventBegin(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0)); 1505 /* 1506 Smooth aggregation on the prolongator 1507 1508 P_{i} := (I - 1.4/emax D^{-1}A) P_i\{i-1} 1509 */ 1510 PetscCall(PetscLogEventBegin(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0)); 1511 PetscCall(MatMatMult(Amat, Prol, MAT_INITIAL_MATRIX, PETSC_CURRENT, &tMat)); 1512 PetscCall(PetscLogEventEnd(petsc_gamg_setup_matmat_events[pc_gamg->current_level][2], 0, 0, 0, 0)); 1513 PetscCall(MatProductClear(tMat)); 1514 PetscCall(MatDiagonalScale(tMat, diag, NULL)); 1515 1516 /* TODO: Document the 1.4 and don't hardwire it in this routine */ 1517 alpha = -1.4 / emax; 1518 PetscCall(MatAYPX(tMat, alpha, Prol, SUBSET_NONZERO_PATTERN)); 1519 PetscCall(MatDestroy(&Prol)); 1520 Prol = tMat; 1521 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPTSM], 0, 0, 0, 0)); 1522 } 1523 PetscCall(VecDestroy(&diag)); 1524 } 1525 PetscCall(PetscLogEventEnd(petsc_gamg_setup_events[GAMG_OPT], 0, 0, 0, 0)); 1526 PetscCall(MatViewFromOptions(Prol, NULL, "-pc_gamg_agg_view_prolongation")); 1527 *a_P = Prol; 1528 PetscFunctionReturn(PETSC_SUCCESS); 1529 } 1530 1531 /*MC 1532 PCGAMGAGG - Smooth aggregation, {cite}`vanek1996algebraic`, {cite}`vanek2001convergence`, variant of PETSc's algebraic multigrid (`PCGAMG`) preconditioner 1533 1534 Options Database Keys: 1535 + -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation to construct prolongation 1536 . -pc_gamg_aggressive_coarsening <n,default=1> - number of aggressive coarsening (MIS-2) levels from finest. 1537 . -pc_gamg_aggressive_square_graph <bool,default=false> - Use square graph (A'A) or MIS-k (k=2) for aggressive coarsening 1538 . -pc_gamg_mis_k_minimum_degree_ordering <bool,default=true> - Use minimum degree ordering in greedy MIS algorithm 1539 . -pc_gamg_pc_gamg_asm_hem_aggs <n,default=0> - Number of HEM aggregation steps for ASM smoother 1540 - -pc_gamg_aggressive_mis_k <n,default=2> - Number (k) distance in MIS coarsening (>2 is 'aggressive') 1541 1542 Level: intermediate 1543 1544 Notes: 1545 To obtain good performance for `PCGAMG` for vector valued problems you must 1546 call `MatSetBlockSize()` to indicate the number of degrees of freedom per grid point. 1547 Call `MatSetNearNullSpace()` (or `PCSetCoordinates()` if solving the equations of elasticity) to indicate the near null space of the operator 1548 1549 The many options for `PCMG` and `PCGAMG` such as controlling the smoothers on each level etc. also work for `PCGAMGAGG` 1550 1551 .seealso: `PCGAMG`, [the Users Manual section on PCGAMG](sec_amg), [the Users Manual section on PCMG](sec_mg), [](ch_ksp), `PCCreate()`, `PCSetType()`, 1552 `MatSetBlockSize()`, `PCMGType`, `PCSetCoordinates()`, `MatSetNearNullSpace()`, `PCGAMGSetType()`, 1553 `PCGAMGAGG`, `PCGAMGGEO`, `PCGAMGCLASSICAL`, `PCGAMGSetProcEqLim()`, `PCGAMGSetCoarseEqLim()`, `PCGAMGSetRepartition()`, `PCGAMGRegister()`, 1554 `PCGAMGSetReuseInterpolation()`, `PCGAMGASMSetUseAggs()`, `PCGAMGSetParallelCoarseGridSolve()`, `PCGAMGSetNlevels()`, `PCGAMGSetThreshold()`, 1555 `PCGAMGGetType()`, `PCGAMGSetUseSAEstEig()` 1556 M*/ 1557 PetscErrorCode PCCreateGAMG_AGG(PC pc) 1558 { 1559 PC_MG *mg = (PC_MG *)pc->data; 1560 PC_GAMG *pc_gamg = (PC_GAMG *)mg->innerctx; 1561 PC_GAMG_AGG *pc_gamg_agg; 1562 1563 PetscFunctionBegin; 1564 /* create sub context for SA */ 1565 PetscCall(PetscNew(&pc_gamg_agg)); 1566 pc_gamg->subctx = pc_gamg_agg; 1567 1568 pc_gamg->ops->setfromoptions = PCSetFromOptions_GAMG_AGG; 1569 pc_gamg->ops->destroy = PCDestroy_GAMG_AGG; 1570 /* reset does not do anything; setup not virtual */ 1571 1572 /* set internal function pointers */ 1573 pc_gamg->ops->creategraph = PCGAMGCreateGraph_AGG; 1574 pc_gamg->ops->coarsen = PCGAMGCoarsen_AGG; 1575 pc_gamg->ops->prolongator = PCGAMGConstructProlongator_AGG; 1576 pc_gamg->ops->optprolongator = PCGAMGOptimizeProlongator_AGG; 1577 pc_gamg->ops->createdefaultdata = PCSetData_AGG; 1578 pc_gamg->ops->view = PCView_GAMG_AGG; 1579 1580 pc_gamg_agg->nsmooths = 1; 1581 pc_gamg_agg->aggressive_coarsening_levels = 1; 1582 pc_gamg_agg->use_aggressive_square_graph = PETSC_TRUE; 1583 pc_gamg_agg->use_minimum_degree_ordering = PETSC_FALSE; 1584 pc_gamg_agg->use_low_mem_filter = PETSC_FALSE; 1585 pc_gamg_agg->aggressive_mis_k = 2; 1586 pc_gamg_agg->graph_symmetrize = PETSC_TRUE; 1587 1588 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetNSmooths_C", PCGAMGSetNSmooths_AGG)); 1589 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveLevels_C", PCGAMGSetAggressiveLevels_AGG)); 1590 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetAggressiveSquareGraph_C", PCGAMGSetAggressiveSquareGraph_AGG)); 1591 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetMinDegreeOrdering_C", PCGAMGMISkSetMinDegreeOrdering_AGG)); 1592 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetLowMemoryFilter_C", PCGAMGSetLowMemoryFilter_AGG)); 1593 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGMISkSetAggressive_C", PCGAMGMISkSetAggressive_AGG)); 1594 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGAMGSetGraphSymmetrize_C", PCGAMGSetGraphSymmetrize_AGG)); 1595 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_AGG)); 1596 PetscFunctionReturn(PETSC_SUCCESS); 1597 } 1598