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