1 /* 2 Creates a matrix class for using the Neumann-Neumann type preconditioners. 3 This stores the matrices in globally unassembled form. Each processor 4 assembles only its local Neumann problem and the parallel matrix vector 5 product is handled "implicitly". 6 7 Currently this allows for only one subdomain per processor. 8 */ 9 10 #include <petsc/private/matisimpl.h> /*I "petscmat.h" I*/ 11 #include <petsc/private/sfimpl.h> 12 #include <petsc/private/vecimpl.h> 13 #include <petsc/private/hashseti.h> 14 15 #define MATIS_MAX_ENTRIES_INSERTION 2048 16 static PetscErrorCode MatSetValuesLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode); 17 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat, PetscInt, const PetscInt *, PetscInt, const PetscInt *, const PetscScalar *, InsertMode); 18 static PetscErrorCode MatISSetUpScatters_Private(Mat); 19 20 static PetscErrorCode MatISContainerDestroyPtAP_Private(void *ptr) { 21 MatISPtAP ptap = (MatISPtAP)ptr; 22 23 PetscFunctionBegin; 24 PetscCall(MatDestroySubMatrices(ptap->ris1 ? 2 : 1, &ptap->lP)); 25 PetscCall(ISDestroy(&ptap->cis0)); 26 PetscCall(ISDestroy(&ptap->cis1)); 27 PetscCall(ISDestroy(&ptap->ris0)); 28 PetscCall(ISDestroy(&ptap->ris1)); 29 PetscCall(PetscFree(ptap)); 30 PetscFunctionReturn(0); 31 } 32 33 static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C) { 34 MatISPtAP ptap; 35 Mat_IS *matis = (Mat_IS *)A->data; 36 Mat lA, lC; 37 MatReuse reuse; 38 IS ris[2], cis[2]; 39 PetscContainer c; 40 PetscInt n; 41 42 PetscFunctionBegin; 43 PetscCall(PetscObjectQuery((PetscObject)C, "_MatIS_PtAP", (PetscObject *)&c)); 44 PetscCheck(c, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing PtAP information"); 45 PetscCall(PetscContainerGetPointer(c, (void **)&ptap)); 46 ris[0] = ptap->ris0; 47 ris[1] = ptap->ris1; 48 cis[0] = ptap->cis0; 49 cis[1] = ptap->cis1; 50 n = ptap->ris1 ? 2 : 1; 51 reuse = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX; 52 PetscCall(MatCreateSubMatrices(P, n, ris, cis, reuse, &ptap->lP)); 53 54 PetscCall(MatISGetLocalMat(A, &lA)); 55 PetscCall(MatISGetLocalMat(C, &lC)); 56 if (ptap->ris1) { /* unsymmetric A mapping */ 57 Mat lPt; 58 59 PetscCall(MatTranspose(ptap->lP[1], MAT_INITIAL_MATRIX, &lPt)); 60 PetscCall(MatMatMatMult(lPt, lA, ptap->lP[0], reuse, ptap->fill, &lC)); 61 if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)(A), "_MatIS_PtAP_l2l", (PetscObject)lPt)); 62 PetscCall(MatDestroy(&lPt)); 63 } else { 64 PetscCall(MatPtAP(lA, ptap->lP[0], reuse, ptap->fill, &lC)); 65 if (matis->storel2l) PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP_l2l", (PetscObject)ptap->lP[0])); 66 } 67 if (reuse == MAT_INITIAL_MATRIX) { 68 PetscCall(MatISSetLocalMat(C, lC)); 69 PetscCall(MatDestroy(&lC)); 70 } 71 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 72 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 73 PetscFunctionReturn(0); 74 } 75 76 static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT, IS *cis) { 77 Mat Po, Pd; 78 IS zd, zo; 79 const PetscInt *garray; 80 PetscInt *aux, i, bs; 81 PetscInt dc, stc, oc, ctd, cto; 82 PetscBool ismpiaij, ismpibaij, isseqaij, isseqbaij; 83 MPI_Comm comm; 84 85 PetscFunctionBegin; 86 PetscValidHeaderSpecific(PT, MAT_CLASSID, 1); 87 PetscValidPointer(cis, 2); 88 PetscCall(PetscObjectGetComm((PetscObject)PT, &comm)); 89 bs = 1; 90 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIAIJ, &ismpiaij)); 91 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATMPIBAIJ, &ismpibaij)); 92 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT, MATSEQAIJ, &isseqaij)); 93 PetscCall(PetscObjectTypeCompare((PetscObject)PT, MATSEQBAIJ, &isseqbaij)); 94 if (isseqaij || isseqbaij) { 95 Pd = PT; 96 Po = NULL; 97 garray = NULL; 98 } else if (ismpiaij) { 99 PetscCall(MatMPIAIJGetSeqAIJ(PT, &Pd, &Po, &garray)); 100 } else if (ismpibaij) { 101 PetscCall(MatMPIBAIJGetSeqBAIJ(PT, &Pd, &Po, &garray)); 102 PetscCall(MatGetBlockSize(PT, &bs)); 103 } else SETERRQ(comm, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)(PT))->type_name); 104 105 /* identify any null columns in Pd or Po */ 106 /* We use a tolerance comparison since it may happen that, with geometric multigrid, 107 some of the columns are not really zero, but very close to */ 108 zo = zd = NULL; 109 if (Po) PetscCall(MatFindNonzeroRowsOrCols_Basic(Po, PETSC_TRUE, PETSC_SMALL, &zo)); 110 PetscCall(MatFindNonzeroRowsOrCols_Basic(Pd, PETSC_TRUE, PETSC_SMALL, &zd)); 111 112 PetscCall(MatGetLocalSize(PT, NULL, &dc)); 113 PetscCall(MatGetOwnershipRangeColumn(PT, &stc, NULL)); 114 if (Po) PetscCall(MatGetLocalSize(Po, NULL, &oc)); 115 else oc = 0; 116 PetscCall(PetscMalloc1((dc + oc) / bs, &aux)); 117 if (zd) { 118 const PetscInt *idxs; 119 PetscInt nz; 120 121 /* this will throw an error if bs is not valid */ 122 PetscCall(ISSetBlockSize(zd, bs)); 123 PetscCall(ISGetLocalSize(zd, &nz)); 124 PetscCall(ISGetIndices(zd, &idxs)); 125 ctd = nz / bs; 126 for (i = 0; i < ctd; i++) aux[i] = (idxs[bs * i] + stc) / bs; 127 PetscCall(ISRestoreIndices(zd, &idxs)); 128 } else { 129 ctd = dc / bs; 130 for (i = 0; i < ctd; i++) aux[i] = i + stc / bs; 131 } 132 if (zo) { 133 const PetscInt *idxs; 134 PetscInt nz; 135 136 /* this will throw an error if bs is not valid */ 137 PetscCall(ISSetBlockSize(zo, bs)); 138 PetscCall(ISGetLocalSize(zo, &nz)); 139 PetscCall(ISGetIndices(zo, &idxs)); 140 cto = nz / bs; 141 for (i = 0; i < cto; i++) aux[i + ctd] = garray[idxs[bs * i] / bs]; 142 PetscCall(ISRestoreIndices(zo, &idxs)); 143 } else { 144 cto = oc / bs; 145 for (i = 0; i < cto; i++) aux[i + ctd] = garray[i]; 146 } 147 PetscCall(ISCreateBlock(comm, bs, ctd + cto, aux, PETSC_OWN_POINTER, cis)); 148 PetscCall(ISDestroy(&zd)); 149 PetscCall(ISDestroy(&zo)); 150 PetscFunctionReturn(0); 151 } 152 153 static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A, Mat P, PetscReal fill, Mat C) { 154 Mat PT, lA; 155 MatISPtAP ptap; 156 ISLocalToGlobalMapping Crl2g, Ccl2g, rl2g, cl2g; 157 PetscContainer c; 158 MatType lmtype; 159 const PetscInt *garray; 160 PetscInt ibs, N, dc; 161 MPI_Comm comm; 162 163 PetscFunctionBegin; 164 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 165 PetscCall(MatSetType(C, MATIS)); 166 PetscCall(MatISGetLocalMat(A, &lA)); 167 PetscCall(MatGetType(lA, &lmtype)); 168 PetscCall(MatISSetLocalMatType(C, lmtype)); 169 PetscCall(MatGetSize(P, NULL, &N)); 170 PetscCall(MatGetLocalSize(P, NULL, &dc)); 171 PetscCall(MatSetSizes(C, dc, dc, N, N)); 172 /* Not sure about this 173 PetscCall(MatGetBlockSizes(P,NULL,&ibs)); 174 PetscCall(MatSetBlockSize(*C,ibs)); 175 */ 176 177 PetscCall(PetscNew(&ptap)); 178 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c)); 179 PetscCall(PetscContainerSetPointer(c, ptap)); 180 PetscCall(PetscContainerSetUserDestroy(c, MatISContainerDestroyPtAP_Private)); 181 PetscCall(PetscObjectCompose((PetscObject)C, "_MatIS_PtAP", (PetscObject)c)); 182 PetscCall(PetscContainerDestroy(&c)); 183 ptap->fill = fill; 184 185 PetscCall(MatISGetLocalToGlobalMapping(A, &rl2g, &cl2g)); 186 187 PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &ibs)); 188 PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &N)); 189 PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &garray)); 190 PetscCall(ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris0)); 191 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &garray)); 192 193 PetscCall(MatCreateSubMatrix(P, ptap->ris0, NULL, MAT_INITIAL_MATRIX, &PT)); 194 PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis0)); 195 PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis0, &Ccl2g)); 196 PetscCall(MatDestroy(&PT)); 197 198 Crl2g = NULL; 199 if (rl2g != cl2g) { /* unsymmetric A mapping */ 200 PetscBool same, lsame = PETSC_FALSE; 201 PetscInt N1, ibs1; 202 203 PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &N1)); 204 PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &ibs1)); 205 PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &garray)); 206 PetscCall(ISCreateBlock(comm, ibs, N / ibs, garray, PETSC_COPY_VALUES, &ptap->ris1)); 207 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &garray)); 208 if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */ 209 const PetscInt *i1, *i2; 210 211 PetscCall(ISBlockGetIndices(ptap->ris0, &i1)); 212 PetscCall(ISBlockGetIndices(ptap->ris1, &i2)); 213 PetscCall(PetscArraycmp(i1, i2, N, &lsame)); 214 } 215 PetscCall(MPIU_Allreduce(&lsame, &same, 1, MPIU_BOOL, MPI_LAND, comm)); 216 if (same) { 217 PetscCall(ISDestroy(&ptap->ris1)); 218 } else { 219 PetscCall(MatCreateSubMatrix(P, ptap->ris1, NULL, MAT_INITIAL_MATRIX, &PT)); 220 PetscCall(MatGetNonzeroColumnsLocal_Private(PT, &ptap->cis1)); 221 PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis1, &Crl2g)); 222 PetscCall(MatDestroy(&PT)); 223 } 224 } 225 /* Not sure about this 226 if (!Crl2g) { 227 PetscCall(MatGetBlockSize(C,&ibs)); 228 PetscCall(ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs)); 229 } 230 */ 231 PetscCall(MatSetLocalToGlobalMapping(C, Crl2g ? Crl2g : Ccl2g, Ccl2g)); 232 PetscCall(ISLocalToGlobalMappingDestroy(&Crl2g)); 233 PetscCall(ISLocalToGlobalMappingDestroy(&Ccl2g)); 234 235 C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ; 236 PetscFunctionReturn(0); 237 } 238 239 /* ----------------------------------------- */ 240 static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C) { 241 Mat_Product *product = C->product; 242 Mat A = product->A, P = product->B; 243 PetscReal fill = product->fill; 244 245 PetscFunctionBegin; 246 PetscCall(MatPtAPSymbolic_IS_XAIJ(A, P, fill, C)); 247 C->ops->productnumeric = MatProductNumeric_PtAP; 248 PetscFunctionReturn(0); 249 } 250 251 static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C) { 252 PetscFunctionBegin; 253 C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ; 254 PetscFunctionReturn(0); 255 } 256 257 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C) { 258 Mat_Product *product = C->product; 259 260 PetscFunctionBegin; 261 if (product->type == MATPRODUCT_PtAP) PetscCall(MatProductSetFromOptions_IS_XAIJ_PtAP(C)); 262 PetscFunctionReturn(0); 263 } 264 265 /* ----------------------------------------- */ 266 static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr) { 267 MatISLocalFields lf = (MatISLocalFields)ptr; 268 PetscInt i; 269 270 PetscFunctionBegin; 271 for (i = 0; i < lf->nr; i++) PetscCall(ISDestroy(&lf->rf[i])); 272 for (i = 0; i < lf->nc; i++) PetscCall(ISDestroy(&lf->cf[i])); 273 PetscCall(PetscFree2(lf->rf, lf->cf)); 274 PetscCall(PetscFree(lf)); 275 PetscFunctionReturn(0); 276 } 277 278 static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat) { 279 Mat B, lB; 280 281 PetscFunctionBegin; 282 if (reuse != MAT_REUSE_MATRIX) { 283 ISLocalToGlobalMapping rl2g, cl2g; 284 PetscInt bs; 285 IS is; 286 287 PetscCall(MatGetBlockSize(A, &bs)); 288 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->n / bs, 0, 1, &is)); 289 if (bs > 1) { 290 IS is2; 291 PetscInt i, *aux; 292 293 PetscCall(ISGetLocalSize(is, &i)); 294 PetscCall(ISGetIndices(is, (const PetscInt **)&aux)); 295 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2)); 296 PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux)); 297 PetscCall(ISDestroy(&is)); 298 is = is2; 299 } 300 PetscCall(ISSetIdentity(is)); 301 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g)); 302 PetscCall(ISDestroy(&is)); 303 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->n / bs, 0, 1, &is)); 304 if (bs > 1) { 305 IS is2; 306 PetscInt i, *aux; 307 308 PetscCall(ISGetLocalSize(is, &i)); 309 PetscCall(ISGetIndices(is, (const PetscInt **)&aux)); 310 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, i, aux, PETSC_COPY_VALUES, &is2)); 311 PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux)); 312 PetscCall(ISDestroy(&is)); 313 is = is2; 314 } 315 PetscCall(ISSetIdentity(is)); 316 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g)); 317 PetscCall(ISDestroy(&is)); 318 PetscCall(MatCreateIS(PetscObjectComm((PetscObject)A), bs, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N, rl2g, cl2g, &B)); 319 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 320 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 321 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &lB)); 322 if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 323 } else { 324 B = *newmat; 325 PetscCall(PetscObjectReference((PetscObject)A)); 326 lB = A; 327 } 328 PetscCall(MatISSetLocalMat(B, lB)); 329 PetscCall(MatDestroy(&lB)); 330 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 331 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 332 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B)); 333 PetscFunctionReturn(0); 334 } 335 336 static PetscErrorCode MatISScaleDisassembling_Private(Mat A) { 337 Mat_IS *matis = (Mat_IS *)(A->data); 338 PetscScalar *aa; 339 const PetscInt *ii, *jj; 340 PetscInt i, n, m; 341 PetscInt *ecount, **eneighs; 342 PetscBool flg; 343 344 PetscFunctionBegin; 345 PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg)); 346 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure"); 347 PetscCall(ISLocalToGlobalMappingGetNodeInfo(matis->rmapping, &n, &ecount, &eneighs)); 348 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected %" PetscInt_FMT " != %" PetscInt_FMT, m, n); 349 PetscCall(MatSeqAIJGetArray(matis->A, &aa)); 350 for (i = 0; i < n; i++) { 351 if (ecount[i] > 1) { 352 PetscInt j; 353 354 for (j = ii[i]; j < ii[i + 1]; j++) { 355 PetscInt i2 = jj[j], p, p2; 356 PetscReal scal = 0.0; 357 358 for (p = 0; p < ecount[i]; p++) { 359 for (p2 = 0; p2 < ecount[i2]; p2++) { 360 if (eneighs[i][p] == eneighs[i2][p2]) { 361 scal += 1.0; 362 break; 363 } 364 } 365 } 366 if (scal) aa[j] /= scal; 367 } 368 } 369 } 370 PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping, &n, &ecount, &eneighs)); 371 PetscCall(MatSeqAIJRestoreArray(matis->A, &aa)); 372 PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &m, &ii, &jj, &flg)); 373 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure"); 374 PetscFunctionReturn(0); 375 } 376 377 typedef enum { 378 MAT_IS_DISASSEMBLE_L2G_NATURAL, 379 MAT_IS_DISASSEMBLE_L2G_MAT, 380 MAT_IS_DISASSEMBLE_L2G_ND 381 } MatISDisassemblel2gType; 382 383 static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g) { 384 Mat Ad, Ao; 385 IS is, ndmap, ndsub; 386 MPI_Comm comm; 387 const PetscInt *garray, *ndmapi; 388 PetscInt bs, i, cnt, nl, *ncount, *ndmapc; 389 PetscBool ismpiaij, ismpibaij; 390 const char *const MatISDisassemblel2gTypes[] = {"NATURAL", "MAT", "ND", "MatISDisassemblel2gType", "MAT_IS_DISASSEMBLE_L2G_", NULL}; 391 MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL; 392 MatPartitioning part; 393 PetscSF sf; 394 PetscObject dm; 395 396 PetscFunctionBegin; 397 PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "MatIS l2g disassembling options", "Mat"); 398 PetscCall(PetscOptionsEnum("-mat_is_disassemble_l2g_type", "Type of local-to-global mapping to be used for disassembling", "MatISDisassemblel2gType", MatISDisassemblel2gTypes, (PetscEnum)mode, (PetscEnum *)&mode, NULL)); 399 PetscOptionsEnd(); 400 if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) { 401 PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL)); 402 PetscFunctionReturn(0); 403 } 404 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 405 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij)); 406 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij)); 407 PetscCall(MatGetBlockSize(A, &bs)); 408 switch (mode) { 409 case MAT_IS_DISASSEMBLE_L2G_ND: 410 PetscCall(MatPartitioningCreate(comm, &part)); 411 PetscCall(MatPartitioningSetAdjacency(part, A)); 412 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part, ((PetscObject)A)->prefix)); 413 PetscCall(MatPartitioningSetFromOptions(part)); 414 PetscCall(MatPartitioningApplyND(part, &ndmap)); 415 PetscCall(MatPartitioningDestroy(&part)); 416 PetscCall(ISBuildTwoSided(ndmap, NULL, &ndsub)); 417 PetscCall(MatMPIAIJSetUseScalableIncreaseOverlap(A, PETSC_TRUE)); 418 PetscCall(MatIncreaseOverlap(A, 1, &ndsub, 1)); 419 PetscCall(ISLocalToGlobalMappingCreateIS(ndsub, l2g)); 420 421 /* it may happen that a separator node is not properly shared */ 422 PetscCall(ISLocalToGlobalMappingGetNodeInfo(*l2g, &nl, &ncount, NULL)); 423 PetscCall(PetscSFCreate(comm, &sf)); 424 PetscCall(ISLocalToGlobalMappingGetIndices(*l2g, &garray)); 425 PetscCall(PetscSFSetGraphLayout(sf, A->rmap, nl, NULL, PETSC_OWN_POINTER, garray)); 426 PetscCall(ISLocalToGlobalMappingRestoreIndices(*l2g, &garray)); 427 PetscCall(PetscCalloc1(A->rmap->n, &ndmapc)); 428 PetscCall(PetscSFReduceBegin(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE)); 429 PetscCall(PetscSFReduceEnd(sf, MPIU_INT, ncount, ndmapc, MPI_REPLACE)); 430 PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(*l2g, NULL, &ncount, NULL)); 431 PetscCall(ISGetIndices(ndmap, &ndmapi)); 432 for (i = 0, cnt = 0; i < A->rmap->n; i++) 433 if (ndmapi[i] < 0 && ndmapc[i] < 2) cnt++; 434 435 PetscCall(MPIU_Allreduce(&cnt, &i, 1, MPIU_INT, MPI_MAX, comm)); 436 if (i) { /* we detected isolated separator nodes */ 437 Mat A2, A3; 438 IS *workis, is2; 439 PetscScalar *vals; 440 PetscInt gcnt = i, *dnz, *onz, j, *lndmapi; 441 ISLocalToGlobalMapping ll2g; 442 PetscBool flg; 443 const PetscInt *ii, *jj; 444 445 /* communicate global id of separators */ 446 MatPreallocateBegin(comm, A->rmap->n, A->cmap->n, dnz, onz); 447 for (i = 0, cnt = 0; i < A->rmap->n; i++) dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1; 448 449 PetscCall(PetscMalloc1(nl, &lndmapi)); 450 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE)); 451 452 /* compute adjacency of isolated separators node */ 453 PetscCall(PetscMalloc1(gcnt, &workis)); 454 for (i = 0, cnt = 0; i < A->rmap->n; i++) { 455 if (ndmapi[i] < 0 && ndmapc[i] < 2) PetscCall(ISCreateStride(comm, 1, i + A->rmap->rstart, 1, &workis[cnt++])); 456 } 457 for (i = cnt; i < gcnt; i++) PetscCall(ISCreateStride(comm, 0, 0, 1, &workis[i])); 458 for (i = 0; i < gcnt; i++) { 459 PetscCall(PetscObjectSetName((PetscObject)workis[i], "ISOLATED")); 460 PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators")); 461 } 462 463 /* no communications since all the ISes correspond to locally owned rows */ 464 PetscCall(MatIncreaseOverlap(A, gcnt, workis, 1)); 465 466 /* end communicate global id of separators */ 467 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dnz, lndmapi, MPI_REPLACE)); 468 469 /* communicate new layers : create a matrix and transpose it */ 470 PetscCall(PetscArrayzero(dnz, A->rmap->n)); 471 PetscCall(PetscArrayzero(onz, A->rmap->n)); 472 for (i = 0, j = 0; i < A->rmap->n; i++) { 473 if (ndmapi[i] < 0 && ndmapc[i] < 2) { 474 const PetscInt *idxs; 475 PetscInt s; 476 477 PetscCall(ISGetLocalSize(workis[j], &s)); 478 PetscCall(ISGetIndices(workis[j], &idxs)); 479 PetscCall(MatPreallocateSet(i + A->rmap->rstart, s, idxs, dnz, onz)); 480 j++; 481 } 482 } 483 PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt); 484 485 for (i = 0; i < gcnt; i++) { 486 PetscCall(PetscObjectSetName((PetscObject)workis[i], "EXTENDED")); 487 PetscCall(ISViewFromOptions(workis[i], NULL, "-view_isolated_separators")); 488 } 489 490 for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j, dnz[i] + onz[i]); 491 PetscCall(PetscMalloc1(j, &vals)); 492 for (i = 0; i < j; i++) vals[i] = 1.0; 493 494 PetscCall(MatCreate(comm, &A2)); 495 PetscCall(MatSetType(A2, MATMPIAIJ)); 496 PetscCall(MatSetSizes(A2, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 497 PetscCall(MatMPIAIJSetPreallocation(A2, 0, dnz, 0, onz)); 498 PetscCall(MatSetOption(A2, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 499 for (i = 0, j = 0; i < A2->rmap->n; i++) { 500 PetscInt row = i + A2->rmap->rstart, s = dnz[i] + onz[i]; 501 const PetscInt *idxs; 502 503 if (s) { 504 PetscCall(ISGetIndices(workis[j], &idxs)); 505 PetscCall(MatSetValues(A2, 1, &row, s, idxs, vals, INSERT_VALUES)); 506 PetscCall(ISRestoreIndices(workis[j], &idxs)); 507 j++; 508 } 509 } 510 PetscCheck(j == cnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT, j, cnt); 511 PetscCall(PetscFree(vals)); 512 PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY)); 513 PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY)); 514 PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2)); 515 516 /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */ 517 for (i = 0, j = 0; i < nl; i++) 518 if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i]; 519 PetscCall(ISCreateGeneral(comm, j, lndmapi, PETSC_USE_POINTER, &is)); 520 PetscCall(MatMPIAIJGetLocalMatCondensed(A2, MAT_INITIAL_MATRIX, &is, NULL, &A3)); 521 PetscCall(ISDestroy(&is)); 522 PetscCall(MatDestroy(&A2)); 523 524 /* extend local to global map to include connected isolated separators */ 525 PetscCall(PetscObjectQuery((PetscObject)A3, "_petsc_GetLocalMatCondensed_iscol", (PetscObject *)&is)); 526 PetscCheck(is, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing column map"); 527 PetscCall(ISLocalToGlobalMappingCreateIS(is, &ll2g)); 528 PetscCall(MatGetRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg)); 529 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure"); 530 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ii[i], jj, PETSC_COPY_VALUES, &is)); 531 PetscCall(MatRestoreRowIJ(A3, 0, PETSC_FALSE, PETSC_FALSE, &i, &ii, &jj, &flg)); 532 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure"); 533 PetscCall(ISLocalToGlobalMappingApplyIS(ll2g, is, &is2)); 534 PetscCall(ISDestroy(&is)); 535 PetscCall(ISLocalToGlobalMappingDestroy(&ll2g)); 536 537 /* add new nodes to the local-to-global map */ 538 PetscCall(ISLocalToGlobalMappingDestroy(l2g)); 539 PetscCall(ISExpand(ndsub, is2, &is)); 540 PetscCall(ISDestroy(&is2)); 541 PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g)); 542 PetscCall(ISDestroy(&is)); 543 544 PetscCall(MatDestroy(&A3)); 545 PetscCall(PetscFree(lndmapi)); 546 MatPreallocateEnd(dnz, onz); 547 for (i = 0; i < gcnt; i++) PetscCall(ISDestroy(&workis[i])); 548 PetscCall(PetscFree(workis)); 549 } 550 PetscCall(ISRestoreIndices(ndmap, &ndmapi)); 551 PetscCall(PetscSFDestroy(&sf)); 552 PetscCall(PetscFree(ndmapc)); 553 PetscCall(ISDestroy(&ndmap)); 554 PetscCall(ISDestroy(&ndsub)); 555 PetscCall(ISLocalToGlobalMappingSetBlockSize(*l2g, bs)); 556 PetscCall(ISLocalToGlobalMappingViewFromOptions(*l2g, NULL, "-matis_nd_l2g_view")); 557 break; 558 case MAT_IS_DISASSEMBLE_L2G_NATURAL: 559 PetscCall(PetscObjectQuery((PetscObject)A, "__PETSc_dm", (PetscObject *)&dm)); 560 if (dm) { /* if a matrix comes from a DM, most likely we can use the l2gmap if any */ 561 PetscCall(MatGetLocalToGlobalMapping(A, l2g, NULL)); 562 PetscCall(PetscObjectReference((PetscObject)*l2g)); 563 if (*l2g) PetscFunctionReturn(0); 564 } 565 if (ismpiaij) { 566 PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray)); 567 } else if (ismpibaij) { 568 PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray)); 569 } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name); 570 if (A->rmap->n) { 571 PetscInt dc, oc, stc, *aux; 572 573 PetscCall(MatGetLocalSize(Ad, NULL, &dc)); 574 PetscCall(MatGetLocalSize(Ao, NULL, &oc)); 575 PetscCheck(!oc || garray, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "garray not present"); 576 PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL)); 577 PetscCall(PetscMalloc1((dc + oc) / bs, &aux)); 578 for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs; 579 for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = (ismpiaij ? garray[i * bs] / bs : garray[i]); 580 PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is)); 581 } else { 582 PetscCall(ISCreateBlock(comm, 1, 0, NULL, PETSC_OWN_POINTER, &is)); 583 } 584 PetscCall(ISLocalToGlobalMappingCreateIS(is, l2g)); 585 PetscCall(ISDestroy(&is)); 586 break; 587 default: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unsupported l2g disassembling type %d", mode); 588 } 589 PetscFunctionReturn(0); 590 } 591 592 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat) { 593 Mat lA, Ad, Ao, B = NULL; 594 ISLocalToGlobalMapping rl2g, cl2g; 595 IS is; 596 MPI_Comm comm; 597 void *ptrs[2]; 598 const char *names[2] = {"_convert_csr_aux", "_convert_csr_data"}; 599 const PetscInt *garray; 600 PetscScalar *dd, *od, *aa, *data; 601 const PetscInt *di, *dj, *oi, *oj; 602 const PetscInt *odi, *odj, *ooi, *ooj; 603 PetscInt *aux, *ii, *jj; 604 PetscInt bs, lc, dr, dc, oc, str, stc, nnz, i, jd, jo, cum; 605 PetscBool flg, ismpiaij, ismpibaij, was_inplace = PETSC_FALSE; 606 PetscMPIInt size; 607 608 PetscFunctionBegin; 609 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 610 PetscCallMPI(MPI_Comm_size(comm, &size)); 611 if (size == 1) { 612 PetscCall(MatConvert_SeqXAIJ_IS(A, type, reuse, newmat)); 613 PetscFunctionReturn(0); 614 } 615 if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) { 616 PetscCall(MatMPIXAIJComputeLocalToGlobalMapping_Private(A, &rl2g)); 617 PetscCall(MatCreate(comm, &B)); 618 PetscCall(MatSetType(B, MATIS)); 619 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 620 PetscCall(MatSetLocalToGlobalMapping(B, rl2g, rl2g)); 621 PetscCall(MatGetBlockSize(A, &bs)); 622 PetscCall(MatSetBlockSize(B, bs)); 623 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 624 if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE; 625 reuse = MAT_REUSE_MATRIX; 626 } 627 if (reuse == MAT_REUSE_MATRIX) { 628 Mat *newlA, lA; 629 IS rows, cols; 630 const PetscInt *ridx, *cidx; 631 PetscInt rbs, cbs, nr, nc; 632 633 if (!B) B = *newmat; 634 PetscCall(MatISGetLocalToGlobalMapping(B, &rl2g, &cl2g)); 635 PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g, &ridx)); 636 PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g, &cidx)); 637 PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &nr)); 638 PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &nc)); 639 PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &rbs)); 640 PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &cbs)); 641 PetscCall(ISCreateBlock(comm, rbs, nr / rbs, ridx, PETSC_USE_POINTER, &rows)); 642 if (rl2g != cl2g) { 643 PetscCall(ISCreateBlock(comm, cbs, nc / cbs, cidx, PETSC_USE_POINTER, &cols)); 644 } else { 645 PetscCall(PetscObjectReference((PetscObject)rows)); 646 cols = rows; 647 } 648 PetscCall(MatISGetLocalMat(B, &lA)); 649 PetscCall(MatCreateSubMatrices(A, 1, &rows, &cols, MAT_INITIAL_MATRIX, &newlA)); 650 PetscCall(MatConvert(newlA[0], MATSEQAIJ, MAT_INPLACE_MATRIX, &newlA[0])); 651 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g, &ridx)); 652 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g, &cidx)); 653 PetscCall(ISDestroy(&rows)); 654 PetscCall(ISDestroy(&cols)); 655 if (!lA->preallocated) { /* first time */ 656 PetscCall(MatDuplicate(newlA[0], MAT_COPY_VALUES, &lA)); 657 PetscCall(MatISSetLocalMat(B, lA)); 658 PetscCall(PetscObjectDereference((PetscObject)lA)); 659 } 660 PetscCall(MatCopy(newlA[0], lA, SAME_NONZERO_PATTERN)); 661 PetscCall(MatDestroySubMatrices(1, &newlA)); 662 PetscCall(MatISScaleDisassembling_Private(B)); 663 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 664 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 665 if (was_inplace) PetscCall(MatHeaderReplace(A, &B)); 666 else *newmat = B; 667 PetscFunctionReturn(0); 668 } 669 /* rectangular case, just compress out the column space */ 670 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij)); 671 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &ismpibaij)); 672 if (ismpiaij) { 673 bs = 1; 674 PetscCall(MatMPIAIJGetSeqAIJ(A, &Ad, &Ao, &garray)); 675 } else if (ismpibaij) { 676 PetscCall(MatGetBlockSize(A, &bs)); 677 PetscCall(MatMPIBAIJGetSeqBAIJ(A, &Ad, &Ao, &garray)); 678 PetscCall(MatConvert(Ad, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ad)); 679 PetscCall(MatConvert(Ao, MATSEQAIJ, MAT_INITIAL_MATRIX, &Ao)); 680 } else SETERRQ(comm, PETSC_ERR_SUP, "Type %s", ((PetscObject)A)->type_name); 681 PetscCall(MatSeqAIJGetArray(Ad, &dd)); 682 PetscCall(MatSeqAIJGetArray(Ao, &od)); 683 PetscCheck(garray, comm, PETSC_ERR_ARG_WRONGSTATE, "garray not present"); 684 685 /* access relevant information from MPIAIJ */ 686 PetscCall(MatGetOwnershipRange(A, &str, NULL)); 687 PetscCall(MatGetOwnershipRangeColumn(A, &stc, NULL)); 688 PetscCall(MatGetLocalSize(A, &dr, &dc)); 689 PetscCall(MatGetLocalSize(Ao, NULL, &oc)); 690 PetscCall(MatGetRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &di, &dj, &flg)); 691 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure"); 692 PetscCall(MatGetRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &oi, &oj, &flg)); 693 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get IJ structure"); 694 nnz = di[dr] + oi[dr]; 695 /* store original pointers to be restored later */ 696 odi = di; 697 odj = dj; 698 ooi = oi; 699 ooj = oj; 700 701 /* generate l2g maps for rows and cols */ 702 PetscCall(ISCreateStride(comm, dr / bs, str / bs, 1, &is)); 703 if (bs > 1) { 704 IS is2; 705 706 PetscCall(ISGetLocalSize(is, &i)); 707 PetscCall(ISGetIndices(is, (const PetscInt **)&aux)); 708 PetscCall(ISCreateBlock(comm, bs, i, aux, PETSC_COPY_VALUES, &is2)); 709 PetscCall(ISRestoreIndices(is, (const PetscInt **)&aux)); 710 PetscCall(ISDestroy(&is)); 711 is = is2; 712 } 713 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g)); 714 PetscCall(ISDestroy(&is)); 715 if (dr) { 716 PetscCall(PetscMalloc1((dc + oc) / bs, &aux)); 717 for (i = 0; i < dc / bs; i++) aux[i] = i + stc / bs; 718 for (i = 0; i < oc / bs; i++) aux[i + dc / bs] = garray[i]; 719 PetscCall(ISCreateBlock(comm, bs, (dc + oc) / bs, aux, PETSC_OWN_POINTER, &is)); 720 lc = dc + oc; 721 } else { 722 PetscCall(ISCreateBlock(comm, bs, 0, NULL, PETSC_OWN_POINTER, &is)); 723 lc = 0; 724 } 725 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g)); 726 PetscCall(ISDestroy(&is)); 727 728 /* create MATIS object */ 729 PetscCall(MatCreate(comm, &B)); 730 PetscCall(MatSetSizes(B, dr, dc, PETSC_DECIDE, PETSC_DECIDE)); 731 PetscCall(MatSetType(B, MATIS)); 732 PetscCall(MatSetBlockSize(B, bs)); 733 PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g)); 734 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 735 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 736 737 /* merge local matrices */ 738 PetscCall(PetscMalloc1(nnz + dr + 1, &aux)); 739 PetscCall(PetscMalloc1(nnz, &data)); 740 ii = aux; 741 jj = aux + dr + 1; 742 aa = data; 743 *ii = *(di++) + *(oi++); 744 for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) { 745 for (; jd < *di; jd++) { 746 *jj++ = *dj++; 747 *aa++ = *dd++; 748 } 749 for (; jo < *oi; jo++) { 750 *jj++ = *oj++ + dc; 751 *aa++ = *od++; 752 } 753 *(++ii) = *(di++) + *(oi++); 754 } 755 for (; cum < dr; cum++) *(++ii) = nnz; 756 757 PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_FALSE, PETSC_FALSE, &i, &odi, &odj, &flg)); 758 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure"); 759 PetscCall(MatRestoreRowIJ(Ao, 0, PETSC_FALSE, PETSC_FALSE, &i, &ooi, &ooj, &flg)); 760 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot restore IJ structure"); 761 PetscCall(MatSeqAIJRestoreArray(Ad, &dd)); 762 PetscCall(MatSeqAIJRestoreArray(Ao, &od)); 763 764 ii = aux; 765 jj = aux + dr + 1; 766 aa = data; 767 PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, lc, ii, jj, aa, &lA)); 768 769 /* create containers to destroy the data */ 770 ptrs[0] = aux; 771 ptrs[1] = data; 772 for (i = 0; i < 2; i++) { 773 PetscContainer c; 774 775 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c)); 776 PetscCall(PetscContainerSetPointer(c, ptrs[i])); 777 PetscCall(PetscContainerSetUserDestroy(c, PetscContainerUserDestroyDefault)); 778 PetscCall(PetscObjectCompose((PetscObject)lA, names[i], (PetscObject)c)); 779 PetscCall(PetscContainerDestroy(&c)); 780 } 781 if (ismpibaij) { /* destroy converted local matrices */ 782 PetscCall(MatDestroy(&Ad)); 783 PetscCall(MatDestroy(&Ao)); 784 } 785 786 /* finalize matrix */ 787 PetscCall(MatISSetLocalMat(B, lA)); 788 PetscCall(MatDestroy(&lA)); 789 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 790 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 791 if (reuse == MAT_INPLACE_MATRIX) { 792 PetscCall(MatHeaderReplace(A, &B)); 793 } else *newmat = B; 794 PetscFunctionReturn(0); 795 } 796 797 PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A, MatType type, MatReuse reuse, Mat *newmat) { 798 Mat **nest, *snest, **rnest, lA, B; 799 IS *iscol, *isrow, *islrow, *islcol; 800 ISLocalToGlobalMapping rl2g, cl2g; 801 MPI_Comm comm; 802 PetscInt *lr, *lc, *l2gidxs; 803 PetscInt i, j, nr, nc, rbs, cbs; 804 PetscBool convert, lreuse, *istrans; 805 806 PetscFunctionBegin; 807 PetscCall(MatNestGetSubMats(A, &nr, &nc, &nest)); 808 lreuse = PETSC_FALSE; 809 rnest = NULL; 810 if (reuse == MAT_REUSE_MATRIX) { 811 PetscBool ismatis, isnest; 812 813 PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis)); 814 PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)(*newmat))->type_name); 815 PetscCall(MatISGetLocalMat(*newmat, &lA)); 816 PetscCall(PetscObjectTypeCompare((PetscObject)lA, MATNEST, &isnest)); 817 if (isnest) { 818 PetscCall(MatNestGetSubMats(lA, &i, &j, &rnest)); 819 lreuse = (PetscBool)(i == nr && j == nc); 820 if (!lreuse) rnest = NULL; 821 } 822 } 823 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 824 PetscCall(PetscCalloc2(nr, &lr, nc, &lc)); 825 PetscCall(PetscCalloc6(nr, &isrow, nc, &iscol, nr, &islrow, nc, &islcol, nr * nc, &snest, nr * nc, &istrans)); 826 PetscCall(MatNestGetISs(A, isrow, iscol)); 827 for (i = 0; i < nr; i++) { 828 for (j = 0; j < nc; j++) { 829 PetscBool ismatis; 830 PetscInt l1, l2, lb1, lb2, ij = i * nc + j; 831 832 /* Null matrix pointers are allowed in MATNEST */ 833 if (!nest[i][j]) continue; 834 835 /* Nested matrices should be of type MATIS */ 836 PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATTRANSPOSEMAT, &istrans[ij])); 837 if (istrans[ij]) { 838 Mat T, lT; 839 PetscCall(MatTransposeGetMat(nest[i][j], &T)); 840 PetscCall(PetscObjectTypeCompare((PetscObject)T, MATIS, &ismatis)); 841 PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") (transposed) is not of type MATIS", i, j); 842 PetscCall(MatISGetLocalMat(T, &lT)); 843 PetscCall(MatCreateTranspose(lT, &snest[ij])); 844 } else { 845 PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j], MATIS, &ismatis)); 846 PetscCheck(ismatis, comm, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") is not of type MATIS", i, j); 847 PetscCall(MatISGetLocalMat(nest[i][j], &snest[ij])); 848 } 849 850 /* Check compatibility of local sizes */ 851 PetscCall(MatGetSize(snest[ij], &l1, &l2)); 852 PetscCall(MatGetBlockSizes(snest[ij], &lb1, &lb2)); 853 if (!l1 || !l2) continue; 854 PetscCheck(!lr[i] || l1 == lr[i], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lr[i], l1); 855 PetscCheck(!lc[j] || l2 == lc[j], PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, lc[j], l2); 856 lr[i] = l1; 857 lc[j] = l2; 858 859 /* check compatibilty for local matrix reusage */ 860 if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE; 861 } 862 } 863 864 if (PetscDefined(USE_DEBUG)) { 865 /* Check compatibility of l2g maps for rows */ 866 for (i = 0; i < nr; i++) { 867 rl2g = NULL; 868 for (j = 0; j < nc; j++) { 869 PetscInt n1, n2; 870 871 if (!nest[i][j]) continue; 872 if (istrans[i * nc + j]) { 873 Mat T; 874 875 PetscCall(MatTransposeGetMat(nest[i][j], &T)); 876 PetscCall(MatISGetLocalToGlobalMapping(T, NULL, &cl2g)); 877 } else { 878 PetscCall(MatISGetLocalToGlobalMapping(nest[i][j], &cl2g, NULL)); 879 } 880 PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1)); 881 if (!n1) continue; 882 if (!rl2g) { 883 rl2g = cl2g; 884 } else { 885 const PetscInt *idxs1, *idxs2; 886 PetscBool same; 887 888 PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2)); 889 PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, i, j, n1, n2); 890 PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1)); 891 PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2)); 892 PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same)); 893 PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1)); 894 PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2)); 895 PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap", i, j); 896 } 897 } 898 } 899 /* Check compatibility of l2g maps for columns */ 900 for (i = 0; i < nc; i++) { 901 rl2g = NULL; 902 for (j = 0; j < nr; j++) { 903 PetscInt n1, n2; 904 905 if (!nest[j][i]) continue; 906 if (istrans[j * nc + i]) { 907 Mat T; 908 909 PetscCall(MatTransposeGetMat(nest[j][i], &T)); 910 PetscCall(MatISGetLocalToGlobalMapping(T, &cl2g, NULL)); 911 } else { 912 PetscCall(MatISGetLocalToGlobalMapping(nest[j][i], NULL, &cl2g)); 913 } 914 PetscCall(ISLocalToGlobalMappingGetSize(cl2g, &n1)); 915 if (!n1) continue; 916 if (!rl2g) { 917 rl2g = cl2g; 918 } else { 919 const PetscInt *idxs1, *idxs2; 920 PetscBool same; 921 922 PetscCall(ISLocalToGlobalMappingGetSize(rl2g, &n2)); 923 PetscCheck(n1 == n2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT, j, i, n1, n2); 924 PetscCall(ISLocalToGlobalMappingGetIndices(cl2g, &idxs1)); 925 PetscCall(ISLocalToGlobalMappingGetIndices(rl2g, &idxs2)); 926 PetscCall(PetscArraycmp(idxs1, idxs2, n1, &same)); 927 PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g, &idxs1)); 928 PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g, &idxs2)); 929 PetscCheck(same, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap", j, i); 930 } 931 } 932 } 933 } 934 935 B = NULL; 936 if (reuse != MAT_REUSE_MATRIX) { 937 PetscInt stl; 938 939 /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */ 940 for (i = 0, stl = 0; i < nr; i++) stl += lr[i]; 941 PetscCall(PetscMalloc1(stl, &l2gidxs)); 942 for (i = 0, stl = 0; i < nr; i++) { 943 Mat usedmat; 944 Mat_IS *matis; 945 const PetscInt *idxs; 946 947 /* local IS for local NEST */ 948 PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i])); 949 950 /* l2gmap */ 951 j = 0; 952 usedmat = nest[i][j]; 953 while (!usedmat && j < nc - 1) usedmat = nest[i][++j]; 954 PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid row mat"); 955 956 if (istrans[i * nc + j]) { 957 Mat T; 958 PetscCall(MatTransposeGetMat(usedmat, &T)); 959 usedmat = T; 960 } 961 matis = (Mat_IS *)(usedmat->data); 962 PetscCall(ISGetIndices(isrow[i], &idxs)); 963 if (istrans[i * nc + j]) { 964 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE)); 965 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE)); 966 } else { 967 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE)); 968 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE)); 969 } 970 PetscCall(ISRestoreIndices(isrow[i], &idxs)); 971 stl += lr[i]; 972 } 973 PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &rl2g)); 974 975 /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */ 976 for (i = 0, stl = 0; i < nc; i++) stl += lc[i]; 977 PetscCall(PetscMalloc1(stl, &l2gidxs)); 978 for (i = 0, stl = 0; i < nc; i++) { 979 Mat usedmat; 980 Mat_IS *matis; 981 const PetscInt *idxs; 982 983 /* local IS for local NEST */ 984 PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i])); 985 986 /* l2gmap */ 987 j = 0; 988 usedmat = nest[j][i]; 989 while (!usedmat && j < nr - 1) usedmat = nest[++j][i]; 990 PetscCheck(usedmat, comm, PETSC_ERR_SUP, "Cannot find valid column mat"); 991 if (istrans[j * nc + i]) { 992 Mat T; 993 PetscCall(MatTransposeGetMat(usedmat, &T)); 994 usedmat = T; 995 } 996 matis = (Mat_IS *)(usedmat->data); 997 PetscCall(ISGetIndices(iscol[i], &idxs)); 998 if (istrans[j * nc + i]) { 999 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE)); 1000 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE)); 1001 } else { 1002 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE)); 1003 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, idxs, l2gidxs + stl, MPI_REPLACE)); 1004 } 1005 PetscCall(ISRestoreIndices(iscol[i], &idxs)); 1006 stl += lc[i]; 1007 } 1008 PetscCall(ISLocalToGlobalMappingCreate(comm, 1, stl, l2gidxs, PETSC_OWN_POINTER, &cl2g)); 1009 1010 /* Create MATIS */ 1011 PetscCall(MatCreate(comm, &B)); 1012 PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N)); 1013 PetscCall(MatGetBlockSizes(A, &rbs, &cbs)); 1014 PetscCall(MatSetBlockSizes(B, rbs, cbs)); 1015 PetscCall(MatSetType(B, MATIS)); 1016 PetscCall(MatISSetLocalMatType(B, MATNEST)); 1017 { /* hack : avoid setup of scatters */ 1018 Mat_IS *matis = (Mat_IS *)(B->data); 1019 matis->islocalref = PETSC_TRUE; 1020 } 1021 PetscCall(MatSetLocalToGlobalMapping(B, rl2g, cl2g)); 1022 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 1023 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 1024 PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA)); 1025 PetscCall(MatNestSetVecType(lA, VECNEST)); 1026 for (i = 0; i < nr * nc; i++) { 1027 if (istrans[i]) PetscCall(MatDestroy(&snest[i])); 1028 } 1029 PetscCall(MatISSetLocalMat(B, lA)); 1030 PetscCall(MatDestroy(&lA)); 1031 { /* hack : setup of scatters done here */ 1032 Mat_IS *matis = (Mat_IS *)(B->data); 1033 1034 matis->islocalref = PETSC_FALSE; 1035 PetscCall(MatISSetUpScatters_Private(B)); 1036 } 1037 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1038 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 1039 if (reuse == MAT_INPLACE_MATRIX) { 1040 PetscCall(MatHeaderReplace(A, &B)); 1041 } else { 1042 *newmat = B; 1043 } 1044 } else { 1045 if (lreuse) { 1046 PetscCall(MatISGetLocalMat(*newmat, &lA)); 1047 for (i = 0; i < nr; i++) { 1048 for (j = 0; j < nc; j++) { 1049 if (snest[i * nc + j]) { 1050 PetscCall(MatNestSetSubMat(lA, i, j, snest[i * nc + j])); 1051 if (istrans[i * nc + j]) PetscCall(MatDestroy(&snest[i * nc + j])); 1052 } 1053 } 1054 } 1055 } else { 1056 PetscInt stl; 1057 for (i = 0, stl = 0; i < nr; i++) { 1058 PetscCall(ISCreateStride(PETSC_COMM_SELF, lr[i], stl, 1, &islrow[i])); 1059 stl += lr[i]; 1060 } 1061 for (i = 0, stl = 0; i < nc; i++) { 1062 PetscCall(ISCreateStride(PETSC_COMM_SELF, lc[i], stl, 1, &islcol[i])); 1063 stl += lc[i]; 1064 } 1065 PetscCall(MatCreateNest(PETSC_COMM_SELF, nr, islrow, nc, islcol, snest, &lA)); 1066 for (i = 0; i < nr * nc; i++) { 1067 if (istrans[i]) PetscCall(MatDestroy(&snest[i])); 1068 } 1069 PetscCall(MatISSetLocalMat(*newmat, lA)); 1070 PetscCall(MatDestroy(&lA)); 1071 } 1072 PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY)); 1073 PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY)); 1074 } 1075 1076 /* Create local matrix in MATNEST format */ 1077 convert = PETSC_FALSE; 1078 PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-matis_convert_local_nest", &convert, NULL)); 1079 if (convert) { 1080 Mat M; 1081 MatISLocalFields lf; 1082 PetscContainer c; 1083 1084 PetscCall(MatISGetLocalMat(*newmat, &lA)); 1085 PetscCall(MatConvert(lA, MATAIJ, MAT_INITIAL_MATRIX, &M)); 1086 PetscCall(MatISSetLocalMat(*newmat, M)); 1087 PetscCall(MatDestroy(&M)); 1088 1089 /* attach local fields to the matrix */ 1090 PetscCall(PetscNew(&lf)); 1091 PetscCall(PetscMalloc2(nr, &lf->rf, nc, &lf->cf)); 1092 for (i = 0; i < nr; i++) { 1093 PetscInt n, st; 1094 1095 PetscCall(ISGetLocalSize(islrow[i], &n)); 1096 PetscCall(ISStrideGetInfo(islrow[i], &st, NULL)); 1097 PetscCall(ISCreateStride(comm, n, st, 1, &lf->rf[i])); 1098 } 1099 for (i = 0; i < nc; i++) { 1100 PetscInt n, st; 1101 1102 PetscCall(ISGetLocalSize(islcol[i], &n)); 1103 PetscCall(ISStrideGetInfo(islcol[i], &st, NULL)); 1104 PetscCall(ISCreateStride(comm, n, st, 1, &lf->cf[i])); 1105 } 1106 lf->nr = nr; 1107 lf->nc = nc; 1108 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)), &c)); 1109 PetscCall(PetscContainerSetPointer(c, lf)); 1110 PetscCall(PetscContainerSetUserDestroy(c, MatISContainerDestroyFields_Private)); 1111 PetscCall(PetscObjectCompose((PetscObject)(*newmat), "_convert_nest_lfields", (PetscObject)c)); 1112 PetscCall(PetscContainerDestroy(&c)); 1113 } 1114 1115 /* Free workspace */ 1116 for (i = 0; i < nr; i++) PetscCall(ISDestroy(&islrow[i])); 1117 for (i = 0; i < nc; i++) PetscCall(ISDestroy(&islcol[i])); 1118 PetscCall(PetscFree6(isrow, iscol, islrow, islcol, snest, istrans)); 1119 PetscCall(PetscFree2(lr, lc)); 1120 PetscFunctionReturn(0); 1121 } 1122 1123 static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r) { 1124 Mat_IS *matis = (Mat_IS *)A->data; 1125 Vec ll, rr; 1126 const PetscScalar *Y, *X; 1127 PetscScalar *x, *y; 1128 1129 PetscFunctionBegin; 1130 if (l) { 1131 ll = matis->y; 1132 PetscCall(VecGetArrayRead(l, &Y)); 1133 PetscCall(VecGetArray(ll, &y)); 1134 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE)); 1135 } else { 1136 ll = NULL; 1137 } 1138 if (r) { 1139 rr = matis->x; 1140 PetscCall(VecGetArrayRead(r, &X)); 1141 PetscCall(VecGetArray(rr, &x)); 1142 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE)); 1143 } else { 1144 rr = NULL; 1145 } 1146 if (ll) { 1147 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_SCALAR, Y, y, MPI_REPLACE)); 1148 PetscCall(VecRestoreArrayRead(l, &Y)); 1149 PetscCall(VecRestoreArray(ll, &y)); 1150 } 1151 if (rr) { 1152 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_SCALAR, X, x, MPI_REPLACE)); 1153 PetscCall(VecRestoreArrayRead(r, &X)); 1154 PetscCall(VecRestoreArray(rr, &x)); 1155 } 1156 PetscCall(MatDiagonalScale(matis->A, ll, rr)); 1157 PetscFunctionReturn(0); 1158 } 1159 1160 static PetscErrorCode MatGetInfo_IS(Mat A, MatInfoType flag, MatInfo *ginfo) { 1161 Mat_IS *matis = (Mat_IS *)A->data; 1162 MatInfo info; 1163 PetscLogDouble isend[6], irecv[6]; 1164 PetscInt bs; 1165 1166 PetscFunctionBegin; 1167 PetscCall(MatGetBlockSize(A, &bs)); 1168 if (matis->A->ops->getinfo) { 1169 PetscCall(MatGetInfo(matis->A, MAT_LOCAL, &info)); 1170 isend[0] = info.nz_used; 1171 isend[1] = info.nz_allocated; 1172 isend[2] = info.nz_unneeded; 1173 isend[3] = info.memory; 1174 isend[4] = info.mallocs; 1175 } else { 1176 isend[0] = 0.; 1177 isend[1] = 0.; 1178 isend[2] = 0.; 1179 isend[3] = 0.; 1180 isend[4] = 0.; 1181 } 1182 isend[5] = matis->A->num_ass; 1183 if (flag == MAT_LOCAL) { 1184 ginfo->nz_used = isend[0]; 1185 ginfo->nz_allocated = isend[1]; 1186 ginfo->nz_unneeded = isend[2]; 1187 ginfo->memory = isend[3]; 1188 ginfo->mallocs = isend[4]; 1189 ginfo->assemblies = isend[5]; 1190 } else if (flag == MAT_GLOBAL_MAX) { 1191 PetscCall(MPIU_Allreduce(isend, irecv, 6, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A))); 1192 1193 ginfo->nz_used = irecv[0]; 1194 ginfo->nz_allocated = irecv[1]; 1195 ginfo->nz_unneeded = irecv[2]; 1196 ginfo->memory = irecv[3]; 1197 ginfo->mallocs = irecv[4]; 1198 ginfo->assemblies = irecv[5]; 1199 } else if (flag == MAT_GLOBAL_SUM) { 1200 PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A))); 1201 1202 ginfo->nz_used = irecv[0]; 1203 ginfo->nz_allocated = irecv[1]; 1204 ginfo->nz_unneeded = irecv[2]; 1205 ginfo->memory = irecv[3]; 1206 ginfo->mallocs = irecv[4]; 1207 ginfo->assemblies = A->num_ass; 1208 } 1209 ginfo->block_size = bs; 1210 ginfo->fill_ratio_given = 0; 1211 ginfo->fill_ratio_needed = 0; 1212 ginfo->factor_mallocs = 0; 1213 PetscFunctionReturn(0); 1214 } 1215 1216 static PetscErrorCode MatTranspose_IS(Mat A, MatReuse reuse, Mat *B) { 1217 Mat C, lC, lA; 1218 1219 PetscFunctionBegin; 1220 if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B)); 1221 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 1222 ISLocalToGlobalMapping rl2g, cl2g; 1223 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); 1224 PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N)); 1225 PetscCall(MatSetBlockSizes(C, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs))); 1226 PetscCall(MatSetType(C, MATIS)); 1227 PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g)); 1228 PetscCall(MatSetLocalToGlobalMapping(C, cl2g, rl2g)); 1229 } else C = *B; 1230 1231 /* perform local transposition */ 1232 PetscCall(MatISGetLocalMat(A, &lA)); 1233 PetscCall(MatTranspose(lA, MAT_INITIAL_MATRIX, &lC)); 1234 PetscCall(MatSetLocalToGlobalMapping(lC, lA->cmap->mapping, lA->rmap->mapping)); 1235 PetscCall(MatISSetLocalMat(C, lC)); 1236 PetscCall(MatDestroy(&lC)); 1237 1238 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 1239 *B = C; 1240 } else { 1241 PetscCall(MatHeaderMerge(A, &C)); 1242 } 1243 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 1244 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 1245 PetscFunctionReturn(0); 1246 } 1247 1248 static PetscErrorCode MatDiagonalSet_IS(Mat A, Vec D, InsertMode insmode) { 1249 Mat_IS *is = (Mat_IS *)A->data; 1250 1251 PetscFunctionBegin; 1252 if (D) { /* MatShift_IS pass D = NULL */ 1253 PetscCall(VecScatterBegin(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD)); 1254 PetscCall(VecScatterEnd(is->rctx, D, is->y, INSERT_VALUES, SCATTER_FORWARD)); 1255 } 1256 PetscCall(VecPointwiseDivide(is->y, is->y, is->counter)); 1257 PetscCall(MatDiagonalSet(is->A, is->y, insmode)); 1258 PetscFunctionReturn(0); 1259 } 1260 1261 static PetscErrorCode MatShift_IS(Mat A, PetscScalar a) { 1262 Mat_IS *is = (Mat_IS *)A->data; 1263 1264 PetscFunctionBegin; 1265 PetscCall(VecSet(is->y, a)); 1266 PetscCall(MatDiagonalSet_IS(A, NULL, ADD_VALUES)); 1267 PetscFunctionReturn(0); 1268 } 1269 1270 static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) { 1271 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1272 1273 PetscFunctionBegin; 1274 PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of row/column indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT, MATIS_MAX_ENTRIES_INSERTION, m, n); 1275 PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, m, rows, rows_l)); 1276 PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping, n, cols, cols_l)); 1277 PetscCall(MatSetValuesLocal_IS(A, m, rows_l, n, cols_l, values, addv)); 1278 PetscFunctionReturn(0); 1279 } 1280 1281 static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) { 1282 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1283 1284 PetscFunctionBegin; 1285 PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION, PETSC_COMM_SELF, PETSC_ERR_SUP, "Number of row/column block indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT, MATIS_MAX_ENTRIES_INSERTION, m, n); 1286 PetscCall(ISLocalToGlobalMappingApplyBlock(A->rmap->mapping, m, rows, rows_l)); 1287 PetscCall(ISLocalToGlobalMappingApplyBlock(A->cmap->mapping, n, cols, cols_l)); 1288 PetscCall(MatSetValuesBlockedLocal_IS(A, m, rows_l, n, cols_l, values, addv)); 1289 PetscFunctionReturn(0); 1290 } 1291 1292 static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat) { 1293 Mat locmat, newlocmat; 1294 Mat_IS *newmatis; 1295 const PetscInt *idxs; 1296 PetscInt i, m, n; 1297 1298 PetscFunctionBegin; 1299 if (scall == MAT_REUSE_MATRIX) { 1300 PetscBool ismatis; 1301 1302 PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis)); 1303 PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Not of MATIS type"); 1304 newmatis = (Mat_IS *)(*newmat)->data; 1305 PetscCheck(newmatis->getsub_ris, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local row IS"); 1306 PetscCheck(newmatis->getsub_cis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local col IS"); 1307 } 1308 /* irow and icol may not have duplicate entries */ 1309 if (PetscDefined(USE_DEBUG)) { 1310 Vec rtest, ltest; 1311 const PetscScalar *array; 1312 1313 PetscCall(MatCreateVecs(mat, <est, &rtest)); 1314 PetscCall(ISGetLocalSize(irow, &n)); 1315 PetscCall(ISGetIndices(irow, &idxs)); 1316 for (i = 0; i < n; i++) PetscCall(VecSetValue(rtest, idxs[i], 1.0, ADD_VALUES)); 1317 PetscCall(VecAssemblyBegin(rtest)); 1318 PetscCall(VecAssemblyEnd(rtest)); 1319 PetscCall(VecGetLocalSize(rtest, &n)); 1320 PetscCall(VecGetOwnershipRange(rtest, &m, NULL)); 1321 PetscCall(VecGetArrayRead(rtest, &array)); 1322 for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Irow may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i])); 1323 PetscCall(VecRestoreArrayRead(rtest, &array)); 1324 PetscCall(ISRestoreIndices(irow, &idxs)); 1325 PetscCall(ISGetLocalSize(icol, &n)); 1326 PetscCall(ISGetIndices(icol, &idxs)); 1327 for (i = 0; i < n; i++) PetscCall(VecSetValue(ltest, idxs[i], 1.0, ADD_VALUES)); 1328 PetscCall(VecAssemblyBegin(ltest)); 1329 PetscCall(VecAssemblyEnd(ltest)); 1330 PetscCall(VecGetLocalSize(ltest, &n)); 1331 PetscCall(VecGetOwnershipRange(ltest, &m, NULL)); 1332 PetscCall(VecGetArrayRead(ltest, &array)); 1333 for (i = 0; i < n; i++) PetscCheck(array[i] == 0. || array[i] == 1., PETSC_COMM_SELF, PETSC_ERR_SUP, "Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Icol may not have duplicate entries", i + m, (PetscInt)PetscRealPart(array[i])); 1334 PetscCall(VecRestoreArrayRead(ltest, &array)); 1335 PetscCall(ISRestoreIndices(icol, &idxs)); 1336 PetscCall(VecDestroy(&rtest)); 1337 PetscCall(VecDestroy(<est)); 1338 } 1339 if (scall == MAT_INITIAL_MATRIX) { 1340 Mat_IS *matis = (Mat_IS *)mat->data; 1341 ISLocalToGlobalMapping rl2g; 1342 IS is; 1343 PetscInt *lidxs, *lgidxs, *newgidxs; 1344 PetscInt ll, newloc, irbs, icbs, arbs, acbs, rbs, cbs; 1345 PetscBool cong; 1346 MPI_Comm comm; 1347 1348 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 1349 PetscCall(MatGetBlockSizes(mat, &arbs, &acbs)); 1350 PetscCall(ISGetBlockSize(irow, &irbs)); 1351 PetscCall(ISGetBlockSize(icol, &icbs)); 1352 rbs = arbs == irbs ? irbs : 1; 1353 cbs = acbs == icbs ? icbs : 1; 1354 PetscCall(ISGetLocalSize(irow, &m)); 1355 PetscCall(ISGetLocalSize(icol, &n)); 1356 PetscCall(MatCreate(comm, newmat)); 1357 PetscCall(MatSetType(*newmat, MATIS)); 1358 PetscCall(MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE)); 1359 PetscCall(MatSetBlockSizes(*newmat, rbs, cbs)); 1360 /* communicate irow to their owners in the layout */ 1361 PetscCall(ISGetIndices(irow, &idxs)); 1362 PetscCall(PetscLayoutMapLocal(mat->rmap, m, idxs, &ll, &lidxs, &lgidxs)); 1363 PetscCall(ISRestoreIndices(irow, &idxs)); 1364 PetscCall(PetscArrayzero(matis->sf_rootdata, matis->sf->nroots)); 1365 for (i = 0; i < ll; i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i] + 1; 1366 PetscCall(PetscFree(lidxs)); 1367 PetscCall(PetscFree(lgidxs)); 1368 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 1369 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 1370 for (i = 0, newloc = 0; i < matis->sf->nleaves; i++) 1371 if (matis->sf_leafdata[i]) newloc++; 1372 PetscCall(PetscMalloc1(newloc, &newgidxs)); 1373 PetscCall(PetscMalloc1(newloc, &lidxs)); 1374 for (i = 0, newloc = 0; i < matis->sf->nleaves; i++) 1375 if (matis->sf_leafdata[i]) { 1376 lidxs[newloc] = i; 1377 newgidxs[newloc++] = matis->sf_leafdata[i] - 1; 1378 } 1379 PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is)); 1380 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g)); 1381 PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs)); 1382 PetscCall(ISDestroy(&is)); 1383 /* local is to extract local submatrix */ 1384 newmatis = (Mat_IS *)(*newmat)->data; 1385 PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_ris)); 1386 PetscCall(MatHasCongruentLayouts(mat, &cong)); 1387 if (cong && irow == icol && matis->csf == matis->sf) { 1388 PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, rl2g)); 1389 PetscCall(PetscObjectReference((PetscObject)newmatis->getsub_ris)); 1390 newmatis->getsub_cis = newmatis->getsub_ris; 1391 } else { 1392 ISLocalToGlobalMapping cl2g; 1393 1394 /* communicate icol to their owners in the layout */ 1395 PetscCall(ISGetIndices(icol, &idxs)); 1396 PetscCall(PetscLayoutMapLocal(mat->cmap, n, idxs, &ll, &lidxs, &lgidxs)); 1397 PetscCall(ISRestoreIndices(icol, &idxs)); 1398 PetscCall(PetscArrayzero(matis->csf_rootdata, matis->csf->nroots)); 1399 for (i = 0; i < ll; i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i] + 1; 1400 PetscCall(PetscFree(lidxs)); 1401 PetscCall(PetscFree(lgidxs)); 1402 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE)); 1403 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE)); 1404 for (i = 0, newloc = 0; i < matis->csf->nleaves; i++) 1405 if (matis->csf_leafdata[i]) newloc++; 1406 PetscCall(PetscMalloc1(newloc, &newgidxs)); 1407 PetscCall(PetscMalloc1(newloc, &lidxs)); 1408 for (i = 0, newloc = 0; i < matis->csf->nleaves; i++) 1409 if (matis->csf_leafdata[i]) { 1410 lidxs[newloc] = i; 1411 newgidxs[newloc++] = matis->csf_leafdata[i] - 1; 1412 } 1413 PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is)); 1414 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g)); 1415 PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs)); 1416 PetscCall(ISDestroy(&is)); 1417 /* local is to extract local submatrix */ 1418 PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_cis)); 1419 PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, cl2g)); 1420 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 1421 } 1422 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 1423 } else { 1424 PetscCall(MatISGetLocalMat(*newmat, &newlocmat)); 1425 } 1426 PetscCall(MatISGetLocalMat(mat, &locmat)); 1427 newmatis = (Mat_IS *)(*newmat)->data; 1428 PetscCall(MatCreateSubMatrix(locmat, newmatis->getsub_ris, newmatis->getsub_cis, scall, &newlocmat)); 1429 if (scall == MAT_INITIAL_MATRIX) { 1430 PetscCall(MatISSetLocalMat(*newmat, newlocmat)); 1431 PetscCall(MatDestroy(&newlocmat)); 1432 } 1433 PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY)); 1434 PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY)); 1435 PetscFunctionReturn(0); 1436 } 1437 1438 static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str) { 1439 Mat_IS *a = (Mat_IS *)A->data, *b; 1440 PetscBool ismatis; 1441 1442 PetscFunctionBegin; 1443 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis)); 1444 PetscCheck(ismatis, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Need to be implemented"); 1445 b = (Mat_IS *)B->data; 1446 PetscCall(MatCopy(a->A, b->A, str)); 1447 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 1448 PetscFunctionReturn(0); 1449 } 1450 1451 static PetscErrorCode MatMissingDiagonal_IS(Mat A, PetscBool *missing, PetscInt *d) { 1452 Vec v; 1453 const PetscScalar *array; 1454 PetscInt i, n; 1455 1456 PetscFunctionBegin; 1457 *missing = PETSC_FALSE; 1458 PetscCall(MatCreateVecs(A, NULL, &v)); 1459 PetscCall(MatGetDiagonal(A, v)); 1460 PetscCall(VecGetLocalSize(v, &n)); 1461 PetscCall(VecGetArrayRead(v, &array)); 1462 for (i = 0; i < n; i++) 1463 if (array[i] == 0.) break; 1464 PetscCall(VecRestoreArrayRead(v, &array)); 1465 PetscCall(VecDestroy(&v)); 1466 if (i != n) *missing = PETSC_TRUE; 1467 if (d) { 1468 *d = -1; 1469 if (*missing) { 1470 PetscInt rstart; 1471 PetscCall(MatGetOwnershipRange(A, &rstart, NULL)); 1472 *d = i + rstart; 1473 } 1474 } 1475 PetscFunctionReturn(0); 1476 } 1477 1478 static PetscErrorCode MatISSetUpSF_IS(Mat B) { 1479 Mat_IS *matis = (Mat_IS *)(B->data); 1480 const PetscInt *gidxs; 1481 PetscInt nleaves; 1482 1483 PetscFunctionBegin; 1484 if (matis->sf) PetscFunctionReturn(0); 1485 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf)); 1486 PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs)); 1487 PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves)); 1488 PetscCall(PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs)); 1489 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs)); 1490 PetscCall(PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata)); 1491 if (matis->rmapping != matis->cmapping) { /* setup SF for columns */ 1492 PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves)); 1493 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf)); 1494 PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs)); 1495 PetscCall(PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs)); 1496 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs)); 1497 PetscCall(PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata)); 1498 } else { 1499 matis->csf = matis->sf; 1500 matis->csf_leafdata = matis->sf_leafdata; 1501 matis->csf_rootdata = matis->sf_rootdata; 1502 } 1503 PetscFunctionReturn(0); 1504 } 1505 1506 /*@ 1507 MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP. 1508 1509 Collective 1510 1511 Input Parameters: 1512 + A - the matrix 1513 - store - the boolean flag 1514 1515 Level: advanced 1516 1517 Notes: 1518 1519 .seealso: `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()` 1520 @*/ 1521 PetscErrorCode MatISStoreL2L(Mat A, PetscBool store) { 1522 PetscFunctionBegin; 1523 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1524 PetscValidType(A, 1); 1525 PetscValidLogicalCollectiveBool(A, store, 2); 1526 PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store)); 1527 PetscFunctionReturn(0); 1528 } 1529 1530 static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store) { 1531 Mat_IS *matis = (Mat_IS *)(A->data); 1532 1533 PetscFunctionBegin; 1534 matis->storel2l = store; 1535 if (!store) PetscCall(PetscObjectCompose((PetscObject)(A), "_MatIS_PtAP_l2l", NULL)); 1536 PetscFunctionReturn(0); 1537 } 1538 1539 /*@ 1540 MatISFixLocalEmpty - Compress out zero local rows from the local matrices 1541 1542 Collective 1543 1544 Input Parameters: 1545 + A - the matrix 1546 - fix - the boolean flag 1547 1548 Level: advanced 1549 1550 Notes: When fix is true, new local matrices and l2g maps are generated during the final assembly process. 1551 1552 .seealso: `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY` 1553 @*/ 1554 PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix) { 1555 PetscFunctionBegin; 1556 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1557 PetscValidType(A, 1); 1558 PetscValidLogicalCollectiveBool(A, fix, 2); 1559 PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix)); 1560 PetscFunctionReturn(0); 1561 } 1562 1563 static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix) { 1564 Mat_IS *matis = (Mat_IS *)(A->data); 1565 1566 PetscFunctionBegin; 1567 matis->locempty = fix; 1568 PetscFunctionReturn(0); 1569 } 1570 1571 /*@ 1572 MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix. 1573 1574 Collective 1575 1576 Input Parameters: 1577 + B - the matrix 1578 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1579 (same value is used for all local rows) 1580 . d_nnz - array containing the number of nonzeros in the various rows of the 1581 DIAGONAL portion of the local submatrix (possibly different for each row) 1582 or NULL, if d_nz is used to specify the nonzero structure. 1583 The size of this array is equal to the number of local rows, i.e 'm'. 1584 For matrices that will be factored, you must leave room for (and set) 1585 the diagonal entry even if it is zero. 1586 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1587 submatrix (same value is used for all local rows). 1588 - o_nnz - array containing the number of nonzeros in the various rows of the 1589 OFF-DIAGONAL portion of the local submatrix (possibly different for 1590 each row) or NULL, if o_nz is used to specify the nonzero 1591 structure. The size of this array is equal to the number 1592 of local rows, i.e 'm'. 1593 1594 If the *_nnz parameter is given then the *_nz parameter is ignored 1595 1596 Level: intermediate 1597 1598 Notes: 1599 This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition 1600 from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local 1601 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects. 1602 1603 .seealso: `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS` 1604 @*/ 1605 PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[]) { 1606 PetscFunctionBegin; 1607 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1608 PetscValidType(B, 1); 1609 PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz)); 1610 PetscFunctionReturn(0); 1611 } 1612 1613 /* this is used by DMDA */ 1614 PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[]) { 1615 Mat_IS *matis = (Mat_IS *)(B->data); 1616 PetscInt bs, i, nlocalcols; 1617 1618 PetscFunctionBegin; 1619 PetscCall(MatSetUp(B)); 1620 if (!d_nnz) 1621 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz; 1622 else 1623 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i]; 1624 1625 if (!o_nnz) 1626 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz; 1627 else 1628 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i]; 1629 1630 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 1631 PetscCall(MatGetSize(matis->A, NULL, &nlocalcols)); 1632 PetscCall(MatGetBlockSize(matis->A, &bs)); 1633 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 1634 1635 for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols); 1636 PetscCall(MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata)); 1637 #if defined(PETSC_HAVE_HYPRE) 1638 PetscCall(MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL)); 1639 #endif 1640 1641 for (i = 0; i < matis->sf->nleaves / bs; i++) { 1642 PetscInt b; 1643 1644 matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs; 1645 for (b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs); 1646 } 1647 PetscCall(MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata)); 1648 1649 nlocalcols /= bs; 1650 for (i = 0; i < matis->sf->nleaves / bs; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols - i); 1651 PetscCall(MatSeqSBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata)); 1652 1653 /* for other matrix types */ 1654 PetscCall(MatSetUp(matis->A)); 1655 PetscFunctionReturn(0); 1656 } 1657 1658 PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce) { 1659 Mat_IS *matis = (Mat_IS *)(A->data); 1660 PetscInt *my_dnz, *my_onz, *dnz, *onz, *mat_ranges, *row_ownership; 1661 const PetscInt *global_indices_r, *global_indices_c; 1662 PetscInt i, j, bs, rows, cols; 1663 PetscInt lrows, lcols; 1664 PetscInt local_rows, local_cols; 1665 PetscMPIInt size; 1666 PetscBool isdense, issbaij; 1667 1668 PetscFunctionBegin; 1669 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 1670 PetscCall(MatGetSize(A, &rows, &cols)); 1671 PetscCall(MatGetBlockSize(A, &bs)); 1672 PetscCall(MatGetSize(matis->A, &local_rows, &local_cols)); 1673 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isdense)); 1674 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &issbaij)); 1675 PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &global_indices_r)); 1676 if (matis->rmapping != matis->cmapping) { 1677 PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &global_indices_c)); 1678 } else global_indices_c = global_indices_r; 1679 1680 if (issbaij) PetscCall(MatGetRowUpperTriangular(matis->A)); 1681 /* 1682 An SF reduce is needed to sum up properly on shared rows. 1683 Note that generally preallocation is not exact, since it overestimates nonzeros 1684 */ 1685 PetscCall(MatGetLocalSize(A, &lrows, &lcols)); 1686 MatPreallocateBegin(PetscObjectComm((PetscObject)A), lrows, lcols, dnz, onz); 1687 /* All processes need to compute entire row ownership */ 1688 PetscCall(PetscMalloc1(rows, &row_ownership)); 1689 PetscCall(MatGetOwnershipRanges(A, (const PetscInt **)&mat_ranges)); 1690 for (i = 0; i < size; i++) { 1691 for (j = mat_ranges[i]; j < mat_ranges[i + 1]; j++) row_ownership[j] = i; 1692 } 1693 PetscCall(MatGetOwnershipRangesColumn(A, (const PetscInt **)&mat_ranges)); 1694 1695 /* 1696 my_dnz and my_onz contains exact contribution to preallocation from each local mat 1697 then, they will be summed up properly. This way, preallocation is always sufficient 1698 */ 1699 PetscCall(PetscCalloc2(local_rows, &my_dnz, local_rows, &my_onz)); 1700 /* preallocation as a MATAIJ */ 1701 if (isdense) { /* special case for dense local matrices */ 1702 for (i = 0; i < local_rows; i++) { 1703 PetscInt owner = row_ownership[global_indices_r[i]]; 1704 for (j = 0; j < local_cols; j++) { 1705 PetscInt index_col = global_indices_c[j]; 1706 if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */ 1707 my_dnz[i] += 1; 1708 } else { /* offdiag block */ 1709 my_onz[i] += 1; 1710 } 1711 } 1712 } 1713 } else if (matis->A->ops->getrowij) { 1714 const PetscInt *ii, *jj, *jptr; 1715 PetscBool done; 1716 PetscCall(MatGetRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done)); 1717 PetscCheck(done, PetscObjectComm((PetscObject)(matis->A)), PETSC_ERR_PLIB, "Error in MatGetRowIJ"); 1718 jptr = jj; 1719 for (i = 0; i < local_rows; i++) { 1720 PetscInt index_row = global_indices_r[i]; 1721 for (j = 0; j < ii[i + 1] - ii[i]; j++, jptr++) { 1722 PetscInt owner = row_ownership[index_row]; 1723 PetscInt index_col = global_indices_c[*jptr]; 1724 if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */ 1725 my_dnz[i] += 1; 1726 } else { /* offdiag block */ 1727 my_onz[i] += 1; 1728 } 1729 /* same as before, interchanging rows and cols */ 1730 if (issbaij && index_col != index_row) { 1731 owner = row_ownership[index_col]; 1732 if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) { 1733 my_dnz[*jptr] += 1; 1734 } else { 1735 my_onz[*jptr] += 1; 1736 } 1737 } 1738 } 1739 } 1740 PetscCall(MatRestoreRowIJ(matis->A, 0, PETSC_FALSE, PETSC_FALSE, &local_rows, &ii, &jj, &done)); 1741 PetscCheck(done, PetscObjectComm((PetscObject)(matis->A)), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ"); 1742 } else { /* loop over rows and use MatGetRow */ 1743 for (i = 0; i < local_rows; i++) { 1744 const PetscInt *cols; 1745 PetscInt ncols, index_row = global_indices_r[i]; 1746 PetscCall(MatGetRow(matis->A, i, &ncols, &cols, NULL)); 1747 for (j = 0; j < ncols; j++) { 1748 PetscInt owner = row_ownership[index_row]; 1749 PetscInt index_col = global_indices_c[cols[j]]; 1750 if (index_col > mat_ranges[owner] - 1 && index_col < mat_ranges[owner + 1]) { /* diag block */ 1751 my_dnz[i] += 1; 1752 } else { /* offdiag block */ 1753 my_onz[i] += 1; 1754 } 1755 /* same as before, interchanging rows and cols */ 1756 if (issbaij && index_col != index_row) { 1757 owner = row_ownership[index_col]; 1758 if (index_row > mat_ranges[owner] - 1 && index_row < mat_ranges[owner + 1]) { 1759 my_dnz[cols[j]] += 1; 1760 } else { 1761 my_onz[cols[j]] += 1; 1762 } 1763 } 1764 } 1765 PetscCall(MatRestoreRow(matis->A, i, &ncols, &cols, NULL)); 1766 } 1767 } 1768 if (global_indices_c != global_indices_r) PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &global_indices_c)); 1769 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &global_indices_r)); 1770 PetscCall(PetscFree(row_ownership)); 1771 1772 /* Reduce my_dnz and my_onz */ 1773 if (maxreduce) { 1774 PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX)); 1775 PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX)); 1776 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_MAX)); 1777 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_MAX)); 1778 } else { 1779 PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM)); 1780 PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM)); 1781 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_dnz, dnz, MPI_SUM)); 1782 PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, my_onz, onz, MPI_SUM)); 1783 } 1784 PetscCall(PetscFree2(my_dnz, my_onz)); 1785 1786 /* Resize preallocation if overestimated */ 1787 for (i = 0; i < lrows; i++) { 1788 dnz[i] = PetscMin(dnz[i], lcols); 1789 onz[i] = PetscMin(onz[i], cols - lcols); 1790 } 1791 1792 /* Set preallocation */ 1793 PetscCall(MatSetBlockSizesFromMats(B, A, A)); 1794 PetscCall(MatSeqAIJSetPreallocation(B, 0, dnz)); 1795 PetscCall(MatMPIAIJSetPreallocation(B, 0, dnz, 0, onz)); 1796 for (i = 0; i < lrows; i += bs) { 1797 PetscInt b, d = dnz[i], o = onz[i]; 1798 1799 for (b = 1; b < bs; b++) { 1800 d = PetscMax(d, dnz[i + b]); 1801 o = PetscMax(o, onz[i + b]); 1802 } 1803 dnz[i / bs] = PetscMin(d / bs + d % bs, lcols / bs); 1804 onz[i / bs] = PetscMin(o / bs + o % bs, (cols - lcols) / bs); 1805 } 1806 PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, dnz)); 1807 PetscCall(MatMPIBAIJSetPreallocation(B, bs, 0, dnz, 0, onz)); 1808 PetscCall(MatMPISBAIJSetPreallocation(B, bs, 0, dnz, 0, onz)); 1809 MatPreallocateEnd(dnz, onz); 1810 if (issbaij) PetscCall(MatRestoreRowUpperTriangular(matis->A)); 1811 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); 1812 PetscFunctionReturn(0); 1813 } 1814 1815 PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M) { 1816 Mat_IS *matis = (Mat_IS *)(mat->data); 1817 Mat local_mat, MT; 1818 PetscInt rbs, cbs, rows, cols, lrows, lcols; 1819 PetscInt local_rows, local_cols; 1820 PetscBool isseqdense, isseqsbaij, isseqaij, isseqbaij; 1821 PetscMPIInt size; 1822 const PetscScalar *array; 1823 1824 PetscFunctionBegin; 1825 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size)); 1826 if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) { 1827 Mat B; 1828 IS irows = NULL, icols = NULL; 1829 PetscInt rbs, cbs; 1830 1831 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs)); 1832 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs)); 1833 if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */ 1834 IS rows, cols; 1835 const PetscInt *ridxs, *cidxs; 1836 PetscInt i, nw, *work; 1837 1838 PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping, &ridxs)); 1839 PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nw)); 1840 nw = nw / rbs; 1841 PetscCall(PetscCalloc1(nw, &work)); 1842 for (i = 0; i < nw; i++) work[ridxs[i]] += 1; 1843 for (i = 0; i < nw; i++) 1844 if (!work[i] || work[i] > 1) break; 1845 if (i == nw) { 1846 PetscCall(ISCreateBlock(PETSC_COMM_SELF, rbs, nw, ridxs, PETSC_USE_POINTER, &rows)); 1847 PetscCall(ISSetPermutation(rows)); 1848 PetscCall(ISInvertPermutation(rows, PETSC_DECIDE, &irows)); 1849 PetscCall(ISDestroy(&rows)); 1850 } 1851 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping, &ridxs)); 1852 PetscCall(PetscFree(work)); 1853 if (irows && matis->rmapping != matis->cmapping) { 1854 PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping, &cidxs)); 1855 PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nw)); 1856 nw = nw / cbs; 1857 PetscCall(PetscCalloc1(nw, &work)); 1858 for (i = 0; i < nw; i++) work[cidxs[i]] += 1; 1859 for (i = 0; i < nw; i++) 1860 if (!work[i] || work[i] > 1) break; 1861 if (i == nw) { 1862 PetscCall(ISCreateBlock(PETSC_COMM_SELF, cbs, nw, cidxs, PETSC_USE_POINTER, &cols)); 1863 PetscCall(ISSetPermutation(cols)); 1864 PetscCall(ISInvertPermutation(cols, PETSC_DECIDE, &icols)); 1865 PetscCall(ISDestroy(&cols)); 1866 } 1867 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping, &cidxs)); 1868 PetscCall(PetscFree(work)); 1869 } else if (irows) { 1870 PetscCall(PetscObjectReference((PetscObject)irows)); 1871 icols = irows; 1872 } 1873 } else { 1874 PetscCall(PetscObjectQuery((PetscObject)(*M), "_MatIS_IS_XAIJ_irows", (PetscObject *)&irows)); 1875 PetscCall(PetscObjectQuery((PetscObject)(*M), "_MatIS_IS_XAIJ_icols", (PetscObject *)&icols)); 1876 if (irows) PetscCall(PetscObjectReference((PetscObject)irows)); 1877 if (icols) PetscCall(PetscObjectReference((PetscObject)icols)); 1878 } 1879 if (!irows || !icols) { 1880 PetscCall(ISDestroy(&icols)); 1881 PetscCall(ISDestroy(&irows)); 1882 goto general_assembly; 1883 } 1884 PetscCall(MatConvert(matis->A, mtype, MAT_INITIAL_MATRIX, &B)); 1885 if (reuse != MAT_INPLACE_MATRIX) { 1886 PetscCall(MatCreateSubMatrix(B, irows, icols, reuse, M)); 1887 PetscCall(PetscObjectCompose((PetscObject)(*M), "_MatIS_IS_XAIJ_irows", (PetscObject)irows)); 1888 PetscCall(PetscObjectCompose((PetscObject)(*M), "_MatIS_IS_XAIJ_icols", (PetscObject)icols)); 1889 } else { 1890 Mat C; 1891 1892 PetscCall(MatCreateSubMatrix(B, irows, icols, MAT_INITIAL_MATRIX, &C)); 1893 PetscCall(MatHeaderReplace(mat, &C)); 1894 } 1895 PetscCall(MatDestroy(&B)); 1896 PetscCall(ISDestroy(&icols)); 1897 PetscCall(ISDestroy(&irows)); 1898 PetscFunctionReturn(0); 1899 } 1900 general_assembly: 1901 PetscCall(MatGetSize(mat, &rows, &cols)); 1902 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs)); 1903 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs)); 1904 PetscCall(MatGetLocalSize(mat, &lrows, &lcols)); 1905 PetscCall(MatGetSize(matis->A, &local_rows, &local_cols)); 1906 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isseqdense)); 1907 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij)); 1908 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQBAIJ, &isseqbaij)); 1909 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &isseqsbaij)); 1910 PetscCheck(isseqdense || isseqaij || isseqbaij || isseqsbaij, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)(matis->A))->type_name); 1911 if (PetscDefined(USE_DEBUG)) { 1912 PetscBool lb[4], bb[4]; 1913 1914 lb[0] = isseqdense; 1915 lb[1] = isseqaij; 1916 lb[2] = isseqbaij; 1917 lb[3] = isseqsbaij; 1918 PetscCall(MPIU_Allreduce(lb, bb, 4, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat))); 1919 PetscCheck(bb[0] || bb[1] || bb[2] || bb[3], PETSC_COMM_SELF, PETSC_ERR_SUP, "Local matrices must have the same type"); 1920 } 1921 1922 if (reuse != MAT_REUSE_MATRIX) { 1923 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &MT)); 1924 PetscCall(MatSetSizes(MT, lrows, lcols, rows, cols)); 1925 PetscCall(MatSetType(MT, mtype)); 1926 PetscCall(MatSetBlockSizes(MT, rbs, cbs)); 1927 PetscCall(MatISSetMPIXAIJPreallocation_Private(mat, MT, PETSC_FALSE)); 1928 } else { 1929 PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols; 1930 1931 /* some checks */ 1932 MT = *M; 1933 PetscCall(MatGetBlockSizes(MT, &mrbs, &mcbs)); 1934 PetscCall(MatGetSize(MT, &mrows, &mcols)); 1935 PetscCall(MatGetLocalSize(MT, &mlrows, &mlcols)); 1936 PetscCheck(mrows == rows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", rows, mrows); 1937 PetscCheck(mcols == cols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", cols, mcols); 1938 PetscCheck(mlrows == lrows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", lrows, mlrows); 1939 PetscCheck(mlcols == lcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", lcols, mlcols); 1940 PetscCheck(mrbs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", rbs, mrbs); 1941 PetscCheck(mcbs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", cbs, mcbs); 1942 PetscCall(MatZeroEntries(MT)); 1943 } 1944 1945 if (isseqsbaij || isseqbaij) { 1946 PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat)); 1947 isseqaij = PETSC_TRUE; 1948 } else { 1949 PetscCall(PetscObjectReference((PetscObject)matis->A)); 1950 local_mat = matis->A; 1951 } 1952 1953 /* Set values */ 1954 PetscCall(MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping)); 1955 if (isseqdense) { /* special case for dense local matrices */ 1956 PetscInt i, *dummy; 1957 1958 PetscCall(PetscMalloc1(PetscMax(local_rows, local_cols), &dummy)); 1959 for (i = 0; i < PetscMax(local_rows, local_cols); i++) dummy[i] = i; 1960 PetscCall(MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_FALSE)); 1961 PetscCall(MatDenseGetArrayRead(local_mat, &array)); 1962 PetscCall(MatSetValuesLocal(MT, local_rows, dummy, local_cols, dummy, array, ADD_VALUES)); 1963 PetscCall(MatDenseRestoreArrayRead(local_mat, &array)); 1964 PetscCall(PetscFree(dummy)); 1965 } else if (isseqaij) { 1966 const PetscInt *blocks; 1967 PetscInt i, nvtxs, *xadj, *adjncy, nb; 1968 PetscBool done; 1969 PetscScalar *sarray; 1970 1971 PetscCall(MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done)); 1972 PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatGetRowIJ"); 1973 PetscCall(MatSeqAIJGetArray(local_mat, &sarray)); 1974 PetscCall(MatGetVariableBlockSizes(local_mat, &nb, &blocks)); 1975 if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */ 1976 PetscInt sum; 1977 1978 for (i = 0, sum = 0; i < nb; i++) sum += blocks[i]; 1979 if (sum == nvtxs) { 1980 PetscInt r; 1981 1982 for (i = 0, r = 0; i < nb; i++) { 1983 PetscAssert(blocks[i] == xadj[r + 1] - xadj[r], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid block sizes prescribed for block %" PetscInt_FMT ": expected %" PetscInt_FMT ", got %" PetscInt_FMT, i, blocks[i], xadj[r + 1] - xadj[r]); 1984 PetscCall(MatSetValuesLocal(MT, blocks[i], adjncy + xadj[r], blocks[i], adjncy + xadj[r], sarray + xadj[r], ADD_VALUES)); 1985 r += blocks[i]; 1986 } 1987 } else { 1988 for (i = 0; i < nvtxs; i++) PetscCall(MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES)); 1989 } 1990 } else { 1991 for (i = 0; i < nvtxs; i++) PetscCall(MatSetValuesLocal(MT, 1, &i, xadj[i + 1] - xadj[i], adjncy + xadj[i], sarray + xadj[i], ADD_VALUES)); 1992 } 1993 PetscCall(MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nvtxs, (const PetscInt **)&xadj, (const PetscInt **)&adjncy, &done)); 1994 PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ"); 1995 PetscCall(MatSeqAIJRestoreArray(local_mat, &sarray)); 1996 } else { /* very basic values insertion for all other matrix types */ 1997 PetscInt i; 1998 1999 for (i = 0; i < local_rows; i++) { 2000 PetscInt j; 2001 const PetscInt *local_indices_cols; 2002 2003 PetscCall(MatGetRow(local_mat, i, &j, &local_indices_cols, &array)); 2004 PetscCall(MatSetValuesLocal(MT, 1, &i, j, local_indices_cols, array, ADD_VALUES)); 2005 PetscCall(MatRestoreRow(local_mat, i, &j, &local_indices_cols, &array)); 2006 } 2007 } 2008 PetscCall(MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY)); 2009 PetscCall(MatDestroy(&local_mat)); 2010 PetscCall(MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY)); 2011 if (isseqdense) PetscCall(MatSetOption(MT, MAT_ROW_ORIENTED, PETSC_TRUE)); 2012 if (reuse == MAT_INPLACE_MATRIX) { 2013 PetscCall(MatHeaderReplace(mat, &MT)); 2014 } else if (reuse == MAT_INITIAL_MATRIX) { 2015 *M = MT; 2016 } 2017 PetscFunctionReturn(0); 2018 } 2019 2020 /*@ 2021 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 2022 2023 Input Parameters: 2024 + mat - the matrix (should be of type MATIS) 2025 - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 2026 2027 Output Parameter: 2028 . newmat - the matrix in AIJ format 2029 2030 Level: developer 2031 2032 Notes: 2033 This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface. 2034 2035 .seealso: `MATIS`, `MatConvert()` 2036 @*/ 2037 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) { 2038 PetscFunctionBegin; 2039 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2040 PetscValidLogicalCollectiveEnum(mat, reuse, 2); 2041 PetscValidPointer(newmat, 3); 2042 if (reuse == MAT_REUSE_MATRIX) { 2043 PetscValidHeaderSpecific(*newmat, MAT_CLASSID, 3); 2044 PetscCheckSameComm(mat, 1, *newmat, 3); 2045 PetscCheck(mat != *newmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse the same matrix"); 2046 } 2047 PetscUseMethod(mat, "MatISGetMPIXAIJ_C", (Mat, MatType, MatReuse, Mat *), (mat, MATAIJ, reuse, newmat)); 2048 PetscFunctionReturn(0); 2049 } 2050 2051 static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat) { 2052 Mat_IS *matis = (Mat_IS *)(mat->data); 2053 PetscInt rbs, cbs, m, n, M, N; 2054 Mat B, localmat; 2055 2056 PetscFunctionBegin; 2057 PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs)); 2058 PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs)); 2059 PetscCall(MatGetSize(mat, &M, &N)); 2060 PetscCall(MatGetLocalSize(mat, &m, &n)); 2061 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &B)); 2062 PetscCall(MatSetSizes(B, m, n, M, N)); 2063 PetscCall(MatSetBlockSize(B, rbs == cbs ? rbs : 1)); 2064 PetscCall(MatSetType(B, MATIS)); 2065 PetscCall(MatISSetLocalMatType(B, matis->lmattype)); 2066 PetscCall(MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping)); 2067 PetscCall(MatDuplicate(matis->A, op, &localmat)); 2068 PetscCall(MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping)); 2069 PetscCall(MatISSetLocalMat(B, localmat)); 2070 PetscCall(MatDestroy(&localmat)); 2071 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 2072 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 2073 *newmat = B; 2074 PetscFunctionReturn(0); 2075 } 2076 2077 static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg) { 2078 Mat_IS *matis = (Mat_IS *)A->data; 2079 PetscBool local_sym; 2080 2081 PetscFunctionBegin; 2082 PetscCall(MatIsHermitian(matis->A, tol, &local_sym)); 2083 PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 2084 PetscFunctionReturn(0); 2085 } 2086 2087 static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg) { 2088 Mat_IS *matis = (Mat_IS *)A->data; 2089 PetscBool local_sym; 2090 2091 PetscFunctionBegin; 2092 if (matis->rmapping != matis->cmapping) { 2093 *flg = PETSC_FALSE; 2094 PetscFunctionReturn(0); 2095 } 2096 PetscCall(MatIsSymmetric(matis->A, tol, &local_sym)); 2097 PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 2098 PetscFunctionReturn(0); 2099 } 2100 2101 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg) { 2102 Mat_IS *matis = (Mat_IS *)A->data; 2103 PetscBool local_sym; 2104 2105 PetscFunctionBegin; 2106 if (matis->rmapping != matis->cmapping) { 2107 *flg = PETSC_FALSE; 2108 PetscFunctionReturn(0); 2109 } 2110 PetscCall(MatIsStructurallySymmetric(matis->A, &local_sym)); 2111 PetscCall(MPIU_Allreduce(&local_sym, flg, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 2112 PetscFunctionReturn(0); 2113 } 2114 2115 static PetscErrorCode MatDestroy_IS(Mat A) { 2116 Mat_IS *b = (Mat_IS *)A->data; 2117 2118 PetscFunctionBegin; 2119 PetscCall(PetscFree(b->bdiag)); 2120 PetscCall(PetscFree(b->lmattype)); 2121 PetscCall(MatDestroy(&b->A)); 2122 PetscCall(VecScatterDestroy(&b->cctx)); 2123 PetscCall(VecScatterDestroy(&b->rctx)); 2124 PetscCall(VecDestroy(&b->x)); 2125 PetscCall(VecDestroy(&b->y)); 2126 PetscCall(VecDestroy(&b->counter)); 2127 PetscCall(ISDestroy(&b->getsub_ris)); 2128 PetscCall(ISDestroy(&b->getsub_cis)); 2129 if (b->sf != b->csf) { 2130 PetscCall(PetscSFDestroy(&b->csf)); 2131 PetscCall(PetscFree2(b->csf_rootdata, b->csf_leafdata)); 2132 } else b->csf = NULL; 2133 PetscCall(PetscSFDestroy(&b->sf)); 2134 PetscCall(PetscFree2(b->sf_rootdata, b->sf_leafdata)); 2135 PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping)); 2136 PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping)); 2137 PetscCall(MatDestroy(&b->dA)); 2138 PetscCall(MatDestroy(&b->assembledA)); 2139 PetscCall(PetscFree(A->data)); 2140 PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL)); 2141 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL)); 2142 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL)); 2143 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL)); 2144 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL)); 2145 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetMPIXAIJ_C", NULL)); 2146 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL)); 2147 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL)); 2148 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL)); 2149 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL)); 2150 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL)); 2151 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL)); 2152 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL)); 2153 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL)); 2154 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL)); 2155 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL)); 2156 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL)); 2157 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL)); 2158 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL)); 2159 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL)); 2160 PetscFunctionReturn(0); 2161 } 2162 2163 static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y) { 2164 Mat_IS *is = (Mat_IS *)A->data; 2165 PetscScalar zero = 0.0; 2166 2167 PetscFunctionBegin; 2168 /* scatter the global vector x into the local work vector */ 2169 PetscCall(VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD)); 2170 PetscCall(VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD)); 2171 2172 /* multiply the local matrix */ 2173 PetscCall(MatMult(is->A, is->x, is->y)); 2174 2175 /* scatter product back into global memory */ 2176 PetscCall(VecSet(y, zero)); 2177 PetscCall(VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE)); 2178 PetscCall(VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE)); 2179 PetscFunctionReturn(0); 2180 } 2181 2182 static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3) { 2183 Vec temp_vec; 2184 2185 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 2186 if (v3 != v2) { 2187 PetscCall(MatMult(A, v1, v3)); 2188 PetscCall(VecAXPY(v3, 1.0, v2)); 2189 } else { 2190 PetscCall(VecDuplicate(v2, &temp_vec)); 2191 PetscCall(MatMult(A, v1, temp_vec)); 2192 PetscCall(VecAXPY(temp_vec, 1.0, v2)); 2193 PetscCall(VecCopy(temp_vec, v3)); 2194 PetscCall(VecDestroy(&temp_vec)); 2195 } 2196 PetscFunctionReturn(0); 2197 } 2198 2199 static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x) { 2200 Mat_IS *is = (Mat_IS *)A->data; 2201 2202 PetscFunctionBegin; 2203 /* scatter the global vector x into the local work vector */ 2204 PetscCall(VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD)); 2205 PetscCall(VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD)); 2206 2207 /* multiply the local matrix */ 2208 PetscCall(MatMultTranspose(is->A, is->y, is->x)); 2209 2210 /* scatter product back into global vector */ 2211 PetscCall(VecSet(x, 0)); 2212 PetscCall(VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE)); 2213 PetscCall(VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE)); 2214 PetscFunctionReturn(0); 2215 } 2216 2217 static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3) { 2218 Vec temp_vec; 2219 2220 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 2221 if (v3 != v2) { 2222 PetscCall(MatMultTranspose(A, v1, v3)); 2223 PetscCall(VecAXPY(v3, 1.0, v2)); 2224 } else { 2225 PetscCall(VecDuplicate(v2, &temp_vec)); 2226 PetscCall(MatMultTranspose(A, v1, temp_vec)); 2227 PetscCall(VecAXPY(temp_vec, 1.0, v2)); 2228 PetscCall(VecCopy(temp_vec, v3)); 2229 PetscCall(VecDestroy(&temp_vec)); 2230 } 2231 PetscFunctionReturn(0); 2232 } 2233 2234 static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer) { 2235 Mat_IS *a = (Mat_IS *)A->data; 2236 PetscViewer sviewer; 2237 PetscBool isascii, view = PETSC_TRUE; 2238 2239 PetscFunctionBegin; 2240 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 2241 if (isascii) { 2242 PetscViewerFormat format; 2243 2244 PetscCall(PetscViewerGetFormat(viewer, &format)); 2245 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 2246 } 2247 if (!view) PetscFunctionReturn(0); 2248 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 2249 PetscCall(MatView(a->A, sviewer)); 2250 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 2251 PetscCall(PetscViewerFlush(viewer)); 2252 PetscFunctionReturn(0); 2253 } 2254 2255 static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values) { 2256 Mat_IS *is = (Mat_IS *)mat->data; 2257 MPI_Datatype nodeType; 2258 const PetscScalar *lv; 2259 PetscInt bs; 2260 2261 PetscFunctionBegin; 2262 PetscCall(MatGetBlockSize(mat, &bs)); 2263 PetscCall(MatSetBlockSize(is->A, bs)); 2264 PetscCall(MatInvertBlockDiagonal(is->A, &lv)); 2265 if (!is->bdiag) PetscCall(PetscMalloc1(bs * mat->rmap->n, &is->bdiag)); 2266 PetscCallMPI(MPI_Type_contiguous(bs, MPIU_SCALAR, &nodeType)); 2267 PetscCallMPI(MPI_Type_commit(&nodeType)); 2268 PetscCall(PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE)); 2269 PetscCall(PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE)); 2270 PetscCallMPI(MPI_Type_free(&nodeType)); 2271 if (values) *values = is->bdiag; 2272 PetscFunctionReturn(0); 2273 } 2274 2275 static PetscErrorCode MatISSetUpScatters_Private(Mat A) { 2276 Vec cglobal, rglobal; 2277 IS from; 2278 Mat_IS *is = (Mat_IS *)A->data; 2279 PetscScalar sum; 2280 const PetscInt *garray; 2281 PetscInt nr, rbs, nc, cbs; 2282 VecType rtype; 2283 2284 PetscFunctionBegin; 2285 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr)); 2286 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs)); 2287 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc)); 2288 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs)); 2289 PetscCall(VecDestroy(&is->x)); 2290 PetscCall(VecDestroy(&is->y)); 2291 PetscCall(VecDestroy(&is->counter)); 2292 PetscCall(VecScatterDestroy(&is->rctx)); 2293 PetscCall(VecScatterDestroy(&is->cctx)); 2294 PetscCall(MatCreateVecs(is->A, &is->x, &is->y)); 2295 PetscCall(VecBindToCPU(is->y, PETSC_TRUE)); 2296 PetscCall(VecGetRootType_Private(is->y, &rtype)); 2297 PetscCall(PetscFree(A->defaultvectype)); 2298 PetscCall(PetscStrallocpy(rtype, &A->defaultvectype)); 2299 PetscCall(MatCreateVecs(A, &cglobal, &rglobal)); 2300 PetscCall(VecBindToCPU(rglobal, PETSC_TRUE)); 2301 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray)); 2302 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from)); 2303 PetscCall(VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx)); 2304 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray)); 2305 PetscCall(ISDestroy(&from)); 2306 if (is->rmapping != is->cmapping) { 2307 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray)); 2308 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from)); 2309 PetscCall(VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx)); 2310 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray)); 2311 PetscCall(ISDestroy(&from)); 2312 } else { 2313 PetscCall(PetscObjectReference((PetscObject)is->rctx)); 2314 is->cctx = is->rctx; 2315 } 2316 PetscCall(VecDestroy(&cglobal)); 2317 2318 /* interface counter vector (local) */ 2319 PetscCall(VecDuplicate(is->y, &is->counter)); 2320 PetscCall(VecBindToCPU(is->counter, PETSC_TRUE)); 2321 PetscCall(VecSet(is->y, 1.)); 2322 PetscCall(VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE)); 2323 PetscCall(VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE)); 2324 PetscCall(VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD)); 2325 PetscCall(VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD)); 2326 PetscCall(VecBindToCPU(is->y, PETSC_FALSE)); 2327 PetscCall(VecBindToCPU(is->counter, PETSC_FALSE)); 2328 2329 /* special functions for block-diagonal matrices */ 2330 PetscCall(VecSum(rglobal, &sum)); 2331 A->ops->invertblockdiagonal = NULL; 2332 if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS; 2333 PetscCall(VecDestroy(&rglobal)); 2334 2335 /* setup SF for general purpose shared indices based communications */ 2336 PetscCall(MatISSetUpSF_IS(A)); 2337 PetscFunctionReturn(0); 2338 } 2339 2340 static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap) { 2341 IS is; 2342 ISLocalToGlobalMappingType l2gtype; 2343 const PetscInt *idxs; 2344 PetscHSetI ht; 2345 PetscInt *nidxs; 2346 PetscInt i, n, bs, c; 2347 PetscBool flg[] = {PETSC_FALSE, PETSC_FALSE}; 2348 2349 PetscFunctionBegin; 2350 PetscCall(ISLocalToGlobalMappingGetSize(map, &n)); 2351 PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs)); 2352 PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs)); 2353 PetscCall(PetscHSetICreate(&ht)); 2354 PetscCall(PetscMalloc1(n / bs, &nidxs)); 2355 for (i = 0, c = 0; i < n / bs; i++) { 2356 PetscBool missing; 2357 if (idxs[i] < 0) { 2358 flg[0] = PETSC_TRUE; 2359 continue; 2360 } 2361 PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing)); 2362 if (!missing) flg[1] = PETSC_TRUE; 2363 else nidxs[c++] = idxs[i]; 2364 } 2365 PetscCall(PetscHSetIDestroy(&ht)); 2366 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A))); 2367 if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */ 2368 *nmap = NULL; 2369 *lmap = NULL; 2370 PetscCall(PetscFree(nidxs)); 2371 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs)); 2372 PetscFunctionReturn(0); 2373 } 2374 2375 /* New l2g map without negative or repeated indices */ 2376 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is)); 2377 PetscCall(ISLocalToGlobalMappingCreateIS(is, nmap)); 2378 PetscCall(ISDestroy(&is)); 2379 PetscCall(ISLocalToGlobalMappingGetType(map, &l2gtype)); 2380 PetscCall(ISLocalToGlobalMappingSetType(*nmap, l2gtype)); 2381 2382 /* New local l2g map for repeated indices */ 2383 PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs)); 2384 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is)); 2385 PetscCall(ISLocalToGlobalMappingCreateIS(is, lmap)); 2386 PetscCall(ISDestroy(&is)); 2387 2388 PetscCall(PetscFree(nidxs)); 2389 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs)); 2390 PetscFunctionReturn(0); 2391 } 2392 2393 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping) { 2394 Mat_IS *is = (Mat_IS *)A->data; 2395 ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL; 2396 PetscBool cong, freem[] = {PETSC_FALSE, PETSC_FALSE}; 2397 PetscInt nr, rbs, nc, cbs; 2398 2399 PetscFunctionBegin; 2400 if (rmapping) PetscCheckSameComm(A, 1, rmapping, 2); 2401 if (cmapping) PetscCheckSameComm(A, 1, cmapping, 3); 2402 2403 PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping)); 2404 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2405 PetscCall(PetscLayoutSetUp(A->rmap)); 2406 PetscCall(PetscLayoutSetUp(A->cmap)); 2407 PetscCall(MatHasCongruentLayouts(A, &cong)); 2408 2409 /* If NULL, local space matches global space */ 2410 if (!rmapping) { 2411 IS is; 2412 2413 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is)); 2414 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmapping)); 2415 if (A->rmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs)); 2416 PetscCall(ISDestroy(&is)); 2417 freem[0] = PETSC_TRUE; 2418 if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping; 2419 } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */ 2420 PetscCall(MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping)); 2421 if (rmapping == cmapping) { 2422 PetscCall(PetscObjectReference((PetscObject)is->rmapping)); 2423 is->cmapping = is->rmapping; 2424 PetscCall(PetscObjectReference((PetscObject)localrmapping)); 2425 localcmapping = localrmapping; 2426 } 2427 } 2428 if (!cmapping) { 2429 IS is; 2430 2431 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is)); 2432 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmapping)); 2433 if (A->cmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs)); 2434 PetscCall(ISDestroy(&is)); 2435 freem[1] = PETSC_TRUE; 2436 } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */ 2437 PetscCall(MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping)); 2438 } 2439 if (!is->rmapping) { 2440 PetscCall(PetscObjectReference((PetscObject)rmapping)); 2441 is->rmapping = rmapping; 2442 } 2443 if (!is->cmapping) { 2444 PetscCall(PetscObjectReference((PetscObject)cmapping)); 2445 is->cmapping = cmapping; 2446 } 2447 2448 /* Clean up */ 2449 PetscCall(MatDestroy(&is->A)); 2450 if (is->csf != is->sf) { 2451 PetscCall(PetscSFDestroy(&is->csf)); 2452 PetscCall(PetscFree2(is->csf_rootdata, is->csf_leafdata)); 2453 } else is->csf = NULL; 2454 PetscCall(PetscSFDestroy(&is->sf)); 2455 PetscCall(PetscFree2(is->sf_rootdata, is->sf_leafdata)); 2456 PetscCall(PetscFree(is->bdiag)); 2457 2458 /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case 2459 (DOLFIN passes 2 different objects) */ 2460 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr)); 2461 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs)); 2462 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc)); 2463 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs)); 2464 if (is->rmapping != is->cmapping && cong) { 2465 PetscBool same = PETSC_FALSE; 2466 if (nr == nc && cbs == rbs) { 2467 const PetscInt *idxs1, *idxs2; 2468 2469 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1)); 2470 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2)); 2471 PetscCall(PetscArraycmp(idxs1, idxs2, nr / rbs, &same)); 2472 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1)); 2473 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2)); 2474 } 2475 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 2476 if (same) { 2477 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2478 PetscCall(PetscObjectReference((PetscObject)is->rmapping)); 2479 is->cmapping = is->rmapping; 2480 } 2481 } 2482 PetscCall(PetscLayoutSetBlockSize(A->rmap, rbs)); 2483 PetscCall(PetscLayoutSetBlockSize(A->cmap, cbs)); 2484 /* Pass the user defined maps to the layout */ 2485 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping)); 2486 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping)); 2487 if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping)); 2488 if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping)); 2489 2490 /* Create the local matrix A */ 2491 PetscCall(MatCreate(PETSC_COMM_SELF, &is->A)); 2492 PetscCall(MatSetType(is->A, is->lmattype)); 2493 PetscCall(MatSetSizes(is->A, nr, nc, nr, nc)); 2494 PetscCall(MatSetBlockSizes(is->A, rbs, cbs)); 2495 PetscCall(MatSetOptionsPrefix(is->A, "is_")); 2496 PetscCall(MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix)); 2497 PetscCall(PetscLayoutSetUp(is->A->rmap)); 2498 PetscCall(PetscLayoutSetUp(is->A->cmap)); 2499 PetscCall(MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping)); 2500 PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping)); 2501 PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping)); 2502 2503 /* setup scatters and local vectors for MatMult */ 2504 if (!is->islocalref) PetscCall(MatISSetUpScatters_Private(A)); 2505 A->preallocated = PETSC_TRUE; 2506 PetscFunctionReturn(0); 2507 } 2508 2509 static PetscErrorCode MatSetUp_IS(Mat A) { 2510 ISLocalToGlobalMapping rmap, cmap; 2511 2512 PetscFunctionBegin; 2513 PetscCall(MatGetLocalToGlobalMapping(A, &rmap, &cmap)); 2514 if (!rmap && !cmap) PetscCall(MatSetLocalToGlobalMapping(A, NULL, NULL)); 2515 PetscFunctionReturn(0); 2516 } 2517 2518 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) { 2519 Mat_IS *is = (Mat_IS *)mat->data; 2520 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2521 2522 PetscFunctionBegin; 2523 PetscCall(ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l)); 2524 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2525 PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l)); 2526 PetscCall(MatSetValues(is->A, m, rows_l, n, cols_l, values, addv)); 2527 } else { 2528 PetscCall(MatSetValues(is->A, m, rows_l, m, rows_l, values, addv)); 2529 } 2530 PetscFunctionReturn(0); 2531 } 2532 2533 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) { 2534 Mat_IS *is = (Mat_IS *)mat->data; 2535 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2536 2537 PetscFunctionBegin; 2538 PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l)); 2539 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2540 PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l)); 2541 PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv)); 2542 } else { 2543 PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, rows_l, values, addv)); 2544 } 2545 PetscFunctionReturn(0); 2546 } 2547 2548 static PetscErrorCode MatSetValuesLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) { 2549 Mat_IS *is = (Mat_IS *)A->data; 2550 2551 PetscFunctionBegin; 2552 if (is->A->rmap->mapping || is->A->cmap->mapping) { 2553 PetscCall(MatSetValuesLocal(is->A, m, rows, n, cols, values, addv)); 2554 } else { 2555 PetscCall(MatSetValues(is->A, m, rows, n, cols, values, addv)); 2556 } 2557 PetscFunctionReturn(0); 2558 } 2559 2560 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) { 2561 Mat_IS *is = (Mat_IS *)A->data; 2562 2563 PetscFunctionBegin; 2564 if (is->A->rmap->mapping || is->A->cmap->mapping) { 2565 PetscCall(MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv)); 2566 } else { 2567 PetscCall(MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv)); 2568 } 2569 PetscFunctionReturn(0); 2570 } 2571 2572 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, PetscBool columns) { 2573 Mat_IS *is = (Mat_IS *)A->data; 2574 2575 PetscFunctionBegin; 2576 if (!n) { 2577 is->pure_neumann = PETSC_TRUE; 2578 } else { 2579 PetscInt i; 2580 is->pure_neumann = PETSC_FALSE; 2581 2582 if (columns) { 2583 PetscCall(MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL)); 2584 } else { 2585 PetscCall(MatZeroRows(is->A, n, rows, diag, NULL, NULL)); 2586 } 2587 if (diag != 0.) { 2588 const PetscScalar *array; 2589 PetscCall(VecGetArrayRead(is->counter, &array)); 2590 for (i = 0; i < n; i++) PetscCall(MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES)); 2591 PetscCall(VecRestoreArrayRead(is->counter, &array)); 2592 } 2593 PetscCall(MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY)); 2594 PetscCall(MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY)); 2595 } 2596 PetscFunctionReturn(0); 2597 } 2598 2599 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns) { 2600 Mat_IS *matis = (Mat_IS *)A->data; 2601 PetscInt nr, nl, len, i; 2602 PetscInt *lrows; 2603 2604 PetscFunctionBegin; 2605 if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) { 2606 PetscBool cong; 2607 2608 PetscCall(PetscLayoutCompare(A->rmap, A->cmap, &cong)); 2609 cong = (PetscBool)(cong && matis->sf == matis->csf); 2610 PetscCheck(cong || !columns, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS"); 2611 PetscCheck(cong || diag == 0., PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS"); 2612 PetscCheck(cong || !x || !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "A->rmap and A->cmap need to be congruent, and the l2g maps be the same"); 2613 } 2614 /* get locally owned rows */ 2615 PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL)); 2616 /* fix right hand side if needed */ 2617 if (x && b) { 2618 const PetscScalar *xx; 2619 PetscScalar *bb; 2620 2621 PetscCall(VecGetArrayRead(x, &xx)); 2622 PetscCall(VecGetArray(b, &bb)); 2623 for (i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]]; 2624 PetscCall(VecRestoreArrayRead(x, &xx)); 2625 PetscCall(VecRestoreArray(b, &bb)); 2626 } 2627 /* get rows associated to the local matrices */ 2628 PetscCall(MatGetSize(matis->A, &nl, NULL)); 2629 PetscCall(PetscArrayzero(matis->sf_leafdata, nl)); 2630 PetscCall(PetscArrayzero(matis->sf_rootdata, A->rmap->n)); 2631 for (i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1; 2632 PetscCall(PetscFree(lrows)); 2633 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 2634 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 2635 PetscCall(PetscMalloc1(nl, &lrows)); 2636 for (i = 0, nr = 0; i < nl; i++) 2637 if (matis->sf_leafdata[i]) lrows[nr++] = i; 2638 PetscCall(MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns)); 2639 PetscCall(PetscFree(lrows)); 2640 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 2641 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 2642 PetscFunctionReturn(0); 2643 } 2644 2645 static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) { 2646 PetscFunctionBegin; 2647 PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE)); 2648 PetscFunctionReturn(0); 2649 } 2650 2651 static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) { 2652 PetscFunctionBegin; 2653 PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE)); 2654 PetscFunctionReturn(0); 2655 } 2656 2657 static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type) { 2658 Mat_IS *is = (Mat_IS *)A->data; 2659 2660 PetscFunctionBegin; 2661 PetscCall(MatAssemblyBegin(is->A, type)); 2662 PetscFunctionReturn(0); 2663 } 2664 2665 static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type) { 2666 Mat_IS *is = (Mat_IS *)A->data; 2667 PetscBool lnnz; 2668 2669 PetscFunctionBegin; 2670 PetscCall(MatAssemblyEnd(is->A, type)); 2671 /* fix for local empty rows/cols */ 2672 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2673 Mat newlA; 2674 ISLocalToGlobalMapping rl2g, cl2g; 2675 IS nzr, nzc; 2676 PetscInt nr, nc, nnzr, nnzc; 2677 PetscBool lnewl2g, newl2g; 2678 2679 PetscCall(MatGetSize(is->A, &nr, &nc)); 2680 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr)); 2681 if (!nzr) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr)); 2682 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc)); 2683 if (!nzc) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc)); 2684 PetscCall(ISGetSize(nzr, &nnzr)); 2685 PetscCall(ISGetSize(nzc, &nnzc)); 2686 if (nnzr != nr || nnzc != nc) { /* need new global l2g map */ 2687 lnewl2g = PETSC_TRUE; 2688 PetscCallMPI(MPI_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A))); 2689 2690 /* extract valid submatrix */ 2691 PetscCall(MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA)); 2692 } else { /* local matrix fully populated */ 2693 lnewl2g = PETSC_FALSE; 2694 PetscCallMPI(MPI_Allreduce(&lnewl2g, &newl2g, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A))); 2695 PetscCall(PetscObjectReference((PetscObject)is->A)); 2696 newlA = is->A; 2697 } 2698 2699 /* attach new global l2g map if needed */ 2700 if (newl2g) { 2701 IS zr, zc; 2702 const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs; 2703 PetscInt *nidxs, i; 2704 2705 PetscCall(ISComplement(nzr, 0, nr, &zr)); 2706 PetscCall(ISComplement(nzc, 0, nc, &zc)); 2707 PetscCall(PetscMalloc1(PetscMax(nr, nc), &nidxs)); 2708 PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs)); 2709 PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs)); 2710 PetscCall(ISGetIndices(zr, &zridxs)); 2711 PetscCall(ISGetIndices(zc, &zcidxs)); 2712 PetscCall(ISGetLocalSize(zr, &nnzr)); 2713 PetscCall(ISGetLocalSize(zc, &nnzc)); 2714 2715 PetscCall(PetscArraycpy(nidxs, ridxs, nr)); 2716 for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1; 2717 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g)); 2718 PetscCall(PetscArraycpy(nidxs, cidxs, nc)); 2719 for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1; 2720 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g)); 2721 2722 PetscCall(ISRestoreIndices(zr, &zridxs)); 2723 PetscCall(ISRestoreIndices(zc, &zcidxs)); 2724 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs)); 2725 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs)); 2726 PetscCall(ISDestroy(&nzr)); 2727 PetscCall(ISDestroy(&nzc)); 2728 PetscCall(ISDestroy(&zr)); 2729 PetscCall(ISDestroy(&zc)); 2730 PetscCall(PetscFree(nidxs)); 2731 PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g)); 2732 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 2733 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 2734 } 2735 PetscCall(MatISSetLocalMat(A, newlA)); 2736 PetscCall(MatDestroy(&newlA)); 2737 PetscCall(ISDestroy(&nzr)); 2738 PetscCall(ISDestroy(&nzc)); 2739 is->locempty = PETSC_FALSE; 2740 } 2741 lnnz = (PetscBool)(is->A->nonzerostate == is->lnnzstate); 2742 is->lnnzstate = is->A->nonzerostate; 2743 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A))); 2744 if (lnnz) A->nonzerostate++; 2745 PetscFunctionReturn(0); 2746 } 2747 2748 static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local) { 2749 Mat_IS *is = (Mat_IS *)mat->data; 2750 2751 PetscFunctionBegin; 2752 *local = is->A; 2753 PetscFunctionReturn(0); 2754 } 2755 2756 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local) { 2757 PetscFunctionBegin; 2758 *local = NULL; 2759 PetscFunctionReturn(0); 2760 } 2761 2762 /*@ 2763 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 2764 2765 Input Parameter: 2766 . mat - the matrix 2767 2768 Output Parameter: 2769 . local - the local matrix 2770 2771 Level: advanced 2772 2773 Notes: 2774 This can be called if you have precomputed the nonzero structure of the 2775 matrix and want to provide it to the inner matrix object to improve the performance 2776 of the MatSetValues() operation. 2777 2778 Call MatISRestoreLocalMat() when finished with the local matrix. 2779 2780 .seealso: `MATIS` 2781 @*/ 2782 PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local) { 2783 PetscFunctionBegin; 2784 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2785 PetscValidPointer(local, 2); 2786 PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local)); 2787 PetscFunctionReturn(0); 2788 } 2789 2790 /*@ 2791 MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat() 2792 2793 Input Parameter: 2794 . mat - the matrix 2795 2796 Output Parameter: 2797 . local - the local matrix 2798 2799 Level: advanced 2800 2801 .seealso: `MATIS` 2802 @*/ 2803 PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local) { 2804 PetscFunctionBegin; 2805 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2806 PetscValidPointer(local, 2); 2807 PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local)); 2808 PetscFunctionReturn(0); 2809 } 2810 2811 static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype) { 2812 Mat_IS *is = (Mat_IS *)mat->data; 2813 2814 PetscFunctionBegin; 2815 if (is->A) PetscCall(MatSetType(is->A, mtype)); 2816 PetscCall(PetscFree(is->lmattype)); 2817 PetscCall(PetscStrallocpy(mtype, &is->lmattype)); 2818 PetscFunctionReturn(0); 2819 } 2820 2821 /*@ 2822 MatISSetLocalMatType - Specifies the type of local matrix 2823 2824 Input Parameters: 2825 + mat - the matrix 2826 - mtype - the local matrix type 2827 2828 Output Parameter: 2829 2830 Level: advanced 2831 2832 .seealso: `MATIS`, `MatSetType()`, `MatType` 2833 @*/ 2834 PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype) { 2835 PetscFunctionBegin; 2836 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2837 PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype)); 2838 PetscFunctionReturn(0); 2839 } 2840 2841 static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local) { 2842 Mat_IS *is = (Mat_IS *)mat->data; 2843 PetscInt nrows, ncols, orows, ocols; 2844 MatType mtype, otype; 2845 PetscBool sametype = PETSC_TRUE; 2846 2847 PetscFunctionBegin; 2848 if (is->A && !is->islocalref) { 2849 PetscCall(MatGetSize(is->A, &orows, &ocols)); 2850 PetscCall(MatGetSize(local, &nrows, &ncols)); 2851 PetscCheck(orows == nrows && ocols == ncols, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local MATIS matrix should be of size %" PetscInt_FMT "x%" PetscInt_FMT " (you passed a %" PetscInt_FMT "x%" PetscInt_FMT " matrix)", orows, ocols, nrows, ncols); 2852 PetscCall(MatGetType(local, &mtype)); 2853 PetscCall(MatGetType(is->A, &otype)); 2854 PetscCall(PetscStrcmp(mtype, otype, &sametype)); 2855 } 2856 PetscCall(PetscObjectReference((PetscObject)local)); 2857 PetscCall(MatDestroy(&is->A)); 2858 is->A = local; 2859 PetscCall(MatGetType(is->A, &mtype)); 2860 PetscCall(MatISSetLocalMatType(mat, mtype)); 2861 if (!sametype && !is->islocalref) PetscCall(MatISSetUpScatters_Private(mat)); 2862 PetscFunctionReturn(0); 2863 } 2864 2865 /*@ 2866 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 2867 2868 Collective on Mat 2869 2870 Input Parameters: 2871 + mat - the matrix 2872 - local - the local matrix 2873 2874 Output Parameter: 2875 2876 Level: advanced 2877 2878 Notes: 2879 This can be called if you have precomputed the local matrix and 2880 want to provide it to the matrix object MATIS. 2881 2882 .seealso: `MATIS` 2883 @*/ 2884 PetscErrorCode MatISSetLocalMat(Mat mat, Mat local) { 2885 PetscFunctionBegin; 2886 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2887 PetscValidHeaderSpecific(local, MAT_CLASSID, 2); 2888 PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local)); 2889 PetscFunctionReturn(0); 2890 } 2891 2892 static PetscErrorCode MatZeroEntries_IS(Mat A) { 2893 Mat_IS *a = (Mat_IS *)A->data; 2894 2895 PetscFunctionBegin; 2896 PetscCall(MatZeroEntries(a->A)); 2897 PetscFunctionReturn(0); 2898 } 2899 2900 static PetscErrorCode MatScale_IS(Mat A, PetscScalar a) { 2901 Mat_IS *is = (Mat_IS *)A->data; 2902 2903 PetscFunctionBegin; 2904 PetscCall(MatScale(is->A, a)); 2905 PetscFunctionReturn(0); 2906 } 2907 2908 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) { 2909 Mat_IS *is = (Mat_IS *)A->data; 2910 2911 PetscFunctionBegin; 2912 /* get diagonal of the local matrix */ 2913 PetscCall(MatGetDiagonal(is->A, is->y)); 2914 2915 /* scatter diagonal back into global vector */ 2916 PetscCall(VecSet(v, 0)); 2917 PetscCall(VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE)); 2918 PetscCall(VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE)); 2919 PetscFunctionReturn(0); 2920 } 2921 2922 static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg) { 2923 Mat_IS *a = (Mat_IS *)A->data; 2924 2925 PetscFunctionBegin; 2926 PetscCall(MatSetOption(a->A, op, flg)); 2927 PetscFunctionReturn(0); 2928 } 2929 2930 static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str) { 2931 Mat_IS *y = (Mat_IS *)Y->data; 2932 Mat_IS *x; 2933 2934 PetscFunctionBegin; 2935 if (PetscDefined(USE_DEBUG)) { 2936 PetscBool ismatis; 2937 PetscCall(PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis)); 2938 PetscCheck(ismatis, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 2939 } 2940 x = (Mat_IS *)X->data; 2941 PetscCall(MatAXPY(y->A, a, x->A, str)); 2942 PetscFunctionReturn(0); 2943 } 2944 2945 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat) { 2946 Mat lA; 2947 Mat_IS *matis = (Mat_IS *)(A->data); 2948 ISLocalToGlobalMapping rl2g, cl2g; 2949 IS is; 2950 const PetscInt *rg, *rl; 2951 PetscInt nrg; 2952 PetscInt N, M, nrl, i, *idxs; 2953 2954 PetscFunctionBegin; 2955 PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg)); 2956 PetscCall(ISGetLocalSize(row, &nrl)); 2957 PetscCall(ISGetIndices(row, &rl)); 2958 PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg)); 2959 if (PetscDefined(USE_DEBUG)) { 2960 for (i = 0; i < nrl; i++) PetscCheck(rl[i] < nrg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local row index %" PetscInt_FMT " -> %" PetscInt_FMT " greater then maximum possible %" PetscInt_FMT, i, rl[i], nrg); 2961 } 2962 PetscCall(PetscMalloc1(nrg, &idxs)); 2963 /* map from [0,nrl) to row */ 2964 for (i = 0; i < nrl; i++) idxs[i] = rl[i]; 2965 for (i = nrl; i < nrg; i++) idxs[i] = -1; 2966 PetscCall(ISRestoreIndices(row, &rl)); 2967 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg)); 2968 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is)); 2969 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g)); 2970 PetscCall(ISDestroy(&is)); 2971 /* compute new l2g map for columns */ 2972 if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) { 2973 const PetscInt *cg, *cl; 2974 PetscInt ncg; 2975 PetscInt ncl; 2976 2977 PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg)); 2978 PetscCall(ISGetLocalSize(col, &ncl)); 2979 PetscCall(ISGetIndices(col, &cl)); 2980 PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg)); 2981 if (PetscDefined(USE_DEBUG)) { 2982 for (i = 0; i < ncl; i++) PetscCheck(cl[i] < ncg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Local column index %" PetscInt_FMT " -> %" PetscInt_FMT " greater then maximum possible %" PetscInt_FMT, i, cl[i], ncg); 2983 } 2984 PetscCall(PetscMalloc1(ncg, &idxs)); 2985 /* map from [0,ncl) to col */ 2986 for (i = 0; i < ncl; i++) idxs[i] = cl[i]; 2987 for (i = ncl; i < ncg; i++) idxs[i] = -1; 2988 PetscCall(ISRestoreIndices(col, &cl)); 2989 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg)); 2990 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is)); 2991 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g)); 2992 PetscCall(ISDestroy(&is)); 2993 } else { 2994 PetscCall(PetscObjectReference((PetscObject)rl2g)); 2995 cl2g = rl2g; 2996 } 2997 /* create the MATIS submatrix */ 2998 PetscCall(MatGetSize(A, &M, &N)); 2999 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), submat)); 3000 PetscCall(MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N)); 3001 PetscCall(MatSetType(*submat, MATIS)); 3002 matis = (Mat_IS *)((*submat)->data); 3003 matis->islocalref = PETSC_TRUE; 3004 PetscCall(MatSetLocalToGlobalMapping(*submat, rl2g, cl2g)); 3005 PetscCall(MatISGetLocalMat(A, &lA)); 3006 PetscCall(MatISSetLocalMat(*submat, lA)); 3007 PetscCall(MatSetUp(*submat)); 3008 PetscCall(MatAssemblyBegin(*submat, MAT_FINAL_ASSEMBLY)); 3009 PetscCall(MatAssemblyEnd(*submat, MAT_FINAL_ASSEMBLY)); 3010 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 3011 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 3012 3013 /* remove unsupported ops */ 3014 PetscCall(PetscMemzero((*submat)->ops, sizeof(struct _MatOps))); 3015 (*submat)->ops->destroy = MatDestroy_IS; 3016 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 3017 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 3018 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 3019 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 3020 PetscFunctionReturn(0); 3021 } 3022 3023 static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems *PetscOptionsObject) { 3024 Mat_IS *a = (Mat_IS *)A->data; 3025 char type[256]; 3026 PetscBool flg; 3027 3028 PetscFunctionBegin; 3029 PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options"); 3030 PetscCall(PetscOptionsBool("-matis_keepassembled", "Store an assembled version if needed", "MatISKeepAssembled", a->keepassembled, &a->keepassembled, NULL)); 3031 PetscCall(PetscOptionsBool("-matis_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL)); 3032 PetscCall(PetscOptionsBool("-matis_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL)); 3033 PetscCall(PetscOptionsFList("-matis_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg)); 3034 if (flg) PetscCall(MatISSetLocalMatType(A, type)); 3035 if (a->A) PetscCall(MatSetFromOptions(a->A)); 3036 PetscOptionsHeadEnd(); 3037 PetscFunctionReturn(0); 3038 } 3039 3040 /*@ 3041 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 3042 process but not across processes. 3043 3044 Input Parameters: 3045 + comm - MPI communicator that will share the matrix 3046 . bs - block size of the matrix 3047 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 3048 . rmap - local to global map for rows 3049 - cmap - local to global map for cols 3050 3051 Output Parameter: 3052 . A - the resulting matrix 3053 3054 Level: advanced 3055 3056 Notes: 3057 See MATIS for more details. 3058 m and n are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors 3059 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 3060 If rmap (cmap) is NULL, then the local row (column) spaces matches the global space. 3061 3062 .seealso: `MATIS`, `MatSetLocalToGlobalMapping()` 3063 @*/ 3064 PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A) { 3065 PetscFunctionBegin; 3066 PetscCall(MatCreate(comm, A)); 3067 PetscCall(MatSetSizes(*A, m, n, M, N)); 3068 if (bs > 0) PetscCall(MatSetBlockSize(*A, bs)); 3069 PetscCall(MatSetType(*A, MATIS)); 3070 PetscCall(MatSetLocalToGlobalMapping(*A, rmap, cmap)); 3071 PetscFunctionReturn(0); 3072 } 3073 3074 static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has) { 3075 Mat_IS *a = (Mat_IS *)A->data; 3076 MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP}; 3077 3078 PetscFunctionBegin; 3079 *has = PETSC_FALSE; 3080 if (!((void **)A->ops)[op] || !a->A) PetscFunctionReturn(0); 3081 *has = PETSC_TRUE; 3082 for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++) 3083 if (op == tobefiltered[i]) PetscFunctionReturn(0); 3084 PetscCall(MatHasOperation(a->A, op, has)); 3085 PetscFunctionReturn(0); 3086 } 3087 3088 static PetscErrorCode MatSetValuesCOO_IS(Mat A, const PetscScalar v[], InsertMode imode) { 3089 Mat_IS *a = (Mat_IS *)A->data; 3090 3091 PetscFunctionBegin; 3092 PetscCall(MatSetValuesCOO(a->A, v, imode)); 3093 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 3094 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 3095 PetscFunctionReturn(0); 3096 } 3097 3098 static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[]) { 3099 Mat_IS *a = (Mat_IS *)A->data; 3100 3101 PetscFunctionBegin; 3102 PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping"); 3103 if (a->A->rmap->mapping || a->A->cmap->mapping) { 3104 PetscCall(MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j)); 3105 } else { 3106 PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j)); 3107 } 3108 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS)); 3109 A->preallocated = PETSC_TRUE; 3110 PetscFunctionReturn(0); 3111 } 3112 3113 static PetscErrorCode MatSetPreallocationCOO_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[]) { 3114 Mat_IS *a = (Mat_IS *)A->data; 3115 3116 PetscFunctionBegin; 3117 PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping"); 3118 PetscCheck(ncoo <= PETSC_MAX_INT, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support", ncoo); 3119 PetscCall(ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo, coo_i, NULL, coo_i)); 3120 PetscCall(ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo, coo_j, NULL, coo_j)); 3121 PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j)); 3122 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS)); 3123 A->preallocated = PETSC_TRUE; 3124 PetscFunctionReturn(0); 3125 } 3126 3127 static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA) { 3128 Mat_IS *a = (Mat_IS *)A->data; 3129 PetscObjectState Astate, aAstate = PETSC_MIN_INT; 3130 PetscObjectState Annzstate, aAnnzstate = PETSC_MIN_INT; 3131 3132 PetscFunctionBegin; 3133 PetscCall(PetscObjectStateGet((PetscObject)A, &Astate)); 3134 Annzstate = A->nonzerostate; 3135 if (a->assembledA) { 3136 PetscCall(PetscObjectStateGet((PetscObject)a->assembledA, &aAstate)); 3137 aAnnzstate = a->assembledA->nonzerostate; 3138 } 3139 if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA)); 3140 if (Astate != aAstate || !a->assembledA) { 3141 MatType aAtype; 3142 PetscMPIInt size; 3143 PetscInt rbs, cbs, bs; 3144 3145 /* the assembled form is used as temporary storage for parallel operations 3146 like createsubmatrices and the like, do not waste device memory */ 3147 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 3148 PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs)); 3149 PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs)); 3150 bs = rbs == cbs ? rbs : 1; 3151 if (a->assembledA) PetscCall(MatGetType(a->assembledA, &aAtype)); 3152 else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ; 3153 else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ; 3154 3155 PetscCall(MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA)); 3156 PetscCall(PetscObjectStateSet((PetscObject)a->assembledA, Astate)); 3157 a->assembledA->nonzerostate = Annzstate; 3158 } 3159 PetscCall(PetscObjectReference((PetscObject)a->assembledA)); 3160 *tA = a->assembledA; 3161 if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA)); 3162 PetscFunctionReturn(0); 3163 } 3164 3165 static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA) { 3166 PetscFunctionBegin; 3167 PetscCall(MatDestroy(tA)); 3168 *tA = NULL; 3169 PetscFunctionReturn(0); 3170 } 3171 3172 static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA) { 3173 Mat_IS *a = (Mat_IS *)A->data; 3174 PetscObjectState Astate, dAstate = PETSC_MIN_INT; 3175 3176 PetscFunctionBegin; 3177 PetscCall(PetscObjectStateGet((PetscObject)A, &Astate)); 3178 if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA, &dAstate)); 3179 if (Astate != dAstate) { 3180 Mat tA; 3181 MatType ltype; 3182 3183 PetscCall(MatDestroy(&a->dA)); 3184 PetscCall(MatISGetAssembled_Private(A, &tA)); 3185 PetscCall(MatGetDiagonalBlock(tA, &a->dA)); 3186 PetscCall(MatPropagateSymmetryOptions(tA, a->dA)); 3187 PetscCall(MatGetType(a->A, <ype)); 3188 PetscCall(MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA)); 3189 PetscCall(PetscObjectReference((PetscObject)a->dA)); 3190 PetscCall(MatISRestoreAssembled_Private(A, &tA)); 3191 PetscCall(PetscObjectStateSet((PetscObject)a->dA, Astate)); 3192 } 3193 *dA = a->dA; 3194 PetscFunctionReturn(0); 3195 } 3196 3197 static PetscErrorCode MatCreateSubMatrices_IS(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse reuse, Mat *submat[]) { 3198 Mat tA; 3199 3200 PetscFunctionBegin; 3201 PetscCall(MatISGetAssembled_Private(A, &tA)); 3202 PetscCall(MatCreateSubMatrices(tA, n, irow, icol, reuse, submat)); 3203 /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */ 3204 #if 0 3205 { 3206 Mat_IS *a = (Mat_IS*)A->data; 3207 MatType ltype; 3208 VecType vtype; 3209 char *flg; 3210 3211 PetscCall(MatGetType(a->A,<ype)); 3212 PetscCall(MatGetVecType(a->A,&vtype)); 3213 PetscCall(PetscStrstr(vtype,"cuda",&flg)); 3214 if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg)); 3215 if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg)); 3216 if (flg) { 3217 for (PetscInt i = 0; i < n; i++) { 3218 Mat sA = (*submat)[i]; 3219 3220 PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA)); 3221 (*submat)[i] = sA; 3222 } 3223 } 3224 } 3225 #endif 3226 PetscCall(MatISRestoreAssembled_Private(A, &tA)); 3227 PetscFunctionReturn(0); 3228 } 3229 3230 static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov) { 3231 Mat tA; 3232 3233 PetscFunctionBegin; 3234 PetscCall(MatISGetAssembled_Private(A, &tA)); 3235 PetscCall(MatIncreaseOverlap(tA, n, is, ov)); 3236 PetscCall(MatISRestoreAssembled_Private(A, &tA)); 3237 PetscFunctionReturn(0); 3238 } 3239 3240 /*@ 3241 MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the MATIS object 3242 3243 Not Collective 3244 3245 Input Parameter: 3246 . A - the matrix 3247 3248 Output Parameters: 3249 + rmapping - row mapping 3250 - cmapping - column mapping 3251 3252 Notes: The returned map can be different from the one used to construct the MATIS object, since it will not contain negative or repeated indices. 3253 3254 Level: advanced 3255 3256 .seealso: `MatSetLocalToGlobalMapping()` 3257 @*/ 3258 PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping) { 3259 PetscFunctionBegin; 3260 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3261 PetscValidType(A, 1); 3262 if (rmapping) PetscValidPointer(rmapping, 2); 3263 if (cmapping) PetscValidPointer(cmapping, 3); 3264 PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping)); 3265 PetscFunctionReturn(0); 3266 } 3267 3268 static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c) { 3269 Mat_IS *a = (Mat_IS *)A->data; 3270 3271 PetscFunctionBegin; 3272 if (r) *r = a->rmapping; 3273 if (c) *c = a->cmapping; 3274 PetscFunctionReturn(0); 3275 } 3276 3277 /*MC 3278 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 3279 This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector 3280 product is handled "implicitly". 3281 3282 Options Database Keys: 3283 + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 3284 . -matis_fixempty - Fixes local matrices in case of empty local rows/columns. 3285 - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP(). 3286 3287 Notes: 3288 Options prefix for the inner matrix are given by -is_mat_xxx 3289 3290 You must call MatSetLocalToGlobalMapping() before using this matrix type. 3291 3292 You can do matrix preallocation on the local matrix after you obtain it with 3293 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 3294 3295 Level: advanced 3296 3297 .seealso: `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP` 3298 3299 M*/ 3300 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) { 3301 Mat_IS *a; 3302 3303 PetscFunctionBegin; 3304 PetscCall(PetscNewLog(A, &a)); 3305 PetscCall(PetscStrallocpy(MATAIJ, &a->lmattype)); 3306 A->data = (void *)a; 3307 3308 /* matrix ops */ 3309 PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps))); 3310 A->ops->mult = MatMult_IS; 3311 A->ops->multadd = MatMultAdd_IS; 3312 A->ops->multtranspose = MatMultTranspose_IS; 3313 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 3314 A->ops->destroy = MatDestroy_IS; 3315 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 3316 A->ops->setvalues = MatSetValues_IS; 3317 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 3318 A->ops->setvalueslocal = MatSetValuesLocal_IS; 3319 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 3320 A->ops->zerorows = MatZeroRows_IS; 3321 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 3322 A->ops->assemblybegin = MatAssemblyBegin_IS; 3323 A->ops->assemblyend = MatAssemblyEnd_IS; 3324 A->ops->view = MatView_IS; 3325 A->ops->zeroentries = MatZeroEntries_IS; 3326 A->ops->scale = MatScale_IS; 3327 A->ops->getdiagonal = MatGetDiagonal_IS; 3328 A->ops->setoption = MatSetOption_IS; 3329 A->ops->ishermitian = MatIsHermitian_IS; 3330 A->ops->issymmetric = MatIsSymmetric_IS; 3331 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 3332 A->ops->duplicate = MatDuplicate_IS; 3333 A->ops->missingdiagonal = MatMissingDiagonal_IS; 3334 A->ops->copy = MatCopy_IS; 3335 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 3336 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 3337 A->ops->axpy = MatAXPY_IS; 3338 A->ops->diagonalset = MatDiagonalSet_IS; 3339 A->ops->shift = MatShift_IS; 3340 A->ops->transpose = MatTranspose_IS; 3341 A->ops->getinfo = MatGetInfo_IS; 3342 A->ops->diagonalscale = MatDiagonalScale_IS; 3343 A->ops->setfromoptions = MatSetFromOptions_IS; 3344 A->ops->setup = MatSetUp_IS; 3345 A->ops->hasoperation = MatHasOperation_IS; 3346 A->ops->getdiagonalblock = MatGetDiagonalBlock_IS; 3347 A->ops->createsubmatrices = MatCreateSubMatrices_IS; 3348 A->ops->increaseoverlap = MatIncreaseOverlap_IS; 3349 3350 /* special MATIS functions */ 3351 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS)); 3352 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS)); 3353 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS)); 3354 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS)); 3355 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetMPIXAIJ_C", MatConvert_IS_XAIJ)); 3356 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS)); 3357 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS)); 3358 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS)); 3359 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS)); 3360 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ)); 3361 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ)); 3362 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ)); 3363 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ)); 3364 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ)); 3365 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ)); 3366 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ)); 3367 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS)); 3368 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS)); 3369 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATIS)); 3370 PetscFunctionReturn(0); 3371 } 3372