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