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