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