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