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, MPI_C_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 = B; 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 = NULL; 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 MatZeroRowsLocal_SubMat_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 1311 { 1312 PetscInt *rows_l; 1313 Mat_IS *is = (Mat_IS *)A->data; 1314 1315 PetscFunctionBegin; 1316 PetscCall(PetscMalloc1(n, &rows_l)); 1317 PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, n, rows, rows_l)); 1318 PetscCall(MatZeroRowsLocal(is->islocalref, n, rows_l, diag, x, b)); 1319 PetscCall(PetscFree(rows_l)); 1320 PetscFunctionReturn(PETSC_SUCCESS); 1321 } 1322 1323 static PetscErrorCode MatZeroRowsColumnsLocal_SubMat_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 1324 { 1325 PetscInt *rows_l; 1326 Mat_IS *is = (Mat_IS *)A->data; 1327 1328 PetscFunctionBegin; 1329 PetscCall(PetscMalloc1(n, &rows_l)); 1330 PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping, n, rows, rows_l)); 1331 PetscCall(MatZeroRowsColumnsLocal(is->islocalref, n, rows_l, diag, x, b)); 1332 PetscCall(PetscFree(rows_l)); 1333 PetscFunctionReturn(PETSC_SUCCESS); 1334 } 1335 1336 static PetscErrorCode MatCreateSubMatrix_IS(Mat mat, IS irow, IS icol, MatReuse scall, Mat *newmat) 1337 { 1338 Mat locmat, newlocmat; 1339 Mat_IS *newmatis; 1340 const PetscInt *idxs; 1341 PetscInt i, m, n; 1342 1343 PetscFunctionBegin; 1344 if (scall == MAT_REUSE_MATRIX) { 1345 PetscBool ismatis; 1346 1347 PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATIS, &ismatis)); 1348 PetscCheck(ismatis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Not of MATIS type"); 1349 newmatis = (Mat_IS *)(*newmat)->data; 1350 PetscCheck(newmatis->getsub_ris, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local row IS"); 1351 PetscCheck(newmatis->getsub_cis, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_ARG_WRONG, "Cannot reuse matrix! Misses local col IS"); 1352 } 1353 /* irow and icol may not have duplicate entries */ 1354 if (PetscDefined(USE_DEBUG)) { 1355 Vec rtest, ltest; 1356 const PetscScalar *array; 1357 1358 PetscCall(MatCreateVecs(mat, <est, &rtest)); 1359 PetscCall(ISGetLocalSize(irow, &n)); 1360 PetscCall(ISGetIndices(irow, &idxs)); 1361 for (i = 0; i < n; i++) PetscCall(VecSetValue(rtest, idxs[i], 1.0, ADD_VALUES)); 1362 PetscCall(VecAssemblyBegin(rtest)); 1363 PetscCall(VecAssemblyEnd(rtest)); 1364 PetscCall(VecGetLocalSize(rtest, &n)); 1365 PetscCall(VecGetOwnershipRange(rtest, &m, NULL)); 1366 PetscCall(VecGetArrayRead(rtest, &array)); 1367 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])); 1368 PetscCall(VecRestoreArrayRead(rtest, &array)); 1369 PetscCall(ISRestoreIndices(irow, &idxs)); 1370 PetscCall(ISGetLocalSize(icol, &n)); 1371 PetscCall(ISGetIndices(icol, &idxs)); 1372 for (i = 0; i < n; i++) PetscCall(VecSetValue(ltest, idxs[i], 1.0, ADD_VALUES)); 1373 PetscCall(VecAssemblyBegin(ltest)); 1374 PetscCall(VecAssemblyEnd(ltest)); 1375 PetscCall(VecGetLocalSize(ltest, &n)); 1376 PetscCall(VecGetOwnershipRange(ltest, &m, NULL)); 1377 PetscCall(VecGetArrayRead(ltest, &array)); 1378 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])); 1379 PetscCall(VecRestoreArrayRead(ltest, &array)); 1380 PetscCall(ISRestoreIndices(icol, &idxs)); 1381 PetscCall(VecDestroy(&rtest)); 1382 PetscCall(VecDestroy(<est)); 1383 } 1384 if (scall == MAT_INITIAL_MATRIX) { 1385 Mat_IS *matis = (Mat_IS *)mat->data; 1386 ISLocalToGlobalMapping rl2g; 1387 IS is; 1388 PetscInt *lidxs, *lgidxs, *newgidxs; 1389 PetscInt ll, newloc, irbs, icbs, arbs, acbs, rbs, cbs; 1390 PetscBool cong; 1391 MPI_Comm comm; 1392 1393 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm)); 1394 PetscCall(MatGetBlockSizes(mat, &arbs, &acbs)); 1395 PetscCall(ISGetBlockSize(irow, &irbs)); 1396 PetscCall(ISGetBlockSize(icol, &icbs)); 1397 rbs = arbs == irbs ? irbs : 1; 1398 cbs = acbs == icbs ? icbs : 1; 1399 PetscCall(ISGetLocalSize(irow, &m)); 1400 PetscCall(ISGetLocalSize(icol, &n)); 1401 PetscCall(MatCreate(comm, newmat)); 1402 PetscCall(MatSetType(*newmat, MATIS)); 1403 PetscCall(MatISSetAllowRepeated(*newmat, matis->allow_repeated)); 1404 PetscCall(MatSetSizes(*newmat, m, n, PETSC_DECIDE, PETSC_DECIDE)); 1405 PetscCall(MatSetBlockSizes(*newmat, rbs, cbs)); 1406 /* communicate irow to their owners in the layout */ 1407 PetscCall(ISGetIndices(irow, &idxs)); 1408 PetscCall(PetscLayoutMapLocal(mat->rmap, m, idxs, &ll, &lidxs, &lgidxs)); 1409 PetscCall(ISRestoreIndices(irow, &idxs)); 1410 PetscCall(PetscArrayzero(matis->sf_rootdata, matis->sf->nroots)); 1411 for (i = 0; i < ll; i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i] + 1; 1412 PetscCall(PetscFree(lidxs)); 1413 PetscCall(PetscFree(lgidxs)); 1414 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 1415 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 1416 for (i = 0, newloc = 0; i < matis->sf->nleaves; i++) 1417 if (matis->sf_leafdata[i]) newloc++; 1418 PetscCall(PetscMalloc1(newloc, &newgidxs)); 1419 PetscCall(PetscMalloc1(newloc, &lidxs)); 1420 for (i = 0, newloc = 0; i < matis->sf->nleaves; i++) 1421 if (matis->sf_leafdata[i]) { 1422 lidxs[newloc] = i; 1423 newgidxs[newloc++] = matis->sf_leafdata[i] - 1; 1424 } 1425 PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is)); 1426 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g)); 1427 PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g, rbs)); 1428 PetscCall(ISDestroy(&is)); 1429 /* local is to extract local submatrix */ 1430 newmatis = (Mat_IS *)(*newmat)->data; 1431 PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_ris)); 1432 PetscCall(MatHasCongruentLayouts(mat, &cong)); 1433 if (cong && irow == icol && matis->csf == matis->sf) { 1434 PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, rl2g)); 1435 PetscCall(PetscObjectReference((PetscObject)newmatis->getsub_ris)); 1436 newmatis->getsub_cis = newmatis->getsub_ris; 1437 } else { 1438 ISLocalToGlobalMapping cl2g; 1439 1440 /* communicate icol to their owners in the layout */ 1441 PetscCall(ISGetIndices(icol, &idxs)); 1442 PetscCall(PetscLayoutMapLocal(mat->cmap, n, idxs, &ll, &lidxs, &lgidxs)); 1443 PetscCall(ISRestoreIndices(icol, &idxs)); 1444 PetscCall(PetscArrayzero(matis->csf_rootdata, matis->csf->nroots)); 1445 for (i = 0; i < ll; i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i] + 1; 1446 PetscCall(PetscFree(lidxs)); 1447 PetscCall(PetscFree(lgidxs)); 1448 PetscCall(PetscSFBcastBegin(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE)); 1449 PetscCall(PetscSFBcastEnd(matis->csf, MPIU_INT, matis->csf_rootdata, matis->csf_leafdata, MPI_REPLACE)); 1450 for (i = 0, newloc = 0; i < matis->csf->nleaves; i++) 1451 if (matis->csf_leafdata[i]) newloc++; 1452 PetscCall(PetscMalloc1(newloc, &newgidxs)); 1453 PetscCall(PetscMalloc1(newloc, &lidxs)); 1454 for (i = 0, newloc = 0; i < matis->csf->nleaves; i++) 1455 if (matis->csf_leafdata[i]) { 1456 lidxs[newloc] = i; 1457 newgidxs[newloc++] = matis->csf_leafdata[i] - 1; 1458 } 1459 PetscCall(ISCreateGeneral(comm, newloc, newgidxs, PETSC_OWN_POINTER, &is)); 1460 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g)); 1461 PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g, cbs)); 1462 PetscCall(ISDestroy(&is)); 1463 /* local is to extract local submatrix */ 1464 PetscCall(ISCreateGeneral(comm, newloc, lidxs, PETSC_OWN_POINTER, &newmatis->getsub_cis)); 1465 PetscCall(MatSetLocalToGlobalMapping(*newmat, rl2g, cl2g)); 1466 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 1467 } 1468 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 1469 } else { 1470 PetscCall(MatISGetLocalMat(*newmat, &newlocmat)); 1471 } 1472 PetscCall(MatISGetLocalMat(mat, &locmat)); 1473 newmatis = (Mat_IS *)(*newmat)->data; 1474 PetscCall(MatCreateSubMatrix(locmat, newmatis->getsub_ris, newmatis->getsub_cis, scall, &newlocmat)); 1475 if (scall == MAT_INITIAL_MATRIX) { 1476 PetscCall(MatISSetLocalMat(*newmat, newlocmat)); 1477 PetscCall(MatDestroy(&newlocmat)); 1478 } 1479 PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY)); 1480 PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY)); 1481 PetscFunctionReturn(PETSC_SUCCESS); 1482 } 1483 1484 static PetscErrorCode MatCopy_IS(Mat A, Mat B, MatStructure str) 1485 { 1486 Mat_IS *a = (Mat_IS *)A->data, *b; 1487 PetscBool ismatis; 1488 1489 PetscFunctionBegin; 1490 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATIS, &ismatis)); 1491 PetscCheck(ismatis, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Need to be implemented"); 1492 b = (Mat_IS *)B->data; 1493 PetscCall(MatCopy(a->A, b->A, str)); 1494 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 1495 PetscFunctionReturn(PETSC_SUCCESS); 1496 } 1497 1498 static PetscErrorCode MatMissingDiagonal_IS(Mat A, PetscBool *missing, PetscInt *d) 1499 { 1500 Vec v; 1501 const PetscScalar *array; 1502 PetscInt i, n; 1503 1504 PetscFunctionBegin; 1505 *missing = PETSC_FALSE; 1506 PetscCall(MatCreateVecs(A, NULL, &v)); 1507 PetscCall(MatGetDiagonal(A, v)); 1508 PetscCall(VecGetLocalSize(v, &n)); 1509 PetscCall(VecGetArrayRead(v, &array)); 1510 for (i = 0; i < n; i++) 1511 if (array[i] == 0.) break; 1512 PetscCall(VecRestoreArrayRead(v, &array)); 1513 PetscCall(VecDestroy(&v)); 1514 if (i != n) *missing = PETSC_TRUE; 1515 if (d) { 1516 *d = -1; 1517 if (*missing) { 1518 PetscInt rstart; 1519 PetscCall(MatGetOwnershipRange(A, &rstart, NULL)); 1520 *d = i + rstart; 1521 } 1522 } 1523 PetscFunctionReturn(PETSC_SUCCESS); 1524 } 1525 1526 static PetscErrorCode MatISSetUpSF_IS(Mat B) 1527 { 1528 Mat_IS *matis = (Mat_IS *)B->data; 1529 const PetscInt *gidxs; 1530 PetscInt nleaves; 1531 1532 PetscFunctionBegin; 1533 if (matis->sf) PetscFunctionReturn(PETSC_SUCCESS); 1534 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->sf)); 1535 PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping, &gidxs)); 1536 PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nleaves)); 1537 PetscCall(PetscSFSetGraphLayout(matis->sf, B->rmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs)); 1538 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping, &gidxs)); 1539 PetscCall(PetscMalloc2(matis->sf->nroots, &matis->sf_rootdata, matis->sf->nleaves, &matis->sf_leafdata)); 1540 if (matis->rmapping != matis->cmapping) { /* setup SF for columns */ 1541 PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nleaves)); 1542 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B), &matis->csf)); 1543 PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping, &gidxs)); 1544 PetscCall(PetscSFSetGraphLayout(matis->csf, B->cmap, nleaves, NULL, PETSC_OWN_POINTER, gidxs)); 1545 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping, &gidxs)); 1546 PetscCall(PetscMalloc2(matis->csf->nroots, &matis->csf_rootdata, matis->csf->nleaves, &matis->csf_leafdata)); 1547 } else { 1548 matis->csf = matis->sf; 1549 matis->csf_leafdata = matis->sf_leafdata; 1550 matis->csf_rootdata = matis->sf_rootdata; 1551 } 1552 PetscFunctionReturn(PETSC_SUCCESS); 1553 } 1554 1555 /*@ 1556 MatISGetAllowRepeated - Get the flag to allow repeated entries in the local to global map 1557 1558 Not Collective 1559 1560 Input Parameter: 1561 . A - the matrix 1562 1563 Output Parameter: 1564 . flg - the boolean flag 1565 1566 Level: intermediate 1567 1568 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISSetAllowRepeated()` 1569 @*/ 1570 PetscErrorCode MatISGetAllowRepeated(Mat A, PetscBool *flg) 1571 { 1572 PetscBool ismatis; 1573 1574 PetscFunctionBegin; 1575 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1576 PetscAssertPointer(flg, 2); 1577 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis)); 1578 PetscCheck(ismatis, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)A)->type_name); 1579 *flg = ((Mat_IS *)A->data)->allow_repeated; 1580 PetscFunctionReturn(PETSC_SUCCESS); 1581 } 1582 1583 /*@ 1584 MatISSetAllowRepeated - Set the flag to allow repeated entries in the local to global map 1585 1586 Logically Collective 1587 1588 Input Parameters: 1589 + A - the matrix 1590 - flg - the boolean flag 1591 1592 Level: intermediate 1593 1594 Notes: 1595 The default value is `PETSC_FALSE`. 1596 When called AFTER calling `MatSetLocalToGlobalMapping()` it will recreate the local matrices 1597 if `flg` is different from the previously set value. 1598 Specifically, when `flg` is true it will just recreate the local matrices, while if 1599 `flg` is false will assemble the local matrices summing up repeated entries. 1600 1601 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatSetLocalToGlobalMapping()`, `MatISGetAllowRepeated()` 1602 @*/ 1603 PetscErrorCode MatISSetAllowRepeated(Mat A, PetscBool flg) 1604 { 1605 PetscFunctionBegin; 1606 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1607 PetscValidType(A, 1); 1608 PetscValidLogicalCollectiveBool(A, flg, 2); 1609 PetscTryMethod(A, "MatISSetAllowRepeated_C", (Mat, PetscBool), (A, flg)); 1610 PetscFunctionReturn(PETSC_SUCCESS); 1611 } 1612 1613 static PetscErrorCode MatISSetAllowRepeated_IS(Mat A, PetscBool flg) 1614 { 1615 Mat_IS *matis = (Mat_IS *)A->data; 1616 Mat lA = NULL; 1617 ISLocalToGlobalMapping lrmap, lcmap; 1618 1619 PetscFunctionBegin; 1620 if (flg == matis->allow_repeated) PetscFunctionReturn(PETSC_SUCCESS); 1621 if (!matis->A) { /* matrix has not been preallocated yet */ 1622 matis->allow_repeated = flg; 1623 PetscFunctionReturn(PETSC_SUCCESS); 1624 } 1625 PetscCheck(!matis->islocalref, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented for local references"); 1626 if (matis->allow_repeated) { /* we will assemble the old local matrix if needed */ 1627 lA = matis->A; 1628 PetscCall(PetscObjectReference((PetscObject)lA)); 1629 } 1630 /* In case flg is True, we only recreate the local matrix */ 1631 matis->allow_repeated = flg; 1632 PetscCall(MatSetLocalToGlobalMapping(A, A->rmap->mapping, A->cmap->mapping)); 1633 if (lA) { /* assemble previous local matrix if needed */ 1634 Mat nA = matis->A; 1635 1636 PetscCall(MatGetLocalToGlobalMapping(nA, &lrmap, &lcmap)); 1637 if (!lrmap && !lcmap) { 1638 PetscCall(MatISSetLocalMat(A, lA)); 1639 } else { 1640 Mat P = NULL, R = NULL; 1641 MatProductType ptype; 1642 1643 if (lrmap == lcmap) { 1644 ptype = MATPRODUCT_PtAP; 1645 PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P)); 1646 PetscCall(MatProductCreate(lA, P, NULL, &nA)); 1647 } else { 1648 if (lcmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lcmap, nA, PETSC_TRUE, PETSC_FALSE, NULL, &P)); 1649 if (lrmap) PetscCall(MatCreateFromISLocalToGlobalMapping(lrmap, nA, PETSC_FALSE, PETSC_TRUE, NULL, &R)); 1650 if (R && P) { 1651 ptype = MATPRODUCT_ABC; 1652 PetscCall(MatProductCreate(R, lA, P, &nA)); 1653 } else if (R) { 1654 ptype = MATPRODUCT_AB; 1655 PetscCall(MatProductCreate(R, lA, NULL, &nA)); 1656 } else { 1657 ptype = MATPRODUCT_AB; 1658 PetscCall(MatProductCreate(lA, P, NULL, &nA)); 1659 } 1660 } 1661 PetscCall(MatProductSetType(nA, ptype)); 1662 PetscCall(MatProductSetFromOptions(nA)); 1663 PetscCall(MatProductSymbolic(nA)); 1664 PetscCall(MatProductNumeric(nA)); 1665 PetscCall(MatProductClear(nA)); 1666 PetscCall(MatConvert(nA, matis->lmattype, MAT_INPLACE_MATRIX, &nA)); 1667 PetscCall(MatISSetLocalMat(A, nA)); 1668 PetscCall(MatDestroy(&nA)); 1669 PetscCall(MatDestroy(&P)); 1670 PetscCall(MatDestroy(&R)); 1671 } 1672 } 1673 PetscCall(MatDestroy(&lA)); 1674 PetscFunctionReturn(PETSC_SUCCESS); 1675 } 1676 1677 /*@ 1678 MatISStoreL2L - Store local-to-local operators during the Galerkin process of computing `MatPtAP()` 1679 1680 Logically Collective 1681 1682 Input Parameters: 1683 + A - the matrix 1684 - store - the boolean flag 1685 1686 Level: advanced 1687 1688 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()` 1689 @*/ 1690 PetscErrorCode MatISStoreL2L(Mat A, PetscBool store) 1691 { 1692 PetscFunctionBegin; 1693 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1694 PetscValidType(A, 1); 1695 PetscValidLogicalCollectiveBool(A, store, 2); 1696 PetscTryMethod(A, "MatISStoreL2L_C", (Mat, PetscBool), (A, store)); 1697 PetscFunctionReturn(PETSC_SUCCESS); 1698 } 1699 1700 static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store) 1701 { 1702 Mat_IS *matis = (Mat_IS *)A->data; 1703 1704 PetscFunctionBegin; 1705 matis->storel2l = store; 1706 if (!store) PetscCall(PetscObjectCompose((PetscObject)A, "_MatIS_PtAP_l2l", NULL)); 1707 PetscFunctionReturn(PETSC_SUCCESS); 1708 } 1709 1710 /*@ 1711 MatISFixLocalEmpty - Compress out zero local rows from the local matrices 1712 1713 Logically Collective 1714 1715 Input Parameters: 1716 + A - the matrix 1717 - fix - the boolean flag 1718 1719 Level: advanced 1720 1721 Note: 1722 When `fix` is `PETSC_TRUE`, new local matrices and l2g maps are generated during the final assembly process. 1723 1724 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY` 1725 @*/ 1726 PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix) 1727 { 1728 PetscFunctionBegin; 1729 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1730 PetscValidType(A, 1); 1731 PetscValidLogicalCollectiveBool(A, fix, 2); 1732 PetscTryMethod(A, "MatISFixLocalEmpty_C", (Mat, PetscBool), (A, fix)); 1733 PetscFunctionReturn(PETSC_SUCCESS); 1734 } 1735 1736 static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix) 1737 { 1738 Mat_IS *matis = (Mat_IS *)A->data; 1739 1740 PetscFunctionBegin; 1741 matis->locempty = fix; 1742 PetscFunctionReturn(PETSC_SUCCESS); 1743 } 1744 1745 /*@ 1746 MatISSetPreallocation - Preallocates memory for a `MATIS` parallel matrix. 1747 1748 Collective 1749 1750 Input Parameters: 1751 + B - the matrix 1752 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1753 (same value is used for all local rows) 1754 . d_nnz - array containing the number of nonzeros in the various rows of the 1755 DIAGONAL portion of the local submatrix (possibly different for each row) 1756 or `NULL`, if `d_nz` is used to specify the nonzero structure. 1757 The size of this array is equal to the number of local rows, i.e `m`. 1758 For matrices that will be factored, you must leave room for (and set) 1759 the diagonal entry even if it is zero. 1760 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1761 submatrix (same value is used for all local rows). 1762 - o_nnz - array containing the number of nonzeros in the various rows of the 1763 OFF-DIAGONAL portion of the local submatrix (possibly different for 1764 each row) or `NULL`, if `o_nz` is used to specify the nonzero 1765 structure. The size of this array is equal to the number 1766 of local rows, i.e `m`. 1767 1768 If the *_nnz parameter is given then the *_nz parameter is ignored 1769 1770 Level: intermediate 1771 1772 Note: 1773 This function has the same interface as the `MATMPIAIJ` preallocation routine in order to simplify the transition 1774 from the asssembled format to the unassembled one. It overestimates the preallocation of `MATIS` local 1775 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects. 1776 1777 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS` 1778 @*/ 1779 PetscErrorCode MatISSetPreallocation(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[]) 1780 { 1781 PetscFunctionBegin; 1782 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 1783 PetscValidType(B, 1); 1784 PetscTryMethod(B, "MatISSetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (B, d_nz, d_nnz, o_nz, o_nnz)); 1785 PetscFunctionReturn(PETSC_SUCCESS); 1786 } 1787 1788 static PetscErrorCode MatISSetPreallocation_IS(Mat B, PetscInt d_nz, const PetscInt d_nnz[], PetscInt o_nz, const PetscInt o_nnz[]) 1789 { 1790 Mat_IS *matis = (Mat_IS *)B->data; 1791 PetscInt bs, i, nlocalcols; 1792 1793 PetscFunctionBegin; 1794 PetscCall(MatSetUp(B)); 1795 if (!d_nnz) 1796 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nz; 1797 else 1798 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] = d_nnz[i]; 1799 1800 if (!o_nnz) 1801 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nz; 1802 else 1803 for (i = 0; i < matis->sf->nroots; i++) matis->sf_rootdata[i] += o_nnz[i]; 1804 1805 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 1806 PetscCall(MatGetSize(matis->A, NULL, &nlocalcols)); 1807 PetscCall(MatGetBlockSize(matis->A, &bs)); 1808 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 1809 1810 for (i = 0; i < matis->sf->nleaves; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols); 1811 PetscCall(MatSeqAIJSetPreallocation(matis->A, 0, matis->sf_leafdata)); 1812 #if defined(PETSC_HAVE_HYPRE) 1813 PetscCall(MatHYPRESetPreallocation(matis->A, 0, matis->sf_leafdata, 0, NULL)); 1814 #endif 1815 1816 for (i = 0; i < matis->sf->nleaves / bs; i++) { 1817 PetscInt b; 1818 1819 matis->sf_leafdata[i] = matis->sf_leafdata[i * bs] / bs; 1820 for (b = 1; b < bs; b++) matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i], matis->sf_leafdata[i * bs + b] / bs); 1821 } 1822 PetscCall(MatSeqBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata)); 1823 1824 nlocalcols /= bs; 1825 for (i = 0; i < matis->sf->nleaves / bs; i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i], nlocalcols - i); 1826 PetscCall(MatSeqSBAIJSetPreallocation(matis->A, bs, 0, matis->sf_leafdata)); 1827 1828 /* for other matrix types */ 1829 PetscCall(MatSetUp(matis->A)); 1830 PetscFunctionReturn(PETSC_SUCCESS); 1831 } 1832 1833 PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M) 1834 { 1835 Mat_IS *matis = (Mat_IS *)mat->data; 1836 Mat local_mat = NULL, MT; 1837 PetscInt rbs, cbs, rows, cols, lrows, lcols; 1838 PetscInt local_rows, local_cols; 1839 PetscBool isseqdense, isseqsbaij, isseqaij, isseqbaij; 1840 PetscMPIInt size; 1841 const PetscScalar *array; 1842 1843 PetscFunctionBegin; 1844 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat), &size)); 1845 if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N && !matis->allow_repeated) { 1846 Mat B; 1847 IS irows = NULL, icols = NULL; 1848 PetscInt rbs, cbs; 1849 1850 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs)); 1851 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs)); 1852 if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */ 1853 IS rows, cols; 1854 const PetscInt *ridxs, *cidxs; 1855 PetscInt i, nw; 1856 PetscBT work; 1857 1858 PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping, &ridxs)); 1859 PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping, &nw)); 1860 nw = nw / rbs; 1861 PetscCall(PetscBTCreate(nw, &work)); 1862 for (i = 0; i < nw; i++) PetscCall(PetscBTSet(work, ridxs[i])); 1863 for (i = 0; i < nw; i++) 1864 if (!PetscBTLookup(work, i)) break; 1865 if (i == nw) { 1866 PetscCall(ISCreateBlock(PETSC_COMM_SELF, rbs, nw, ridxs, PETSC_USE_POINTER, &rows)); 1867 PetscCall(ISSetPermutation(rows)); 1868 PetscCall(ISInvertPermutation(rows, PETSC_DECIDE, &irows)); 1869 PetscCall(ISDestroy(&rows)); 1870 } 1871 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping, &ridxs)); 1872 PetscCall(PetscBTDestroy(&work)); 1873 if (irows && matis->rmapping != matis->cmapping) { 1874 PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping, &cidxs)); 1875 PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping, &nw)); 1876 nw = nw / cbs; 1877 PetscCall(PetscBTCreate(nw, &work)); 1878 for (i = 0; i < nw; i++) PetscCall(PetscBTSet(work, cidxs[i])); 1879 for (i = 0; i < nw; i++) 1880 if (!PetscBTLookup(work, i)) break; 1881 if (i == nw) { 1882 PetscCall(ISCreateBlock(PETSC_COMM_SELF, cbs, nw, cidxs, PETSC_USE_POINTER, &cols)); 1883 PetscCall(ISSetPermutation(cols)); 1884 PetscCall(ISInvertPermutation(cols, PETSC_DECIDE, &icols)); 1885 PetscCall(ISDestroy(&cols)); 1886 } 1887 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping, &cidxs)); 1888 PetscCall(PetscBTDestroy(&work)); 1889 } else if (irows) { 1890 PetscCall(PetscObjectReference((PetscObject)irows)); 1891 icols = irows; 1892 } 1893 } else { 1894 PetscCall(PetscObjectQuery((PetscObject)*M, "_MatIS_IS_XAIJ_irows", (PetscObject *)&irows)); 1895 PetscCall(PetscObjectQuery((PetscObject)*M, "_MatIS_IS_XAIJ_icols", (PetscObject *)&icols)); 1896 if (irows) PetscCall(PetscObjectReference((PetscObject)irows)); 1897 if (icols) PetscCall(PetscObjectReference((PetscObject)icols)); 1898 } 1899 if (!irows || !icols) { 1900 PetscCall(ISDestroy(&icols)); 1901 PetscCall(ISDestroy(&irows)); 1902 goto general_assembly; 1903 } 1904 PetscCall(MatConvert(matis->A, mtype, MAT_INITIAL_MATRIX, &B)); 1905 if (reuse != MAT_INPLACE_MATRIX) { 1906 PetscCall(MatCreateSubMatrix(B, irows, icols, reuse, M)); 1907 PetscCall(PetscObjectCompose((PetscObject)*M, "_MatIS_IS_XAIJ_irows", (PetscObject)irows)); 1908 PetscCall(PetscObjectCompose((PetscObject)*M, "_MatIS_IS_XAIJ_icols", (PetscObject)icols)); 1909 } else { 1910 Mat C; 1911 1912 PetscCall(MatCreateSubMatrix(B, irows, icols, MAT_INITIAL_MATRIX, &C)); 1913 PetscCall(MatHeaderReplace(mat, &C)); 1914 } 1915 PetscCall(MatDestroy(&B)); 1916 PetscCall(ISDestroy(&icols)); 1917 PetscCall(ISDestroy(&irows)); 1918 PetscFunctionReturn(PETSC_SUCCESS); 1919 } 1920 general_assembly: 1921 PetscCall(MatGetSize(mat, &rows, &cols)); 1922 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping, &rbs)); 1923 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping, &cbs)); 1924 PetscCall(MatGetLocalSize(mat, &lrows, &lcols)); 1925 PetscCall(MatGetSize(matis->A, &local_rows, &local_cols)); 1926 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQDENSE, &isseqdense)); 1927 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQAIJ, &isseqaij)); 1928 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQBAIJ, &isseqbaij)); 1929 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &isseqsbaij)); 1930 PetscCheck(isseqdense || isseqaij || isseqbaij || isseqsbaij, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)matis->A)->type_name); 1931 if (PetscDefined(USE_DEBUG)) { 1932 PetscBool lb[4], bb[4]; 1933 1934 lb[0] = isseqdense; 1935 lb[1] = isseqaij; 1936 lb[2] = isseqbaij; 1937 lb[3] = isseqsbaij; 1938 PetscCallMPI(MPIU_Allreduce(lb, bb, 4, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)mat))); 1939 PetscCheck(bb[0] || bb[1] || bb[2] || bb[3], PETSC_COMM_SELF, PETSC_ERR_SUP, "Local matrices must have the same type"); 1940 } 1941 1942 if (reuse != MAT_REUSE_MATRIX) { 1943 PetscCount ncoo; 1944 PetscInt *coo_i, *coo_j; 1945 1946 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &MT)); 1947 PetscCall(MatSetSizes(MT, lrows, lcols, rows, cols)); 1948 PetscCall(MatSetType(MT, mtype)); 1949 PetscCall(MatSetBlockSizes(MT, rbs, cbs)); 1950 if (!isseqaij && !isseqdense) { 1951 PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat)); 1952 } else { 1953 PetscCall(PetscObjectReference((PetscObject)matis->A)); 1954 local_mat = matis->A; 1955 } 1956 PetscCall(MatSetLocalToGlobalMapping(MT, matis->rmapping, matis->cmapping)); 1957 if (isseqdense) { 1958 PetscInt nr, nc; 1959 1960 PetscCall(MatGetSize(local_mat, &nr, &nc)); 1961 ncoo = nr * nc; 1962 PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j)); 1963 for (PetscInt j = 0; j < nc; j++) { 1964 for (PetscInt i = 0; i < nr; i++) { 1965 coo_i[j * nr + i] = i; 1966 coo_j[j * nr + i] = j; 1967 } 1968 } 1969 } else { 1970 const PetscInt *ii, *jj; 1971 PetscInt nr; 1972 PetscBool done; 1973 1974 PetscCall(MatGetRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nr, &ii, &jj, &done)); 1975 PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatGetRowIJ"); 1976 ncoo = ii[nr]; 1977 PetscCall(PetscMalloc2(ncoo, &coo_i, ncoo, &coo_j)); 1978 PetscCall(PetscArraycpy(coo_j, jj, ncoo)); 1979 for (PetscInt i = 0; i < nr; i++) { 1980 for (PetscInt j = ii[i]; j < ii[i + 1]; j++) coo_i[j] = i; 1981 } 1982 PetscCall(MatRestoreRowIJ(local_mat, 0, PETSC_FALSE, PETSC_FALSE, &nr, &ii, &jj, &done)); 1983 PetscCheck(done, PetscObjectComm((PetscObject)local_mat), PETSC_ERR_PLIB, "Error in MatRestoreRowIJ"); 1984 } 1985 PetscCall(MatSetPreallocationCOOLocal(MT, ncoo, coo_i, coo_j)); 1986 PetscCall(PetscFree2(coo_i, coo_j)); 1987 } else { 1988 PetscInt mrbs, mcbs, mrows, mcols, mlrows, mlcols; 1989 1990 /* some checks */ 1991 MT = *M; 1992 PetscCall(MatGetBlockSizes(MT, &mrbs, &mcbs)); 1993 PetscCall(MatGetSize(MT, &mrows, &mcols)); 1994 PetscCall(MatGetLocalSize(MT, &mlrows, &mlcols)); 1995 PetscCheck(mrows == rows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", rows, mrows); 1996 PetscCheck(mcols == cols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", cols, mcols); 1997 PetscCheck(mlrows == lrows, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")", lrows, mlrows); 1998 PetscCheck(mlcols == lcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")", lcols, mlcols); 1999 PetscCheck(mrbs == rbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", rbs, mrbs); 2000 PetscCheck(mcbs == cbs, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")", cbs, mcbs); 2001 PetscCall(MatZeroEntries(MT)); 2002 if (!isseqaij && !isseqdense) { 2003 PetscCall(MatConvert(matis->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &local_mat)); 2004 } else { 2005 PetscCall(PetscObjectReference((PetscObject)matis->A)); 2006 local_mat = matis->A; 2007 } 2008 } 2009 2010 /* Set values */ 2011 if (isseqdense) { 2012 PetscCall(MatDenseGetArrayRead(local_mat, &array)); 2013 PetscCall(MatSetValuesCOO(MT, array, INSERT_VALUES)); 2014 PetscCall(MatDenseRestoreArrayRead(local_mat, &array)); 2015 } else { 2016 PetscCall(MatSeqAIJGetArrayRead(local_mat, &array)); 2017 PetscCall(MatSetValuesCOO(MT, array, INSERT_VALUES)); 2018 PetscCall(MatSeqAIJRestoreArrayRead(local_mat, &array)); 2019 } 2020 PetscCall(MatDestroy(&local_mat)); 2021 PetscCall(MatAssemblyBegin(MT, MAT_FINAL_ASSEMBLY)); 2022 PetscCall(MatAssemblyEnd(MT, MAT_FINAL_ASSEMBLY)); 2023 if (reuse == MAT_INPLACE_MATRIX) { 2024 PetscCall(MatHeaderReplace(mat, &MT)); 2025 } else if (reuse == MAT_INITIAL_MATRIX) { 2026 *M = MT; 2027 } 2028 PetscFunctionReturn(PETSC_SUCCESS); 2029 } 2030 2031 static PetscErrorCode MatDuplicate_IS(Mat mat, MatDuplicateOption op, Mat *newmat) 2032 { 2033 Mat_IS *matis = (Mat_IS *)mat->data; 2034 PetscInt rbs, cbs, m, n, M, N; 2035 Mat B, localmat; 2036 2037 PetscFunctionBegin; 2038 PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping, &rbs)); 2039 PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping, &cbs)); 2040 PetscCall(MatGetSize(mat, &M, &N)); 2041 PetscCall(MatGetLocalSize(mat, &m, &n)); 2042 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &B)); 2043 PetscCall(MatSetSizes(B, m, n, M, N)); 2044 PetscCall(MatSetBlockSize(B, rbs == cbs ? rbs : 1)); 2045 PetscCall(MatSetType(B, MATIS)); 2046 PetscCall(MatISSetLocalMatType(B, matis->lmattype)); 2047 PetscCall(MatISSetAllowRepeated(B, matis->allow_repeated)); 2048 PetscCall(MatSetLocalToGlobalMapping(B, mat->rmap->mapping, mat->cmap->mapping)); 2049 PetscCall(MatDuplicate(matis->A, op, &localmat)); 2050 PetscCall(MatSetLocalToGlobalMapping(localmat, matis->A->rmap->mapping, matis->A->cmap->mapping)); 2051 PetscCall(MatISSetLocalMat(B, localmat)); 2052 PetscCall(MatDestroy(&localmat)); 2053 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 2054 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 2055 *newmat = B; 2056 PetscFunctionReturn(PETSC_SUCCESS); 2057 } 2058 2059 static PetscErrorCode MatIsHermitian_IS(Mat A, PetscReal tol, PetscBool *flg) 2060 { 2061 Mat_IS *matis = (Mat_IS *)A->data; 2062 PetscBool local_sym; 2063 2064 PetscFunctionBegin; 2065 PetscCall(MatIsHermitian(matis->A, tol, &local_sym)); 2066 PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 2067 PetscFunctionReturn(PETSC_SUCCESS); 2068 } 2069 2070 static PetscErrorCode MatIsSymmetric_IS(Mat A, PetscReal tol, PetscBool *flg) 2071 { 2072 Mat_IS *matis = (Mat_IS *)A->data; 2073 PetscBool local_sym; 2074 2075 PetscFunctionBegin; 2076 if (matis->rmapping != matis->cmapping) { 2077 *flg = PETSC_FALSE; 2078 PetscFunctionReturn(PETSC_SUCCESS); 2079 } 2080 PetscCall(MatIsSymmetric(matis->A, tol, &local_sym)); 2081 PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 2082 PetscFunctionReturn(PETSC_SUCCESS); 2083 } 2084 2085 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A, PetscBool *flg) 2086 { 2087 Mat_IS *matis = (Mat_IS *)A->data; 2088 PetscBool local_sym; 2089 2090 PetscFunctionBegin; 2091 if (matis->rmapping != matis->cmapping) { 2092 *flg = PETSC_FALSE; 2093 PetscFunctionReturn(PETSC_SUCCESS); 2094 } 2095 PetscCall(MatIsStructurallySymmetric(matis->A, &local_sym)); 2096 PetscCallMPI(MPIU_Allreduce(&local_sym, flg, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 2097 PetscFunctionReturn(PETSC_SUCCESS); 2098 } 2099 2100 static PetscErrorCode MatDestroy_IS(Mat A) 2101 { 2102 Mat_IS *b = (Mat_IS *)A->data; 2103 2104 PetscFunctionBegin; 2105 PetscCall(PetscFree(b->bdiag)); 2106 PetscCall(PetscFree(b->lmattype)); 2107 PetscCall(MatDestroy(&b->A)); 2108 PetscCall(VecScatterDestroy(&b->cctx)); 2109 PetscCall(VecScatterDestroy(&b->rctx)); 2110 PetscCall(VecDestroy(&b->x)); 2111 PetscCall(VecDestroy(&b->y)); 2112 PetscCall(VecDestroy(&b->counter)); 2113 PetscCall(ISDestroy(&b->getsub_ris)); 2114 PetscCall(ISDestroy(&b->getsub_cis)); 2115 if (b->sf != b->csf) { 2116 PetscCall(PetscSFDestroy(&b->csf)); 2117 PetscCall(PetscFree2(b->csf_rootdata, b->csf_leafdata)); 2118 } else b->csf = NULL; 2119 PetscCall(PetscSFDestroy(&b->sf)); 2120 PetscCall(PetscFree2(b->sf_rootdata, b->sf_leafdata)); 2121 PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping)); 2122 PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping)); 2123 PetscCall(MatDestroy(&b->dA)); 2124 PetscCall(MatDestroy(&b->assembledA)); 2125 PetscCall(PetscFree(A->data)); 2126 PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL)); 2127 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", NULL)); 2128 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", NULL)); 2129 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", NULL)); 2130 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", NULL)); 2131 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", NULL)); 2132 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", NULL)); 2133 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", NULL)); 2134 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", NULL)); 2135 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", NULL)); 2136 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", NULL)); 2137 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", NULL)); 2138 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", NULL)); 2139 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", NULL)); 2140 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", NULL)); 2141 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", NULL)); 2142 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", NULL)); 2143 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL)); 2144 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL)); 2145 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", NULL)); 2146 PetscFunctionReturn(PETSC_SUCCESS); 2147 } 2148 2149 static PetscErrorCode MatMult_IS(Mat A, Vec x, Vec y) 2150 { 2151 Mat_IS *is = (Mat_IS *)A->data; 2152 PetscScalar zero = 0.0; 2153 2154 PetscFunctionBegin; 2155 /* scatter the global vector x into the local work vector */ 2156 PetscCall(VecScatterBegin(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD)); 2157 PetscCall(VecScatterEnd(is->cctx, x, is->x, INSERT_VALUES, SCATTER_FORWARD)); 2158 2159 /* multiply the local matrix */ 2160 PetscCall(MatMult(is->A, is->x, is->y)); 2161 2162 /* scatter product back into global memory */ 2163 PetscCall(VecSet(y, zero)); 2164 PetscCall(VecScatterBegin(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE)); 2165 PetscCall(VecScatterEnd(is->rctx, is->y, y, ADD_VALUES, SCATTER_REVERSE)); 2166 PetscFunctionReturn(PETSC_SUCCESS); 2167 } 2168 2169 static PetscErrorCode MatMultAdd_IS(Mat A, Vec v1, Vec v2, Vec v3) 2170 { 2171 Vec temp_vec; 2172 2173 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 2174 if (v3 != v2) { 2175 PetscCall(MatMult(A, v1, v3)); 2176 PetscCall(VecAXPY(v3, 1.0, v2)); 2177 } else { 2178 PetscCall(VecDuplicate(v2, &temp_vec)); 2179 PetscCall(MatMult(A, v1, temp_vec)); 2180 PetscCall(VecAXPY(temp_vec, 1.0, v2)); 2181 PetscCall(VecCopy(temp_vec, v3)); 2182 PetscCall(VecDestroy(&temp_vec)); 2183 } 2184 PetscFunctionReturn(PETSC_SUCCESS); 2185 } 2186 2187 static PetscErrorCode MatMultTranspose_IS(Mat A, Vec y, Vec x) 2188 { 2189 Mat_IS *is = (Mat_IS *)A->data; 2190 2191 PetscFunctionBegin; 2192 /* scatter the global vector x into the local work vector */ 2193 PetscCall(VecScatterBegin(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD)); 2194 PetscCall(VecScatterEnd(is->rctx, y, is->y, INSERT_VALUES, SCATTER_FORWARD)); 2195 2196 /* multiply the local matrix */ 2197 PetscCall(MatMultTranspose(is->A, is->y, is->x)); 2198 2199 /* scatter product back into global vector */ 2200 PetscCall(VecSet(x, 0)); 2201 PetscCall(VecScatterBegin(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE)); 2202 PetscCall(VecScatterEnd(is->cctx, is->x, x, ADD_VALUES, SCATTER_REVERSE)); 2203 PetscFunctionReturn(PETSC_SUCCESS); 2204 } 2205 2206 static PetscErrorCode MatMultTransposeAdd_IS(Mat A, Vec v1, Vec v2, Vec v3) 2207 { 2208 Vec temp_vec; 2209 2210 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 2211 if (v3 != v2) { 2212 PetscCall(MatMultTranspose(A, v1, v3)); 2213 PetscCall(VecAXPY(v3, 1.0, v2)); 2214 } else { 2215 PetscCall(VecDuplicate(v2, &temp_vec)); 2216 PetscCall(MatMultTranspose(A, v1, temp_vec)); 2217 PetscCall(VecAXPY(temp_vec, 1.0, v2)); 2218 PetscCall(VecCopy(temp_vec, v3)); 2219 PetscCall(VecDestroy(&temp_vec)); 2220 } 2221 PetscFunctionReturn(PETSC_SUCCESS); 2222 } 2223 2224 static PetscErrorCode ISLocalToGlobalMappingView_Multi(ISLocalToGlobalMapping mapping, PetscInt lsize, PetscInt gsize, const PetscInt vblocks[], PetscViewer viewer) 2225 { 2226 PetscInt tr[3], n; 2227 const PetscInt *indices; 2228 2229 PetscFunctionBegin; 2230 tr[0] = IS_LTOGM_FILE_CLASSID; 2231 tr[1] = 1; 2232 tr[2] = gsize; 2233 PetscCall(PetscViewerBinaryWrite(viewer, tr, 3, PETSC_INT)); 2234 PetscCall(PetscViewerBinaryWriteAll(viewer, vblocks, lsize, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT)); 2235 PetscCall(ISLocalToGlobalMappingGetSize(mapping, &n)); 2236 PetscCall(ISLocalToGlobalMappingGetIndices(mapping, &indices)); 2237 PetscCall(PetscViewerBinaryWriteAll(viewer, indices, n, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT)); 2238 PetscCall(ISLocalToGlobalMappingRestoreIndices(mapping, &indices)); 2239 PetscFunctionReturn(PETSC_SUCCESS); 2240 } 2241 2242 static PetscErrorCode MatView_IS(Mat A, PetscViewer viewer) 2243 { 2244 Mat_IS *a = (Mat_IS *)A->data; 2245 PetscViewer sviewer; 2246 PetscBool isascii, isbinary, viewl2g = PETSC_FALSE, native; 2247 PetscViewerFormat format; 2248 ISLocalToGlobalMapping rmap, cmap; 2249 2250 PetscFunctionBegin; 2251 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 2252 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 2253 PetscCall(PetscViewerGetFormat(viewer, &format)); 2254 native = (PetscBool)(format == PETSC_VIEWER_NATIVE); 2255 if (native) { 2256 rmap = A->rmap->mapping; 2257 cmap = A->cmap->mapping; 2258 } else { 2259 rmap = a->rmapping; 2260 cmap = a->cmapping; 2261 } 2262 if (isascii) { 2263 if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS); 2264 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_MATLAB) viewl2g = PETSC_TRUE; 2265 } else if (isbinary) { 2266 PetscInt tr[6], nr, nc, lsize = 0; 2267 char lmattype[64] = {'\0'}; 2268 PetscMPIInt size; 2269 PetscBool skipHeader, vbs = PETSC_FALSE; 2270 IS is; 2271 const PetscInt *vblocks = NULL; 2272 2273 PetscCall(PetscViewerSetUp(viewer)); 2274 PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)A)->prefix, "-mat_is_view_variableblocksizes", &vbs, NULL)); 2275 if (vbs) { 2276 PetscCall(MatGetVariableBlockSizes(a->A, &lsize, &vblocks)); 2277 PetscCall(PetscMPIIntCast(lsize, &size)); 2278 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &size, 1, MPI_INT, MPI_SUM, PetscObjectComm((PetscObject)viewer))); 2279 } else { 2280 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size)); 2281 } 2282 tr[0] = MAT_FILE_CLASSID; 2283 tr[1] = A->rmap->N; 2284 tr[2] = A->cmap->N; 2285 tr[3] = -size; /* AIJ stores nnz here */ 2286 tr[4] = (PetscInt)(rmap == cmap); 2287 tr[5] = a->allow_repeated; 2288 PetscCall(PetscSNPrintf(lmattype, sizeof(lmattype), "%s", a->lmattype)); 2289 2290 PetscCall(PetscViewerBinaryWrite(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), PETSC_INT)); 2291 PetscCall(PetscViewerBinaryWrite(viewer, lmattype, sizeof(lmattype), PETSC_CHAR)); 2292 2293 /* first dump l2g info (we need the header for proper loading on different number of processes) */ 2294 PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader)); 2295 PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_FALSE)); 2296 if (vbs) { 2297 PetscCall(ISLocalToGlobalMappingView_Multi(rmap, lsize, size, vblocks, viewer)); 2298 if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView_Multi(cmap, lsize, size, vblocks, viewer)); 2299 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), lsize, vblocks, PETSC_USE_POINTER, &is)); 2300 PetscCall(ISView(is, viewer)); 2301 PetscCall(ISView(is, viewer)); 2302 PetscCall(ISDestroy(&is)); 2303 } else { 2304 PetscCall(ISLocalToGlobalMappingView(rmap, viewer)); 2305 if (cmap != rmap) PetscCall(ISLocalToGlobalMappingView(cmap, viewer)); 2306 2307 /* then the sizes of the local matrices */ 2308 PetscCall(MatGetSize(a->A, &nr, &nc)); 2309 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nr, PETSC_USE_POINTER, &is)); 2310 PetscCall(ISView(is, viewer)); 2311 PetscCall(ISDestroy(&is)); 2312 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)viewer), 1, &nc, PETSC_USE_POINTER, &is)); 2313 PetscCall(ISView(is, viewer)); 2314 PetscCall(ISDestroy(&is)); 2315 } 2316 PetscCall(PetscViewerBinarySetSkipHeader(viewer, skipHeader)); 2317 } 2318 if (format == PETSC_VIEWER_ASCII_MATLAB) { 2319 char name[64]; 2320 PetscMPIInt size, rank; 2321 2322 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size)); 2323 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)viewer), &rank)); 2324 if (size > 1) PetscCall(PetscSNPrintf(name, sizeof(name), "lmat_%d", rank)); 2325 else PetscCall(PetscSNPrintf(name, sizeof(name), "lmat")); 2326 PetscCall(PetscObjectSetName((PetscObject)a->A, name)); 2327 } 2328 2329 /* Dump the local matrices */ 2330 if (isbinary) { /* ViewerGetSubViewer does not work in parallel */ 2331 PetscBool isaij; 2332 PetscInt nr, nc; 2333 Mat lA, B; 2334 Mat_MPIAIJ *b; 2335 2336 /* We create a temporary MPIAIJ matrix that stores the unassembled operator */ 2337 PetscCall(PetscObjectBaseTypeCompare((PetscObject)a->A, MATAIJ, &isaij)); 2338 if (!isaij) PetscCall(MatConvert(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &lA)); 2339 else { 2340 PetscCall(PetscObjectReference((PetscObject)a->A)); 2341 lA = a->A; 2342 } 2343 PetscCall(MatCreate(PetscObjectComm((PetscObject)viewer), &B)); 2344 PetscCall(MatSetType(B, MATMPIAIJ)); 2345 PetscCall(MatGetSize(lA, &nr, &nc)); 2346 PetscCall(MatSetSizes(B, nr, nc, PETSC_DECIDE, PETSC_DECIDE)); 2347 PetscCall(MatMPIAIJSetPreallocation(B, 0, NULL, 0, NULL)); 2348 2349 b = (Mat_MPIAIJ *)B->data; 2350 PetscCall(MatDestroy(&b->A)); 2351 b->A = lA; 2352 2353 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 2354 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 2355 PetscCall(MatView(B, viewer)); 2356 PetscCall(MatDestroy(&B)); 2357 } else { 2358 PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 2359 PetscCall(MatView(a->A, sviewer)); 2360 PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer)); 2361 } 2362 2363 /* with ASCII, we dump the l2gmaps at the end */ 2364 if (viewl2g) { 2365 if (format == PETSC_VIEWER_ASCII_MATLAB) { 2366 PetscCall(PetscObjectSetName((PetscObject)rmap, "row")); 2367 PetscCall(ISLocalToGlobalMappingView(rmap, viewer)); 2368 PetscCall(PetscObjectSetName((PetscObject)cmap, "col")); 2369 PetscCall(ISLocalToGlobalMappingView(cmap, viewer)); 2370 } else { 2371 PetscCall(ISLocalToGlobalMappingView(rmap, viewer)); 2372 PetscCall(ISLocalToGlobalMappingView(cmap, viewer)); 2373 } 2374 } 2375 PetscFunctionReturn(PETSC_SUCCESS); 2376 } 2377 2378 static PetscErrorCode ISLocalToGlobalMappingHasRepeatedLocal_Private(ISLocalToGlobalMapping map, PetscBool *has) 2379 { 2380 const PetscInt *idxs; 2381 PetscHSetI ht; 2382 PetscInt n, bs; 2383 2384 PetscFunctionBegin; 2385 PetscCall(ISLocalToGlobalMappingGetSize(map, &n)); 2386 PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs)); 2387 PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs)); 2388 PetscCall(PetscHSetICreate(&ht)); 2389 *has = PETSC_FALSE; 2390 for (PetscInt i = 0; i < n / bs; i++) { 2391 PetscBool missing = PETSC_TRUE; 2392 if (idxs[i] < 0) continue; 2393 PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing)); 2394 if (!missing) { 2395 *has = PETSC_TRUE; 2396 break; 2397 } 2398 } 2399 PetscCall(PetscHSetIDestroy(&ht)); 2400 PetscFunctionReturn(PETSC_SUCCESS); 2401 } 2402 2403 static PetscErrorCode MatLoad_IS(Mat A, PetscViewer viewer) 2404 { 2405 ISLocalToGlobalMapping rmap, cmap; 2406 MPI_Comm comm = PetscObjectComm((PetscObject)A); 2407 PetscBool isbinary, samel, allow, isbaij; 2408 PetscInt tr[6], M, N, nr, nc, Asize, isn; 2409 const PetscInt *idx; 2410 PetscMPIInt size; 2411 char lmattype[64]; 2412 Mat dA, lA; 2413 IS is; 2414 2415 PetscFunctionBegin; 2416 PetscCheckSameComm(A, 1, viewer, 2); 2417 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 2418 PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Invalid viewer of type %s", ((PetscObject)viewer)->type_name); 2419 2420 PetscCall(PetscViewerBinaryRead(viewer, tr, PETSC_STATIC_ARRAY_LENGTH(tr), NULL, PETSC_INT)); 2421 PetscCheck(tr[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix next in file"); 2422 PetscCheck(tr[1] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file"); 2423 PetscCheck(tr[2] >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file"); 2424 PetscCheck(tr[3] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file"); 2425 PetscCheck(tr[4] == 0 || tr[4] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file"); 2426 PetscCheck(tr[5] == 0 || tr[5] == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a IS matrix next in file"); 2427 M = tr[1]; 2428 N = tr[2]; 2429 Asize = -tr[3]; 2430 samel = (PetscBool)tr[4]; 2431 allow = (PetscBool)tr[5]; 2432 PetscCall(PetscViewerBinaryRead(viewer, lmattype, sizeof(lmattype), NULL, PETSC_CHAR)); 2433 2434 /* if we are loading from a larger set of processes, allow repeated entries */ 2435 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)viewer), &size)); 2436 if (Asize > size) allow = PETSC_TRUE; 2437 2438 /* set global sizes if not set already */ 2439 if (A->rmap->N < 0) A->rmap->N = M; 2440 if (A->cmap->N < 0) A->cmap->N = N; 2441 PetscCall(PetscLayoutSetUp(A->rmap)); 2442 PetscCall(PetscLayoutSetUp(A->cmap)); 2443 PetscCheck(M == A->rmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix rows should be %" PetscInt_FMT ", found %" PetscInt_FMT, M, A->rmap->N); 2444 PetscCheck(N == A->cmap->N, comm, PETSC_ERR_ARG_SIZ, "Matrix columns should be %" PetscInt_FMT ", found %" PetscInt_FMT, N, A->cmap->N); 2445 2446 /* load l2g maps */ 2447 PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &rmap)); 2448 PetscCall(ISLocalToGlobalMappingLoad(rmap, viewer)); 2449 if (!samel) { 2450 PetscCall(ISLocalToGlobalMappingCreate(comm, 0, 0, NULL, PETSC_USE_POINTER, &cmap)); 2451 PetscCall(ISLocalToGlobalMappingLoad(cmap, viewer)); 2452 } else { 2453 PetscCall(PetscObjectReference((PetscObject)rmap)); 2454 cmap = rmap; 2455 } 2456 2457 /* load sizes of local matrices */ 2458 PetscCall(ISCreate(comm, &is)); 2459 PetscCall(ISSetType(is, ISGENERAL)); 2460 PetscCall(ISLoad(is, viewer)); 2461 PetscCall(ISGetLocalSize(is, &isn)); 2462 PetscCall(ISGetIndices(is, &idx)); 2463 nr = 0; 2464 for (PetscInt i = 0; i < isn; i++) nr += idx[i]; 2465 PetscCall(ISRestoreIndices(is, &idx)); 2466 PetscCall(ISDestroy(&is)); 2467 PetscCall(ISCreate(comm, &is)); 2468 PetscCall(ISSetType(is, ISGENERAL)); 2469 PetscCall(ISLoad(is, viewer)); 2470 PetscCall(ISGetLocalSize(is, &isn)); 2471 PetscCall(ISGetIndices(is, &idx)); 2472 nc = 0; 2473 for (PetscInt i = 0; i < isn; i++) nc += idx[i]; 2474 PetscCall(ISRestoreIndices(is, &idx)); 2475 PetscCall(ISDestroy(&is)); 2476 2477 /* now load the unassembled operator */ 2478 PetscCall(MatCreate(comm, &dA)); 2479 PetscCall(MatSetType(dA, MATMPIAIJ)); 2480 PetscCall(MatSetSizes(dA, nr, nc, PETSC_DECIDE, PETSC_DECIDE)); 2481 PetscCall(MatLoad(dA, viewer)); 2482 PetscCall(MatMPIAIJGetSeqAIJ(dA, &lA, NULL, NULL)); 2483 PetscCall(PetscObjectReference((PetscObject)lA)); 2484 PetscCall(MatDestroy(&dA)); 2485 2486 /* and convert to the desired format */ 2487 PetscCall(PetscStrcmpAny(lmattype, &isbaij, MATSBAIJ, MATSEQSBAIJ, "")); 2488 if (isbaij) PetscCall(MatSetOption(lA, MAT_SYMMETRIC, PETSC_TRUE)); 2489 PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA)); 2490 2491 /* check if we actually have repeated entries */ 2492 if (allow) { 2493 PetscBool rhas, chas, hasrepeated; 2494 2495 PetscCall(ISLocalToGlobalMappingHasRepeatedLocal_Private(rmap, &rhas)); 2496 if (rmap != cmap) PetscCall(ISLocalToGlobalMappingHasRepeatedLocal_Private(cmap, &chas)); 2497 else chas = rhas; 2498 hasrepeated = (PetscBool)(rhas || chas); 2499 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &hasrepeated, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A))); 2500 if (!hasrepeated) allow = PETSC_FALSE; 2501 } 2502 2503 /* assemble the MATIS object */ 2504 PetscCall(MatISSetAllowRepeated(A, allow)); 2505 PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap)); 2506 PetscCall(MatISSetLocalMat(A, lA)); 2507 PetscCall(MatDestroy(&lA)); 2508 PetscCall(ISLocalToGlobalMappingDestroy(&rmap)); 2509 PetscCall(ISLocalToGlobalMappingDestroy(&cmap)); 2510 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 2511 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 2512 PetscFunctionReturn(PETSC_SUCCESS); 2513 } 2514 2515 static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat, const PetscScalar **values) 2516 { 2517 Mat_IS *is = (Mat_IS *)mat->data; 2518 MPI_Datatype nodeType; 2519 const PetscScalar *lv; 2520 PetscInt bs; 2521 PetscMPIInt mbs; 2522 2523 PetscFunctionBegin; 2524 PetscCall(MatGetBlockSize(mat, &bs)); 2525 PetscCall(MatSetBlockSize(is->A, bs)); 2526 PetscCall(MatInvertBlockDiagonal(is->A, &lv)); 2527 if (!is->bdiag) PetscCall(PetscMalloc1(bs * mat->rmap->n, &is->bdiag)); 2528 PetscCall(PetscMPIIntCast(bs, &mbs)); 2529 PetscCallMPI(MPI_Type_contiguous(mbs, MPIU_SCALAR, &nodeType)); 2530 PetscCallMPI(MPI_Type_commit(&nodeType)); 2531 PetscCall(PetscSFReduceBegin(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE)); 2532 PetscCall(PetscSFReduceEnd(is->sf, nodeType, lv, is->bdiag, MPI_REPLACE)); 2533 PetscCallMPI(MPI_Type_free(&nodeType)); 2534 if (values) *values = is->bdiag; 2535 PetscFunctionReturn(PETSC_SUCCESS); 2536 } 2537 2538 static PetscErrorCode MatISSetUpScatters_Private(Mat A) 2539 { 2540 Vec cglobal, rglobal; 2541 IS from; 2542 Mat_IS *is = (Mat_IS *)A->data; 2543 PetscScalar sum; 2544 const PetscInt *garray; 2545 PetscInt nr, rbs, nc, cbs; 2546 VecType rtype; 2547 2548 PetscFunctionBegin; 2549 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr)); 2550 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs)); 2551 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc)); 2552 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs)); 2553 PetscCall(VecDestroy(&is->x)); 2554 PetscCall(VecDestroy(&is->y)); 2555 PetscCall(VecDestroy(&is->counter)); 2556 PetscCall(VecScatterDestroy(&is->rctx)); 2557 PetscCall(VecScatterDestroy(&is->cctx)); 2558 PetscCall(MatCreateVecs(is->A, &is->x, &is->y)); 2559 PetscCall(VecBindToCPU(is->y, PETSC_TRUE)); 2560 PetscCall(VecGetRootType_Private(is->y, &rtype)); 2561 PetscCall(PetscFree(A->defaultvectype)); 2562 PetscCall(PetscStrallocpy(rtype, &A->defaultvectype)); 2563 PetscCall(MatCreateVecs(A, &cglobal, &rglobal)); 2564 PetscCall(VecBindToCPU(rglobal, PETSC_TRUE)); 2565 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &garray)); 2566 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), rbs, nr / rbs, garray, PETSC_USE_POINTER, &from)); 2567 PetscCall(VecScatterCreate(rglobal, from, is->y, NULL, &is->rctx)); 2568 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &garray)); 2569 PetscCall(ISDestroy(&from)); 2570 if (is->rmapping != is->cmapping) { 2571 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &garray)); 2572 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), cbs, nc / cbs, garray, PETSC_USE_POINTER, &from)); 2573 PetscCall(VecScatterCreate(cglobal, from, is->x, NULL, &is->cctx)); 2574 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &garray)); 2575 PetscCall(ISDestroy(&from)); 2576 } else { 2577 PetscCall(PetscObjectReference((PetscObject)is->rctx)); 2578 is->cctx = is->rctx; 2579 } 2580 PetscCall(VecDestroy(&cglobal)); 2581 2582 /* interface counter vector (local) */ 2583 PetscCall(VecDuplicate(is->y, &is->counter)); 2584 PetscCall(VecBindToCPU(is->counter, PETSC_TRUE)); 2585 PetscCall(VecSet(is->y, 1.)); 2586 PetscCall(VecScatterBegin(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE)); 2587 PetscCall(VecScatterEnd(is->rctx, is->y, rglobal, ADD_VALUES, SCATTER_REVERSE)); 2588 PetscCall(VecScatterBegin(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD)); 2589 PetscCall(VecScatterEnd(is->rctx, rglobal, is->counter, INSERT_VALUES, SCATTER_FORWARD)); 2590 PetscCall(VecBindToCPU(is->y, PETSC_FALSE)); 2591 PetscCall(VecBindToCPU(is->counter, PETSC_FALSE)); 2592 2593 /* special functions for block-diagonal matrices */ 2594 PetscCall(VecSum(rglobal, &sum)); 2595 A->ops->invertblockdiagonal = NULL; 2596 if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS; 2597 PetscCall(VecDestroy(&rglobal)); 2598 2599 /* setup SF for general purpose shared indices based communications */ 2600 PetscCall(MatISSetUpSF_IS(A)); 2601 PetscFunctionReturn(PETSC_SUCCESS); 2602 } 2603 2604 static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap) 2605 { 2606 Mat_IS *matis = (Mat_IS *)A->data; 2607 IS is; 2608 ISLocalToGlobalMappingType l2gtype; 2609 const PetscInt *idxs; 2610 PetscHSetI ht; 2611 PetscInt *nidxs; 2612 PetscInt i, n, bs, c; 2613 PetscBool flg[] = {PETSC_FALSE, PETSC_FALSE}; 2614 2615 PetscFunctionBegin; 2616 PetscCall(ISLocalToGlobalMappingGetSize(map, &n)); 2617 PetscCall(ISLocalToGlobalMappingGetBlockSize(map, &bs)); 2618 PetscCall(ISLocalToGlobalMappingGetBlockIndices(map, &idxs)); 2619 PetscCall(PetscHSetICreate(&ht)); 2620 PetscCall(PetscMalloc1(n / bs, &nidxs)); 2621 for (i = 0, c = 0; i < n / bs; i++) { 2622 PetscBool missing = PETSC_TRUE; 2623 if (idxs[i] < 0) { 2624 flg[0] = PETSC_TRUE; 2625 continue; 2626 } 2627 if (!matis->allow_repeated) PetscCall(PetscHSetIQueryAdd(ht, idxs[i], &missing)); 2628 if (!missing) flg[1] = PETSC_TRUE; 2629 else nidxs[c++] = idxs[i]; 2630 } 2631 PetscCall(PetscHSetIDestroy(&ht)); 2632 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, flg, 2, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A))); 2633 if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */ 2634 *nmap = NULL; 2635 *lmap = NULL; 2636 PetscCall(PetscFree(nidxs)); 2637 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs)); 2638 PetscFunctionReturn(PETSC_SUCCESS); 2639 } 2640 2641 /* New l2g map without negative indices (and repeated indices if not allowed) */ 2642 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A), bs, c, nidxs, PETSC_USE_POINTER, &is)); 2643 PetscCall(ISLocalToGlobalMappingCreateIS(is, nmap)); 2644 PetscCall(ISDestroy(&is)); 2645 PetscCall(ISLocalToGlobalMappingGetType(map, &l2gtype)); 2646 PetscCall(ISLocalToGlobalMappingSetType(*nmap, l2gtype)); 2647 2648 /* New local l2g map for repeated indices if not allowed */ 2649 PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap, IS_GTOLM_MASK, n / bs, idxs, NULL, nidxs)); 2650 PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, n / bs, nidxs, PETSC_USE_POINTER, &is)); 2651 PetscCall(ISLocalToGlobalMappingCreateIS(is, lmap)); 2652 PetscCall(ISDestroy(&is)); 2653 PetscCall(PetscFree(nidxs)); 2654 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map, &idxs)); 2655 PetscFunctionReturn(PETSC_SUCCESS); 2656 } 2657 2658 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping rmapping, ISLocalToGlobalMapping cmapping) 2659 { 2660 Mat_IS *is = (Mat_IS *)A->data; 2661 ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL; 2662 PetscInt nr, rbs, nc, cbs; 2663 PetscBool cong, freem[] = {PETSC_FALSE, PETSC_FALSE}; 2664 2665 PetscFunctionBegin; 2666 if (rmapping) PetscCheckSameComm(A, 1, rmapping, 2); 2667 if (cmapping) PetscCheckSameComm(A, 1, cmapping, 3); 2668 2669 PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping)); 2670 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2671 PetscCall(PetscLayoutSetUp(A->rmap)); 2672 PetscCall(PetscLayoutSetUp(A->cmap)); 2673 PetscCall(MatHasCongruentLayouts(A, &cong)); 2674 2675 /* If NULL, local space matches global space */ 2676 if (!rmapping) { 2677 IS is; 2678 2679 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->rmap->N, 0, 1, &is)); 2680 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmapping)); 2681 PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping, A->rmap->bs)); 2682 PetscCall(ISDestroy(&is)); 2683 freem[0] = PETSC_TRUE; 2684 if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping; 2685 } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */ 2686 PetscCall(MatISFilterL2GMap(A, rmapping, &is->rmapping, &localrmapping)); 2687 if (rmapping == cmapping) { 2688 PetscCall(PetscObjectReference((PetscObject)is->rmapping)); 2689 is->cmapping = is->rmapping; 2690 PetscCall(PetscObjectReference((PetscObject)localrmapping)); 2691 localcmapping = localrmapping; 2692 } 2693 } 2694 if (!cmapping) { 2695 IS is; 2696 2697 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), A->cmap->N, 0, 1, &is)); 2698 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmapping)); 2699 PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping, A->cmap->bs)); 2700 PetscCall(ISDestroy(&is)); 2701 freem[1] = PETSC_TRUE; 2702 } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */ 2703 PetscCall(MatISFilterL2GMap(A, cmapping, &is->cmapping, &localcmapping)); 2704 } 2705 if (!is->rmapping) { 2706 PetscCall(PetscObjectReference((PetscObject)rmapping)); 2707 is->rmapping = rmapping; 2708 } 2709 if (!is->cmapping) { 2710 PetscCall(PetscObjectReference((PetscObject)cmapping)); 2711 is->cmapping = cmapping; 2712 } 2713 2714 /* Clean up */ 2715 is->lnnzstate = 0; 2716 PetscCall(MatDestroy(&is->dA)); 2717 PetscCall(MatDestroy(&is->assembledA)); 2718 PetscCall(MatDestroy(&is->A)); 2719 if (is->csf != is->sf) { 2720 PetscCall(PetscSFDestroy(&is->csf)); 2721 PetscCall(PetscFree2(is->csf_rootdata, is->csf_leafdata)); 2722 } else is->csf = NULL; 2723 PetscCall(PetscSFDestroy(&is->sf)); 2724 PetscCall(PetscFree2(is->sf_rootdata, is->sf_leafdata)); 2725 PetscCall(PetscFree(is->bdiag)); 2726 2727 /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case 2728 (DOLFIN passes 2 different objects) */ 2729 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping, &nr)); 2730 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping, &rbs)); 2731 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping, &nc)); 2732 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping, &cbs)); 2733 if (is->rmapping != is->cmapping && cong) { 2734 PetscBool same = PETSC_FALSE; 2735 if (nr == nc && cbs == rbs) { 2736 const PetscInt *idxs1, *idxs2; 2737 2738 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping, &idxs1)); 2739 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping, &idxs2)); 2740 PetscCall(PetscArraycmp(idxs1, idxs2, nr / rbs, &same)); 2741 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping, &idxs1)); 2742 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping, &idxs2)); 2743 } 2744 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 2745 if (same) { 2746 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2747 PetscCall(PetscObjectReference((PetscObject)is->rmapping)); 2748 is->cmapping = is->rmapping; 2749 } 2750 } 2751 PetscCall(PetscLayoutSetBlockSize(A->rmap, rbs)); 2752 PetscCall(PetscLayoutSetBlockSize(A->cmap, cbs)); 2753 /* Pass the user defined maps to the layout */ 2754 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap, rmapping)); 2755 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap, cmapping)); 2756 if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping)); 2757 if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping)); 2758 2759 if (!is->islocalref) { 2760 /* Create the local matrix A */ 2761 PetscCall(MatCreate(PETSC_COMM_SELF, &is->A)); 2762 PetscCall(MatSetType(is->A, is->lmattype)); 2763 PetscCall(MatSetSizes(is->A, nr, nc, nr, nc)); 2764 PetscCall(MatSetBlockSizes(is->A, rbs, cbs)); 2765 PetscCall(MatSetOptionsPrefix(is->A, "is_")); 2766 PetscCall(MatAppendOptionsPrefix(is->A, ((PetscObject)A)->prefix)); 2767 PetscCall(PetscLayoutSetUp(is->A->rmap)); 2768 PetscCall(PetscLayoutSetUp(is->A->cmap)); 2769 PetscCall(MatSetLocalToGlobalMapping(is->A, localrmapping, localcmapping)); 2770 PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping)); 2771 PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping)); 2772 /* setup scatters and local vectors for MatMult */ 2773 PetscCall(MatISSetUpScatters_Private(A)); 2774 } 2775 A->preallocated = PETSC_TRUE; 2776 PetscFunctionReturn(PETSC_SUCCESS); 2777 } 2778 2779 static PetscErrorCode MatSetUp_IS(Mat A) 2780 { 2781 Mat_IS *is = (Mat_IS *)A->data; 2782 ISLocalToGlobalMapping rmap, cmap; 2783 2784 PetscFunctionBegin; 2785 if (!is->sf) { 2786 PetscCall(MatGetLocalToGlobalMapping(A, &rmap, &cmap)); 2787 PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap)); 2788 } 2789 PetscFunctionReturn(PETSC_SUCCESS); 2790 } 2791 2792 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2793 { 2794 Mat_IS *is = (Mat_IS *)mat->data; 2795 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2796 2797 PetscFunctionBegin; 2798 PetscCall(ISGlobalToLocalMappingApply(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l)); 2799 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2800 PetscCall(ISGlobalToLocalMappingApply(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l)); 2801 PetscCall(MatSetValues(is->A, m, rows_l, n, cols_l, values, addv)); 2802 } else { 2803 PetscCall(MatSetValues(is->A, m, rows_l, m, rows_l, values, addv)); 2804 } 2805 PetscFunctionReturn(PETSC_SUCCESS); 2806 } 2807 2808 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2809 { 2810 Mat_IS *is = (Mat_IS *)mat->data; 2811 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION], cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2812 2813 PetscFunctionBegin; 2814 PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping, IS_GTOLM_MASK, m, rows, &m, rows_l)); 2815 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2816 PetscCall(ISGlobalToLocalMappingApplyBlock(is->cmapping, IS_GTOLM_MASK, n, cols, &n, cols_l)); 2817 PetscCall(MatSetValuesBlocked(is->A, m, rows_l, n, cols_l, values, addv)); 2818 } else { 2819 PetscCall(MatSetValuesBlocked(is->A, m, rows_l, m, rows_l, values, addv)); 2820 } 2821 PetscFunctionReturn(PETSC_SUCCESS); 2822 } 2823 2824 static PetscErrorCode MatSetValuesLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2825 { 2826 Mat_IS *is = (Mat_IS *)A->data; 2827 2828 PetscFunctionBegin; 2829 if (is->A->rmap->mapping || is->A->cmap->mapping) { 2830 PetscCall(MatSetValuesLocal(is->A, m, rows, n, cols, values, addv)); 2831 } else { 2832 PetscCall(MatSetValues(is->A, m, rows, n, cols, values, addv)); 2833 } 2834 PetscFunctionReturn(PETSC_SUCCESS); 2835 } 2836 2837 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2838 { 2839 Mat_IS *is = (Mat_IS *)A->data; 2840 2841 PetscFunctionBegin; 2842 if (is->A->rmap->mapping || is->A->cmap->mapping) { 2843 PetscCall(MatSetValuesBlockedLocal(is->A, m, rows, n, cols, values, addv)); 2844 } else { 2845 PetscCall(MatSetValuesBlocked(is->A, m, rows, n, cols, values, addv)); 2846 } 2847 PetscFunctionReturn(PETSC_SUCCESS); 2848 } 2849 2850 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, PetscBool columns) 2851 { 2852 Mat_IS *is = (Mat_IS *)A->data; 2853 2854 PetscFunctionBegin; 2855 if (!n) PetscFunctionReturn(PETSC_SUCCESS); 2856 is->pure_neumann = PETSC_FALSE; 2857 2858 if (columns) { 2859 PetscCall(MatZeroRowsColumns(is->A, n, rows, diag, NULL, NULL)); 2860 } else { 2861 PetscCall(MatZeroRows(is->A, n, rows, diag, NULL, NULL)); 2862 } 2863 if (diag != 0.) { 2864 const PetscScalar *array; 2865 2866 PetscCall(VecGetArrayRead(is->counter, &array)); 2867 for (PetscInt i = 0; i < n; i++) PetscCall(MatSetValue(is->A, rows[i], rows[i], diag / (array[rows[i]]), INSERT_VALUES)); 2868 PetscCall(VecRestoreArrayRead(is->counter, &array)); 2869 PetscCall(MatAssemblyBegin(is->A, MAT_FINAL_ASSEMBLY)); 2870 PetscCall(MatAssemblyEnd(is->A, MAT_FINAL_ASSEMBLY)); 2871 } 2872 PetscFunctionReturn(PETSC_SUCCESS); 2873 } 2874 2875 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b, PetscBool columns) 2876 { 2877 Mat_IS *matis = (Mat_IS *)A->data; 2878 PetscInt nr, nl, len; 2879 PetscInt *lrows = NULL; 2880 2881 PetscFunctionBegin; 2882 if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) { 2883 PetscBool cong; 2884 2885 PetscCall(PetscLayoutCompare(A->rmap, A->cmap, &cong)); 2886 cong = (PetscBool)(cong && matis->sf == matis->csf); 2887 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"); 2888 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"); 2889 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"); 2890 } 2891 PetscCall(MatGetSize(matis->A, &nl, NULL)); 2892 /* get locally owned rows */ 2893 PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL)); 2894 /* fix right-hand side if needed */ 2895 if (x && b) { 2896 const PetscScalar *xx; 2897 PetscScalar *bb; 2898 2899 PetscCall(VecGetArrayRead(x, &xx)); 2900 PetscCall(VecGetArray(b, &bb)); 2901 for (PetscInt i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]]; 2902 PetscCall(VecRestoreArrayRead(x, &xx)); 2903 PetscCall(VecRestoreArray(b, &bb)); 2904 } 2905 /* get rows associated to the local matrices */ 2906 PetscCall(PetscArrayzero(matis->sf_leafdata, nl)); 2907 PetscCall(PetscArrayzero(matis->sf_rootdata, A->rmap->n)); 2908 for (PetscInt i = 0; i < len; i++) matis->sf_rootdata[lrows[i]] = 1; 2909 PetscCall(PetscFree(lrows)); 2910 PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 2911 PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE)); 2912 PetscCall(PetscMalloc1(nl, &lrows)); 2913 nr = 0; 2914 for (PetscInt i = 0; i < nl; i++) 2915 if (matis->sf_leafdata[i]) lrows[nr++] = i; 2916 PetscCall(MatISZeroRowsColumnsLocal_Private(A, nr, lrows, diag, columns)); 2917 PetscCall(PetscFree(lrows)); 2918 PetscFunctionReturn(PETSC_SUCCESS); 2919 } 2920 2921 static PetscErrorCode MatZeroRows_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 2922 { 2923 PetscFunctionBegin; 2924 PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_FALSE)); 2925 PetscFunctionReturn(PETSC_SUCCESS); 2926 } 2927 2928 static PetscErrorCode MatZeroRowsColumns_IS(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 2929 { 2930 PetscFunctionBegin; 2931 PetscCall(MatZeroRowsColumns_Private_IS(A, n, rows, diag, x, b, PETSC_TRUE)); 2932 PetscFunctionReturn(PETSC_SUCCESS); 2933 } 2934 2935 static PetscErrorCode MatAssemblyBegin_IS(Mat A, MatAssemblyType type) 2936 { 2937 Mat_IS *is = (Mat_IS *)A->data; 2938 2939 PetscFunctionBegin; 2940 PetscCall(MatAssemblyBegin(is->A, type)); 2941 PetscFunctionReturn(PETSC_SUCCESS); 2942 } 2943 2944 static PetscErrorCode MatAssemblyEnd_IS(Mat A, MatAssemblyType type) 2945 { 2946 Mat_IS *is = (Mat_IS *)A->data; 2947 PetscBool lnnz; 2948 2949 PetscFunctionBegin; 2950 PetscCall(MatAssemblyEnd(is->A, type)); 2951 /* fix for local empty rows/cols */ 2952 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2953 Mat newlA; 2954 ISLocalToGlobalMapping rl2g, cl2g; 2955 IS nzr, nzc; 2956 PetscInt nr, nc, nnzr, nnzc; 2957 PetscBool lnewl2g, newl2g; 2958 2959 PetscCall(MatGetSize(is->A, &nr, &nc)); 2960 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_FALSE, PETSC_SMALL, &nzr)); 2961 if (!nzr) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nr, 0, 1, &nzr)); 2962 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A, PETSC_TRUE, PETSC_SMALL, &nzc)); 2963 if (!nzc) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A), nc, 0, 1, &nzc)); 2964 PetscCall(ISGetSize(nzr, &nnzr)); 2965 PetscCall(ISGetSize(nzc, &nnzc)); 2966 if (nnzr != nr || nnzc != nc) { /* need new global l2g map */ 2967 lnewl2g = PETSC_TRUE; 2968 PetscCallMPI(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A))); 2969 2970 /* extract valid submatrix */ 2971 PetscCall(MatCreateSubMatrix(is->A, nzr, nzc, MAT_INITIAL_MATRIX, &newlA)); 2972 } else { /* local matrix fully populated */ 2973 lnewl2g = PETSC_FALSE; 2974 PetscCallMPI(MPIU_Allreduce(&lnewl2g, &newl2g, 1, MPI_C_BOOL, MPI_LOR, PetscObjectComm((PetscObject)A))); 2975 PetscCall(PetscObjectReference((PetscObject)is->A)); 2976 newlA = is->A; 2977 } 2978 2979 /* attach new global l2g map if needed */ 2980 if (newl2g) { 2981 IS zr, zc; 2982 const PetscInt *ridxs, *cidxs, *zridxs, *zcidxs; 2983 PetscInt *nidxs, i; 2984 2985 PetscCall(ISComplement(nzr, 0, nr, &zr)); 2986 PetscCall(ISComplement(nzc, 0, nc, &zc)); 2987 PetscCall(PetscMalloc1(PetscMax(nr, nc), &nidxs)); 2988 PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping, &ridxs)); 2989 PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping, &cidxs)); 2990 PetscCall(ISGetIndices(zr, &zridxs)); 2991 PetscCall(ISGetIndices(zc, &zcidxs)); 2992 PetscCall(ISGetLocalSize(zr, &nnzr)); 2993 PetscCall(ISGetLocalSize(zc, &nnzc)); 2994 2995 PetscCall(PetscArraycpy(nidxs, ridxs, nr)); 2996 for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1; 2997 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nr, nidxs, PETSC_COPY_VALUES, &rl2g)); 2998 PetscCall(PetscArraycpy(nidxs, cidxs, nc)); 2999 for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1; 3000 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, nc, nidxs, PETSC_COPY_VALUES, &cl2g)); 3001 3002 PetscCall(ISRestoreIndices(zr, &zridxs)); 3003 PetscCall(ISRestoreIndices(zc, &zcidxs)); 3004 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping, &ridxs)); 3005 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping, &cidxs)); 3006 PetscCall(ISDestroy(&nzr)); 3007 PetscCall(ISDestroy(&nzc)); 3008 PetscCall(ISDestroy(&zr)); 3009 PetscCall(ISDestroy(&zc)); 3010 PetscCall(PetscFree(nidxs)); 3011 PetscCall(MatSetLocalToGlobalMapping(A, rl2g, cl2g)); 3012 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 3013 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 3014 } 3015 PetscCall(MatISSetLocalMat(A, newlA)); 3016 PetscCall(MatDestroy(&newlA)); 3017 PetscCall(ISDestroy(&nzr)); 3018 PetscCall(ISDestroy(&nzc)); 3019 is->locempty = PETSC_FALSE; 3020 } 3021 lnnz = (PetscBool)(is->A->nonzerostate == is->lnnzstate); 3022 is->lnnzstate = is->A->nonzerostate; 3023 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &lnnz, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A))); 3024 if (!lnnz) A->nonzerostate++; 3025 PetscFunctionReturn(PETSC_SUCCESS); 3026 } 3027 3028 static PetscErrorCode MatISGetLocalMat_IS(Mat mat, Mat *local) 3029 { 3030 Mat_IS *is = (Mat_IS *)mat->data; 3031 3032 PetscFunctionBegin; 3033 *local = is->A; 3034 PetscFunctionReturn(PETSC_SUCCESS); 3035 } 3036 3037 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat, Mat *local) 3038 { 3039 PetscFunctionBegin; 3040 *local = NULL; 3041 PetscFunctionReturn(PETSC_SUCCESS); 3042 } 3043 3044 /*@ 3045 MatISGetLocalMat - Gets the local matrix stored inside a `MATIS` matrix. 3046 3047 Not Collective. 3048 3049 Input Parameter: 3050 . mat - the matrix 3051 3052 Output Parameter: 3053 . local - the local matrix 3054 3055 Level: intermediate 3056 3057 Notes: 3058 This can be called if you have precomputed the nonzero structure of the 3059 matrix and want to provide it to the inner matrix object to improve the performance 3060 of the `MatSetValues()` operation. 3061 3062 Call `MatISRestoreLocalMat()` when finished with the local matrix. 3063 3064 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISRestoreLocalMat()` 3065 @*/ 3066 PetscErrorCode MatISGetLocalMat(Mat mat, Mat *local) 3067 { 3068 PetscFunctionBegin; 3069 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3070 PetscAssertPointer(local, 2); 3071 PetscUseMethod(mat, "MatISGetLocalMat_C", (Mat, Mat *), (mat, local)); 3072 PetscFunctionReturn(PETSC_SUCCESS); 3073 } 3074 3075 /*@ 3076 MatISRestoreLocalMat - Restores the local matrix obtained with `MatISGetLocalMat()` 3077 3078 Not Collective. 3079 3080 Input Parameters: 3081 + mat - the matrix 3082 - local - the local matrix 3083 3084 Level: intermediate 3085 3086 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISGetLocalMat()` 3087 @*/ 3088 PetscErrorCode MatISRestoreLocalMat(Mat mat, Mat *local) 3089 { 3090 PetscFunctionBegin; 3091 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3092 PetscAssertPointer(local, 2); 3093 PetscUseMethod(mat, "MatISRestoreLocalMat_C", (Mat, Mat *), (mat, local)); 3094 PetscFunctionReturn(PETSC_SUCCESS); 3095 } 3096 3097 static PetscErrorCode MatISSetLocalMatType_IS(Mat mat, MatType mtype) 3098 { 3099 Mat_IS *is = (Mat_IS *)mat->data; 3100 3101 PetscFunctionBegin; 3102 if (is->A) PetscCall(MatSetType(is->A, mtype)); 3103 PetscCall(PetscFree(is->lmattype)); 3104 PetscCall(PetscStrallocpy(mtype, &is->lmattype)); 3105 PetscFunctionReturn(PETSC_SUCCESS); 3106 } 3107 3108 /*@ 3109 MatISSetLocalMatType - Specifies the type of local matrix inside the `MATIS` 3110 3111 Logically Collective. 3112 3113 Input Parameters: 3114 + mat - the matrix 3115 - mtype - the local matrix type 3116 3117 Level: intermediate 3118 3119 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetType()`, `MatType` 3120 @*/ 3121 PetscErrorCode MatISSetLocalMatType(Mat mat, MatType mtype) 3122 { 3123 PetscFunctionBegin; 3124 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3125 PetscUseMethod(mat, "MatISSetLocalMatType_C", (Mat, MatType), (mat, mtype)); 3126 PetscFunctionReturn(PETSC_SUCCESS); 3127 } 3128 3129 static PetscErrorCode MatISSetLocalMat_IS(Mat mat, Mat local) 3130 { 3131 Mat_IS *is = (Mat_IS *)mat->data; 3132 PetscInt nrows, ncols, orows, ocols; 3133 MatType mtype, otype; 3134 PetscBool sametype = PETSC_TRUE; 3135 3136 PetscFunctionBegin; 3137 if (is->A && !is->islocalref) { 3138 PetscCall(MatGetSize(is->A, &orows, &ocols)); 3139 PetscCall(MatGetSize(local, &nrows, &ncols)); 3140 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); 3141 PetscCall(MatGetType(local, &mtype)); 3142 PetscCall(MatGetType(is->A, &otype)); 3143 PetscCall(PetscStrcmp(mtype, otype, &sametype)); 3144 } 3145 PetscCall(PetscObjectReference((PetscObject)local)); 3146 PetscCall(MatDestroy(&is->A)); 3147 is->A = local; 3148 PetscCall(MatGetType(is->A, &mtype)); 3149 PetscCall(MatISSetLocalMatType(mat, mtype)); 3150 if (!sametype && !is->islocalref) PetscCall(MatISSetUpScatters_Private(mat)); 3151 is->lnnzstate = 0; 3152 PetscFunctionReturn(PETSC_SUCCESS); 3153 } 3154 3155 /*@ 3156 MatISSetLocalMat - Replace the local matrix stored inside a `MATIS` object. 3157 3158 Not Collective 3159 3160 Input Parameters: 3161 + mat - the matrix 3162 - local - the local matrix 3163 3164 Level: intermediate 3165 3166 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatISSetLocalMatType`, `MatISGetLocalMat()` 3167 @*/ 3168 PetscErrorCode MatISSetLocalMat(Mat mat, Mat local) 3169 { 3170 PetscFunctionBegin; 3171 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 3172 PetscValidHeaderSpecific(local, MAT_CLASSID, 2); 3173 PetscUseMethod(mat, "MatISSetLocalMat_C", (Mat, Mat), (mat, local)); 3174 PetscFunctionReturn(PETSC_SUCCESS); 3175 } 3176 3177 static PetscErrorCode MatZeroEntries_IS(Mat A) 3178 { 3179 Mat_IS *a = (Mat_IS *)A->data; 3180 3181 PetscFunctionBegin; 3182 PetscCall(MatZeroEntries(a->A)); 3183 PetscFunctionReturn(PETSC_SUCCESS); 3184 } 3185 3186 static PetscErrorCode MatScale_IS(Mat A, PetscScalar a) 3187 { 3188 Mat_IS *is = (Mat_IS *)A->data; 3189 3190 PetscFunctionBegin; 3191 PetscCall(MatScale(is->A, a)); 3192 PetscFunctionReturn(PETSC_SUCCESS); 3193 } 3194 3195 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 3196 { 3197 Mat_IS *is = (Mat_IS *)A->data; 3198 3199 PetscFunctionBegin; 3200 /* get diagonal of the local matrix */ 3201 PetscCall(MatGetDiagonal(is->A, is->y)); 3202 3203 /* scatter diagonal back into global vector */ 3204 PetscCall(VecSet(v, 0)); 3205 PetscCall(VecScatterBegin(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE)); 3206 PetscCall(VecScatterEnd(is->rctx, is->y, v, ADD_VALUES, SCATTER_REVERSE)); 3207 PetscFunctionReturn(PETSC_SUCCESS); 3208 } 3209 3210 static PetscErrorCode MatSetOption_IS(Mat A, MatOption op, PetscBool flg) 3211 { 3212 Mat_IS *a = (Mat_IS *)A->data; 3213 3214 PetscFunctionBegin; 3215 PetscCall(MatSetOption(a->A, op, flg)); 3216 PetscFunctionReturn(PETSC_SUCCESS); 3217 } 3218 3219 static PetscErrorCode MatAXPY_IS(Mat Y, PetscScalar a, Mat X, MatStructure str) 3220 { 3221 Mat_IS *y = (Mat_IS *)Y->data; 3222 Mat_IS *x; 3223 3224 PetscFunctionBegin; 3225 if (PetscDefined(USE_DEBUG)) { 3226 PetscBool ismatis; 3227 PetscCall(PetscObjectTypeCompare((PetscObject)X, MATIS, &ismatis)); 3228 PetscCheck(ismatis, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 3229 } 3230 x = (Mat_IS *)X->data; 3231 PetscCall(MatAXPY(y->A, a, x->A, str)); 3232 PetscFunctionReturn(PETSC_SUCCESS); 3233 } 3234 3235 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A, IS row, IS col, Mat *submat) 3236 { 3237 Mat lA; 3238 Mat_IS *matis = (Mat_IS *)A->data; 3239 ISLocalToGlobalMapping rl2g, cl2g; 3240 IS is; 3241 const PetscInt *rg, *rl; 3242 PetscInt nrg; 3243 PetscInt N, M, nrl, i, *idxs; 3244 3245 PetscFunctionBegin; 3246 PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping, &rg)); 3247 PetscCall(ISGetLocalSize(row, &nrl)); 3248 PetscCall(ISGetIndices(row, &rl)); 3249 PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping, &nrg)); 3250 if (PetscDefined(USE_DEBUG)) { 3251 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); 3252 } 3253 PetscCall(PetscMalloc1(nrg, &idxs)); 3254 /* map from [0,nrl) to row */ 3255 for (i = 0; i < nrl; i++) idxs[i] = rl[i]; 3256 for (i = nrl; i < nrg; i++) idxs[i] = -1; 3257 PetscCall(ISRestoreIndices(row, &rl)); 3258 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping, &rg)); 3259 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), nrg, idxs, PETSC_OWN_POINTER, &is)); 3260 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g)); 3261 PetscCall(ISDestroy(&is)); 3262 /* compute new l2g map for columns */ 3263 if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) { 3264 const PetscInt *cg, *cl; 3265 PetscInt ncg; 3266 PetscInt ncl; 3267 3268 PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping, &cg)); 3269 PetscCall(ISGetLocalSize(col, &ncl)); 3270 PetscCall(ISGetIndices(col, &cl)); 3271 PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping, &ncg)); 3272 if (PetscDefined(USE_DEBUG)) { 3273 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); 3274 } 3275 PetscCall(PetscMalloc1(ncg, &idxs)); 3276 /* map from [0,ncl) to col */ 3277 for (i = 0; i < ncl; i++) idxs[i] = cl[i]; 3278 for (i = ncl; i < ncg; i++) idxs[i] = -1; 3279 PetscCall(ISRestoreIndices(col, &cl)); 3280 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping, &cg)); 3281 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A), ncg, idxs, PETSC_OWN_POINTER, &is)); 3282 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g)); 3283 PetscCall(ISDestroy(&is)); 3284 } else { 3285 PetscCall(PetscObjectReference((PetscObject)rl2g)); 3286 cl2g = rl2g; 3287 } 3288 3289 /* create the MATIS submatrix */ 3290 PetscCall(MatGetSize(A, &M, &N)); 3291 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), submat)); 3292 PetscCall(MatSetSizes(*submat, PETSC_DECIDE, PETSC_DECIDE, M, N)); 3293 PetscCall(MatSetType(*submat, MATIS)); 3294 matis = (Mat_IS *)((*submat)->data); 3295 matis->islocalref = A; 3296 PetscCall(MatSetLocalToGlobalMapping(*submat, rl2g, cl2g)); 3297 PetscCall(MatISGetLocalMat(A, &lA)); 3298 PetscCall(MatISSetLocalMat(*submat, lA)); 3299 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 3300 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 3301 3302 /* remove unsupported ops */ 3303 PetscCall(PetscMemzero((*submat)->ops, sizeof(struct _MatOps))); 3304 (*submat)->ops->destroy = MatDestroy_IS; 3305 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 3306 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 3307 (*submat)->ops->zerorowslocal = MatZeroRowsLocal_SubMat_IS; 3308 (*submat)->ops->zerorowscolumnslocal = MatZeroRowsColumnsLocal_SubMat_IS; 3309 (*submat)->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 3310 PetscFunctionReturn(PETSC_SUCCESS); 3311 } 3312 3313 static PetscErrorCode MatSetFromOptions_IS(Mat A, PetscOptionItems PetscOptionsObject) 3314 { 3315 Mat_IS *a = (Mat_IS *)A->data; 3316 char type[256]; 3317 PetscBool flg; 3318 3319 PetscFunctionBegin; 3320 PetscOptionsHeadBegin(PetscOptionsObject, "MATIS options"); 3321 PetscCall(PetscOptionsDeprecated("-matis_keepassembled", "-mat_is_keepassembled", "3.21", NULL)); 3322 PetscCall(PetscOptionsDeprecated("-matis_fixempty", "-mat_is_fixempty", "3.21", NULL)); 3323 PetscCall(PetscOptionsDeprecated("-matis_storel2l", "-mat_is_storel2l", "3.21", NULL)); 3324 PetscCall(PetscOptionsDeprecated("-matis_localmat_type", "-mat_is_localmat_type", "3.21", NULL)); 3325 PetscCall(PetscOptionsBool("-mat_is_keepassembled", "Store an assembled version if needed", NULL, a->keepassembled, &a->keepassembled, NULL)); 3326 PetscCall(PetscOptionsBool("-mat_is_fixempty", "Fix local matrices in case of empty local rows/columns", "MatISFixLocalEmpty", a->locempty, &a->locempty, NULL)); 3327 PetscCall(PetscOptionsBool("-mat_is_storel2l", "Store local-to-local matrices generated from PtAP operations", "MatISStoreL2L", a->storel2l, &a->storel2l, NULL)); 3328 PetscCall(PetscOptionsBool("-mat_is_allow_repeated", "Allow local repeated entries", "MatISSetAllowRepeated", a->allow_repeated, &a->allow_repeated, NULL)); 3329 PetscCall(PetscOptionsFList("-mat_is_localmat_type", "Matrix type", "MatISSetLocalMatType", MatList, a->lmattype, type, 256, &flg)); 3330 if (flg) PetscCall(MatISSetLocalMatType(A, type)); 3331 if (a->A) PetscCall(MatSetFromOptions(a->A)); 3332 PetscOptionsHeadEnd(); 3333 PetscFunctionReturn(PETSC_SUCCESS); 3334 } 3335 3336 /*@ 3337 MatCreateIS - Creates a "process" unassembled matrix. 3338 3339 Collective. 3340 3341 Input Parameters: 3342 + comm - MPI communicator that will share the matrix 3343 . bs - block size of the matrix 3344 . m - local size of left vector used in matrix vector products 3345 . n - local size of right vector used in matrix vector products 3346 . M - global size of left vector used in matrix vector products 3347 . N - global size of right vector used in matrix vector products 3348 . rmap - local to global map for rows 3349 - cmap - local to global map for cols 3350 3351 Output Parameter: 3352 . A - the resulting matrix 3353 3354 Level: intermediate 3355 3356 Notes: 3357 `m` and `n` are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors 3358 used in `MatMult()` operations. The local sizes of `rmap` and `cmap` define the size of the local matrices. 3359 3360 If `rmap` (`cmap`) is `NULL`, then the local row (column) spaces matches the global space. 3361 3362 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()` 3363 @*/ 3364 PetscErrorCode MatCreateIS(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt M, PetscInt N, ISLocalToGlobalMapping rmap, ISLocalToGlobalMapping cmap, Mat *A) 3365 { 3366 PetscFunctionBegin; 3367 PetscCall(MatCreate(comm, A)); 3368 PetscCall(MatSetSizes(*A, m, n, M, N)); 3369 if (bs > 0) PetscCall(MatSetBlockSize(*A, bs)); 3370 PetscCall(MatSetType(*A, MATIS)); 3371 PetscCall(MatSetLocalToGlobalMapping(*A, rmap, cmap)); 3372 PetscFunctionReturn(PETSC_SUCCESS); 3373 } 3374 3375 static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has) 3376 { 3377 Mat_IS *a = (Mat_IS *)A->data; 3378 MatOperation tobefiltered[] = {MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP}; 3379 3380 PetscFunctionBegin; 3381 *has = PETSC_FALSE; 3382 if (!((void **)A->ops)[op] || !a->A) PetscFunctionReturn(PETSC_SUCCESS); 3383 *has = PETSC_TRUE; 3384 for (PetscInt i = 0; i < (PetscInt)PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++) 3385 if (op == tobefiltered[i]) PetscFunctionReturn(PETSC_SUCCESS); 3386 PetscCall(MatHasOperation(a->A, op, has)); 3387 PetscFunctionReturn(PETSC_SUCCESS); 3388 } 3389 3390 static PetscErrorCode MatSetValuesCOO_IS(Mat A, const PetscScalar v[], InsertMode imode) 3391 { 3392 Mat_IS *a = (Mat_IS *)A->data; 3393 3394 PetscFunctionBegin; 3395 PetscCall(MatSetValuesCOO(a->A, v, imode)); 3396 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 3397 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 3398 PetscFunctionReturn(PETSC_SUCCESS); 3399 } 3400 3401 static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[]) 3402 { 3403 Mat_IS *a = (Mat_IS *)A->data; 3404 3405 PetscFunctionBegin; 3406 PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping"); 3407 if (a->A->rmap->mapping || a->A->cmap->mapping) { 3408 PetscCall(MatSetPreallocationCOOLocal(a->A, ncoo, coo_i, coo_j)); 3409 } else { 3410 PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j)); 3411 } 3412 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS)); 3413 A->preallocated = PETSC_TRUE; 3414 PetscFunctionReturn(PETSC_SUCCESS); 3415 } 3416 3417 static PetscErrorCode MatSetPreallocationCOO_IS(Mat A, PetscCount ncoo, PetscInt coo_i[], PetscInt coo_j[]) 3418 { 3419 Mat_IS *a = (Mat_IS *)A->data; 3420 PetscInt ncoo_i; 3421 3422 PetscFunctionBegin; 3423 PetscCheck(a->A, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to provide l2g map first via MatSetLocalToGlobalMapping"); 3424 PetscCall(PetscIntCast(ncoo, &ncoo_i)); 3425 PetscCall(ISGlobalToLocalMappingApply(a->rmapping, IS_GTOLM_MASK, ncoo_i, coo_i, NULL, coo_i)); 3426 PetscCall(ISGlobalToLocalMappingApply(a->cmapping, IS_GTOLM_MASK, ncoo_i, coo_j, NULL, coo_j)); 3427 PetscCall(MatSetPreallocationCOO(a->A, ncoo, coo_i, coo_j)); 3428 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", MatSetValuesCOO_IS)); 3429 A->preallocated = PETSC_TRUE; 3430 PetscFunctionReturn(PETSC_SUCCESS); 3431 } 3432 3433 static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA) 3434 { 3435 Mat_IS *a = (Mat_IS *)A->data; 3436 PetscObjectState Astate, aAstate = PETSC_INT_MIN; 3437 PetscObjectState Annzstate, aAnnzstate = PETSC_INT_MIN; 3438 3439 PetscFunctionBegin; 3440 PetscCall(PetscObjectStateGet((PetscObject)A, &Astate)); 3441 Annzstate = A->nonzerostate; 3442 if (a->assembledA) { 3443 PetscCall(PetscObjectStateGet((PetscObject)a->assembledA, &aAstate)); 3444 aAnnzstate = a->assembledA->nonzerostate; 3445 } 3446 if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA)); 3447 if (Astate != aAstate || !a->assembledA) { 3448 MatType aAtype; 3449 PetscMPIInt size; 3450 PetscInt rbs, cbs, bs; 3451 3452 /* the assembled form is used as temporary storage for parallel operations 3453 like createsubmatrices and the like, do not waste device memory */ 3454 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 3455 PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping, &cbs)); 3456 PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping, &rbs)); 3457 bs = rbs == cbs ? rbs : 1; 3458 if (a->assembledA) PetscCall(MatGetType(a->assembledA, &aAtype)); 3459 else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ; 3460 else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ; 3461 3462 PetscCall(MatConvert(A, aAtype, a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &a->assembledA)); 3463 PetscCall(PetscObjectStateSet((PetscObject)a->assembledA, Astate)); 3464 a->assembledA->nonzerostate = Annzstate; 3465 } 3466 PetscCall(PetscObjectReference((PetscObject)a->assembledA)); 3467 *tA = a->assembledA; 3468 if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA)); 3469 PetscFunctionReturn(PETSC_SUCCESS); 3470 } 3471 3472 static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA) 3473 { 3474 PetscFunctionBegin; 3475 PetscCall(MatDestroy(tA)); 3476 *tA = NULL; 3477 PetscFunctionReturn(PETSC_SUCCESS); 3478 } 3479 3480 static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA) 3481 { 3482 Mat_IS *a = (Mat_IS *)A->data; 3483 PetscObjectState Astate, dAstate = PETSC_INT_MIN; 3484 3485 PetscFunctionBegin; 3486 PetscCall(PetscObjectStateGet((PetscObject)A, &Astate)); 3487 if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA, &dAstate)); 3488 if (Astate != dAstate) { 3489 Mat tA; 3490 MatType ltype; 3491 3492 PetscCall(MatDestroy(&a->dA)); 3493 PetscCall(MatISGetAssembled_Private(A, &tA)); 3494 PetscCall(MatGetDiagonalBlock(tA, &a->dA)); 3495 PetscCall(MatPropagateSymmetryOptions(tA, a->dA)); 3496 PetscCall(MatGetType(a->A, <ype)); 3497 PetscCall(MatConvert(a->dA, ltype, MAT_INPLACE_MATRIX, &a->dA)); 3498 PetscCall(PetscObjectReference((PetscObject)a->dA)); 3499 PetscCall(MatISRestoreAssembled_Private(A, &tA)); 3500 PetscCall(PetscObjectStateSet((PetscObject)a->dA, Astate)); 3501 } 3502 *dA = a->dA; 3503 PetscFunctionReturn(PETSC_SUCCESS); 3504 } 3505 3506 static PetscErrorCode MatCreateSubMatrices_IS(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse reuse, Mat *submat[]) 3507 { 3508 Mat tA; 3509 3510 PetscFunctionBegin; 3511 PetscCall(MatISGetAssembled_Private(A, &tA)); 3512 PetscCall(MatCreateSubMatrices(tA, n, irow, icol, reuse, submat)); 3513 /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */ 3514 #if 0 3515 { 3516 Mat_IS *a = (Mat_IS*)A->data; 3517 MatType ltype; 3518 VecType vtype; 3519 char *flg; 3520 3521 PetscCall(MatGetType(a->A,<ype)); 3522 PetscCall(MatGetVecType(a->A,&vtype)); 3523 PetscCall(PetscStrstr(vtype,"cuda",&flg)); 3524 if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg)); 3525 if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg)); 3526 if (flg) { 3527 for (PetscInt i = 0; i < n; i++) { 3528 Mat sA = (*submat)[i]; 3529 3530 PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA)); 3531 (*submat)[i] = sA; 3532 } 3533 } 3534 } 3535 #endif 3536 PetscCall(MatISRestoreAssembled_Private(A, &tA)); 3537 PetscFunctionReturn(PETSC_SUCCESS); 3538 } 3539 3540 static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov) 3541 { 3542 Mat tA; 3543 3544 PetscFunctionBegin; 3545 PetscCall(MatISGetAssembled_Private(A, &tA)); 3546 PetscCall(MatIncreaseOverlap(tA, n, is, ov)); 3547 PetscCall(MatISRestoreAssembled_Private(A, &tA)); 3548 PetscFunctionReturn(PETSC_SUCCESS); 3549 } 3550 3551 /*@ 3552 MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the `MATIS` object 3553 3554 Not Collective 3555 3556 Input Parameter: 3557 . A - the matrix 3558 3559 Output Parameters: 3560 + rmapping - row mapping 3561 - cmapping - column mapping 3562 3563 Level: advanced 3564 3565 Note: 3566 The returned map can be different from the one used to construct the `MATIS` object, since it will not contain negative or repeated indices. 3567 3568 .seealso: [](ch_matrices), `Mat`, `MATIS`, `MatSetLocalToGlobalMapping()` 3569 @*/ 3570 PetscErrorCode MatISGetLocalToGlobalMapping(Mat A, ISLocalToGlobalMapping *rmapping, ISLocalToGlobalMapping *cmapping) 3571 { 3572 PetscFunctionBegin; 3573 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3574 PetscValidType(A, 1); 3575 if (rmapping) PetscAssertPointer(rmapping, 2); 3576 if (cmapping) PetscAssertPointer(cmapping, 3); 3577 PetscUseMethod(A, "MatISGetLocalToGlobalMapping_C", (Mat, ISLocalToGlobalMapping *, ISLocalToGlobalMapping *), (A, rmapping, cmapping)); 3578 PetscFunctionReturn(PETSC_SUCCESS); 3579 } 3580 3581 static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A, ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c) 3582 { 3583 Mat_IS *a = (Mat_IS *)A->data; 3584 3585 PetscFunctionBegin; 3586 if (r) *r = a->rmapping; 3587 if (c) *c = a->cmapping; 3588 PetscFunctionReturn(PETSC_SUCCESS); 3589 } 3590 3591 static PetscErrorCode MatSetBlockSizes_IS(Mat A, PetscInt rbs, PetscInt cbs) 3592 { 3593 Mat_IS *a = (Mat_IS *)A->data; 3594 3595 PetscFunctionBegin; 3596 if (a->A) PetscCall(MatSetBlockSizes(a->A, rbs, cbs)); 3597 PetscFunctionReturn(PETSC_SUCCESS); 3598 } 3599 3600 /*MC 3601 MATIS - MATIS = "is" - A matrix type to be used for non-overlapping domain decomposition methods (e.g. `PCBDDC` or `KSPFETIDP`). 3602 This stores the matrices in globally unassembled form and the parallel matrix vector product is handled "implicitly". 3603 3604 Options Database Keys: 3605 + -mat_type is - Set the matrix type to `MATIS`. 3606 . -mat_is_allow_repeated - Allow repeated entries in the local part of the local to global maps. 3607 . -mat_is_fixempty - Fix local matrices in case of empty local rows/columns. 3608 - -mat_is_storel2l - Store the local-to-local operators generated by the Galerkin process of `MatPtAP()`. 3609 3610 Level: intermediate 3611 3612 Notes: 3613 Options prefix for the inner matrix are given by `-is_mat_xxx` 3614 3615 You must call `MatSetLocalToGlobalMapping()` before using this matrix type. 3616 3617 You can do matrix preallocation on the local matrix after you obtain it with 3618 `MatISGetLocalMat()`; otherwise, you could use `MatISSetPreallocation()` or `MatXAIJSetPreallocation()` 3619 3620 .seealso: [](ch_matrices), `Mat`, `MATIS`, `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP` 3621 M*/ 3622 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 3623 { 3624 Mat_IS *a; 3625 3626 PetscFunctionBegin; 3627 PetscCall(PetscNew(&a)); 3628 PetscCall(PetscStrallocpy(MATAIJ, &a->lmattype)); 3629 A->data = (void *)a; 3630 3631 /* matrix ops */ 3632 PetscCall(PetscMemzero(A->ops, sizeof(struct _MatOps))); 3633 A->ops->mult = MatMult_IS; 3634 A->ops->multadd = MatMultAdd_IS; 3635 A->ops->multtranspose = MatMultTranspose_IS; 3636 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 3637 A->ops->destroy = MatDestroy_IS; 3638 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 3639 A->ops->setvalues = MatSetValues_IS; 3640 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 3641 A->ops->setvalueslocal = MatSetValuesLocal_IS; 3642 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 3643 A->ops->zerorows = MatZeroRows_IS; 3644 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 3645 A->ops->assemblybegin = MatAssemblyBegin_IS; 3646 A->ops->assemblyend = MatAssemblyEnd_IS; 3647 A->ops->view = MatView_IS; 3648 A->ops->load = MatLoad_IS; 3649 A->ops->zeroentries = MatZeroEntries_IS; 3650 A->ops->scale = MatScale_IS; 3651 A->ops->getdiagonal = MatGetDiagonal_IS; 3652 A->ops->setoption = MatSetOption_IS; 3653 A->ops->ishermitian = MatIsHermitian_IS; 3654 A->ops->issymmetric = MatIsSymmetric_IS; 3655 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 3656 A->ops->duplicate = MatDuplicate_IS; 3657 A->ops->missingdiagonal = MatMissingDiagonal_IS; 3658 A->ops->copy = MatCopy_IS; 3659 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 3660 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 3661 A->ops->axpy = MatAXPY_IS; 3662 A->ops->diagonalset = MatDiagonalSet_IS; 3663 A->ops->shift = MatShift_IS; 3664 A->ops->transpose = MatTranspose_IS; 3665 A->ops->getinfo = MatGetInfo_IS; 3666 A->ops->diagonalscale = MatDiagonalScale_IS; 3667 A->ops->setfromoptions = MatSetFromOptions_IS; 3668 A->ops->setup = MatSetUp_IS; 3669 A->ops->hasoperation = MatHasOperation_IS; 3670 A->ops->getdiagonalblock = MatGetDiagonalBlock_IS; 3671 A->ops->createsubmatrices = MatCreateSubMatrices_IS; 3672 A->ops->increaseoverlap = MatIncreaseOverlap_IS; 3673 A->ops->setblocksizes = MatSetBlockSizes_IS; 3674 3675 /* special MATIS functions */ 3676 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMatType_C", MatISSetLocalMatType_IS)); 3677 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalMat_C", MatISGetLocalMat_IS)); 3678 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISRestoreLocalMat_C", MatISRestoreLocalMat_IS)); 3679 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetLocalMat_C", MatISSetLocalMat_IS)); 3680 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetPreallocation_C", MatISSetPreallocation_IS)); 3681 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISSetAllowRepeated_C", MatISSetAllowRepeated_IS)); 3682 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISStoreL2L_C", MatISStoreL2L_IS)); 3683 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISFixLocalEmpty_C", MatISFixLocalEmpty_IS)); 3684 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatISGetLocalToGlobalMapping_C", MatISGetLocalToGlobalMapping_IS)); 3685 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpiaij_C", MatConvert_IS_XAIJ)); 3686 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpibaij_C", MatConvert_IS_XAIJ)); 3687 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_mpisbaij_C", MatConvert_IS_XAIJ)); 3688 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqaij_C", MatConvert_IS_XAIJ)); 3689 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqbaij_C", MatConvert_IS_XAIJ)); 3690 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_seqsbaij_C", MatConvert_IS_XAIJ)); 3691 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_is_aij_C", MatConvert_IS_XAIJ)); 3692 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOOLocal_C", MatSetPreallocationCOOLocal_IS)); 3693 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_IS)); 3694 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATIS)); 3695 PetscFunctionReturn(PETSC_SUCCESS); 3696 } 3697