1 #include <../src/mat/impls/nest/matnestimpl.h> /*I "petscmat.h" I*/ 2 #include <../src/mat/impls/aij/seq/aij.h> 3 #include <petscsf.h> 4 5 static PetscErrorCode MatSetUp_NestIS_Private(Mat, PetscInt, const IS[], PetscInt, const IS[]); 6 static PetscErrorCode MatCreateVecs_Nest(Mat, Vec *, Vec *); 7 static PetscErrorCode MatReset_Nest(Mat); 8 9 PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat, MatType, MatReuse, Mat *); 10 11 /* private functions */ 12 static PetscErrorCode MatNestGetSizes_Private(Mat A, PetscInt *m, PetscInt *n, PetscInt *M, PetscInt *N) 13 { 14 Mat_Nest *bA = (Mat_Nest *)A->data; 15 PetscInt i, j; 16 17 PetscFunctionBegin; 18 *m = *n = *M = *N = 0; 19 for (i = 0; i < bA->nr; i++) { /* rows */ 20 PetscInt sm, sM; 21 PetscCall(ISGetLocalSize(bA->isglobal.row[i], &sm)); 22 PetscCall(ISGetSize(bA->isglobal.row[i], &sM)); 23 *m += sm; 24 *M += sM; 25 } 26 for (j = 0; j < bA->nc; j++) { /* cols */ 27 PetscInt sn, sN; 28 PetscCall(ISGetLocalSize(bA->isglobal.col[j], &sn)); 29 PetscCall(ISGetSize(bA->isglobal.col[j], &sN)); 30 *n += sn; 31 *N += sN; 32 } 33 PetscFunctionReturn(PETSC_SUCCESS); 34 } 35 36 /* operations */ 37 static PetscErrorCode MatMult_Nest(Mat A, Vec x, Vec y) 38 { 39 Mat_Nest *bA = (Mat_Nest *)A->data; 40 Vec *bx = bA->right, *by = bA->left; 41 PetscInt i, j, nr = bA->nr, nc = bA->nc; 42 43 PetscFunctionBegin; 44 for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(y, bA->isglobal.row[i], &by[i])); 45 for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(x, bA->isglobal.col[i], &bx[i])); 46 for (i = 0; i < nr; i++) { 47 PetscCall(VecZeroEntries(by[i])); 48 for (j = 0; j < nc; j++) { 49 if (!bA->m[i][j]) continue; 50 /* y[i] <- y[i] + A[i][j] * x[j] */ 51 PetscCall(MatMultAdd(bA->m[i][j], bx[j], by[i], by[i])); 52 } 53 } 54 for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(y, bA->isglobal.row[i], &by[i])); 55 for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.col[i], &bx[i])); 56 PetscFunctionReturn(PETSC_SUCCESS); 57 } 58 59 static PetscErrorCode MatMultAdd_Nest(Mat A, Vec x, Vec y, Vec z) 60 { 61 Mat_Nest *bA = (Mat_Nest *)A->data; 62 Vec *bx = bA->right, *bz = bA->left; 63 PetscInt i, j, nr = bA->nr, nc = bA->nc; 64 65 PetscFunctionBegin; 66 for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(z, bA->isglobal.row[i], &bz[i])); 67 for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(x, bA->isglobal.col[i], &bx[i])); 68 for (i = 0; i < nr; i++) { 69 if (y != z) { 70 Vec by; 71 PetscCall(VecGetSubVector(y, bA->isglobal.row[i], &by)); 72 PetscCall(VecCopy(by, bz[i])); 73 PetscCall(VecRestoreSubVector(y, bA->isglobal.row[i], &by)); 74 } 75 for (j = 0; j < nc; j++) { 76 if (!bA->m[i][j]) continue; 77 /* y[i] <- y[i] + A[i][j] * x[j] */ 78 PetscCall(MatMultAdd(bA->m[i][j], bx[j], bz[i], bz[i])); 79 } 80 } 81 for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(z, bA->isglobal.row[i], &bz[i])); 82 for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.col[i], &bx[i])); 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 86 typedef struct { 87 Mat *workC; /* array of Mat with specific containers depending on the underlying MatMatMult implementation */ 88 PetscScalar *tarray; /* buffer for storing all temporary products A[i][j] B[j] */ 89 PetscInt *dm, *dn, k; /* displacements and number of submatrices */ 90 } Nest_Dense; 91 92 static PetscErrorCode MatProductNumeric_Nest_Dense(Mat C) 93 { 94 Mat_Nest *bA; 95 Nest_Dense *contents; 96 Mat viewB, viewC, productB, workC; 97 const PetscScalar *barray; 98 PetscScalar *carray; 99 PetscInt i, j, M, N, nr, nc, ldb, ldc; 100 Mat A, B; 101 102 PetscFunctionBegin; 103 MatCheckProduct(C, 1); 104 A = C->product->A; 105 B = C->product->B; 106 PetscCall(MatGetSize(B, NULL, &N)); 107 if (!N) { 108 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 109 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 110 PetscFunctionReturn(PETSC_SUCCESS); 111 } 112 contents = (Nest_Dense *)C->product->data; 113 PetscCheck(contents, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty"); 114 bA = (Mat_Nest *)A->data; 115 nr = bA->nr; 116 nc = bA->nc; 117 PetscCall(MatDenseGetLDA(B, &ldb)); 118 PetscCall(MatDenseGetLDA(C, &ldc)); 119 PetscCall(MatZeroEntries(C)); 120 PetscCall(MatDenseGetArrayRead(B, &barray)); 121 PetscCall(MatDenseGetArray(C, &carray)); 122 for (i = 0; i < nr; i++) { 123 PetscCall(ISGetSize(bA->isglobal.row[i], &M)); 124 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), contents->dm[i + 1] - contents->dm[i], PETSC_DECIDE, M, N, carray + contents->dm[i], &viewC)); 125 PetscCall(MatDenseSetLDA(viewC, ldc)); 126 for (j = 0; j < nc; j++) { 127 if (!bA->m[i][j]) continue; 128 PetscCall(ISGetSize(bA->isglobal.col[j], &M)); 129 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), contents->dn[j + 1] - contents->dn[j], PETSC_DECIDE, M, N, (PetscScalar *)(barray + contents->dn[j]), &viewB)); 130 PetscCall(MatDenseSetLDA(viewB, ldb)); 131 132 /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */ 133 workC = contents->workC[i * nc + j]; 134 productB = workC->product->B; 135 workC->product->B = viewB; /* use newly created dense matrix viewB */ 136 PetscCall(MatProductNumeric(workC)); 137 PetscCall(MatDestroy(&viewB)); 138 workC->product->B = productB; /* resume original B */ 139 140 /* C[i] <- workC + C[i] */ 141 PetscCall(MatAXPY(viewC, 1.0, contents->workC[i * nc + j], SAME_NONZERO_PATTERN)); 142 } 143 PetscCall(MatDestroy(&viewC)); 144 } 145 PetscCall(MatDenseRestoreArray(C, &carray)); 146 PetscCall(MatDenseRestoreArrayRead(B, &barray)); 147 148 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 149 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 150 PetscFunctionReturn(PETSC_SUCCESS); 151 } 152 153 static PetscErrorCode MatNest_DenseDestroy(void *ctx) 154 { 155 Nest_Dense *contents = (Nest_Dense *)ctx; 156 PetscInt i; 157 158 PetscFunctionBegin; 159 PetscCall(PetscFree(contents->tarray)); 160 for (i = 0; i < contents->k; i++) PetscCall(MatDestroy(contents->workC + i)); 161 PetscCall(PetscFree3(contents->dm, contents->dn, contents->workC)); 162 PetscCall(PetscFree(contents)); 163 PetscFunctionReturn(PETSC_SUCCESS); 164 } 165 166 static PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C) 167 { 168 Mat_Nest *bA; 169 Mat viewB, workC; 170 const PetscScalar *barray; 171 PetscInt i, j, M, N, m, n, nr, nc, maxm = 0, ldb; 172 Nest_Dense *contents = NULL; 173 PetscBool cisdense; 174 Mat A, B; 175 PetscReal fill; 176 177 PetscFunctionBegin; 178 MatCheckProduct(C, 1); 179 PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty"); 180 A = C->product->A; 181 B = C->product->B; 182 fill = C->product->fill; 183 bA = (Mat_Nest *)A->data; 184 nr = bA->nr; 185 nc = bA->nc; 186 PetscCall(MatGetLocalSize(C, &m, &n)); 187 PetscCall(MatGetSize(C, &M, &N)); 188 if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) { 189 PetscCall(MatGetLocalSize(B, NULL, &n)); 190 PetscCall(MatGetSize(B, NULL, &N)); 191 PetscCall(MatGetLocalSize(A, &m, NULL)); 192 PetscCall(MatGetSize(A, &M, NULL)); 193 PetscCall(MatSetSizes(C, m, n, M, N)); 194 } 195 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATMPIDENSE, MATSEQDENSECUDA, MATMPIDENSECUDA, "")); 196 if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)B)->type_name)); 197 PetscCall(MatSetUp(C)); 198 if (!N) { 199 C->ops->productnumeric = MatProductNumeric_Nest_Dense; 200 PetscFunctionReturn(PETSC_SUCCESS); 201 } 202 203 PetscCall(PetscNew(&contents)); 204 C->product->data = contents; 205 C->product->destroy = MatNest_DenseDestroy; 206 PetscCall(PetscCalloc3(nr + 1, &contents->dm, nc + 1, &contents->dn, nr * nc, &contents->workC)); 207 contents->k = nr * nc; 208 for (i = 0; i < nr; i++) { 209 PetscCall(ISGetLocalSize(bA->isglobal.row[i], contents->dm + i + 1)); 210 maxm = PetscMax(maxm, contents->dm[i + 1]); 211 contents->dm[i + 1] += contents->dm[i]; 212 } 213 for (i = 0; i < nc; i++) { 214 PetscCall(ISGetLocalSize(bA->isglobal.col[i], contents->dn + i + 1)); 215 contents->dn[i + 1] += contents->dn[i]; 216 } 217 PetscCall(PetscMalloc1(maxm * N, &contents->tarray)); 218 PetscCall(MatDenseGetLDA(B, &ldb)); 219 PetscCall(MatGetSize(B, NULL, &N)); 220 PetscCall(MatDenseGetArrayRead(B, &barray)); 221 /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */ 222 for (j = 0; j < nc; j++) { 223 PetscCall(ISGetSize(bA->isglobal.col[j], &M)); 224 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), contents->dn[j + 1] - contents->dn[j], PETSC_DECIDE, M, N, (PetscScalar *)(barray + contents->dn[j]), &viewB)); 225 PetscCall(MatDenseSetLDA(viewB, ldb)); 226 for (i = 0; i < nr; i++) { 227 if (!bA->m[i][j]) continue; 228 /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */ 229 230 PetscCall(MatProductCreate(bA->m[i][j], viewB, NULL, &contents->workC[i * nc + j])); 231 workC = contents->workC[i * nc + j]; 232 PetscCall(MatProductSetType(workC, MATPRODUCT_AB)); 233 PetscCall(MatProductSetAlgorithm(workC, "default")); 234 PetscCall(MatProductSetFill(workC, fill)); 235 PetscCall(MatProductSetFromOptions(workC)); 236 PetscCall(MatProductSymbolic(workC)); 237 238 /* since tarray will be shared by all Mat */ 239 PetscCall(MatSeqDenseSetPreallocation(workC, contents->tarray)); 240 PetscCall(MatMPIDenseSetPreallocation(workC, contents->tarray)); 241 } 242 PetscCall(MatDestroy(&viewB)); 243 } 244 PetscCall(MatDenseRestoreArrayRead(B, &barray)); 245 246 C->ops->productnumeric = MatProductNumeric_Nest_Dense; 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 250 static PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C) 251 { 252 Mat_Product *product = C->product; 253 254 PetscFunctionBegin; 255 if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_Nest_Dense; 256 PetscFunctionReturn(PETSC_SUCCESS); 257 } 258 259 static PetscErrorCode MatMultTranspose_Nest(Mat A, Vec x, Vec y) 260 { 261 Mat_Nest *bA = (Mat_Nest *)A->data; 262 Vec *bx = bA->left, *by = bA->right; 263 PetscInt i, j, nr = bA->nr, nc = bA->nc; 264 265 PetscFunctionBegin; 266 for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(x, bA->isglobal.row[i], &bx[i])); 267 for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(y, bA->isglobal.col[i], &by[i])); 268 for (j = 0; j < nc; j++) { 269 PetscCall(VecZeroEntries(by[j])); 270 for (i = 0; i < nr; i++) { 271 if (!bA->m[i][j]) continue; 272 /* y[j] <- y[j] + (A[i][j])^T * x[i] */ 273 PetscCall(MatMultTransposeAdd(bA->m[i][j], bx[i], by[j], by[j])); 274 } 275 } 276 for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i])); 277 for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(y, bA->isglobal.col[i], &by[i])); 278 PetscFunctionReturn(PETSC_SUCCESS); 279 } 280 281 static PetscErrorCode MatMultTransposeAdd_Nest(Mat A, Vec x, Vec y, Vec z) 282 { 283 Mat_Nest *bA = (Mat_Nest *)A->data; 284 Vec *bx = bA->left, *bz = bA->right; 285 PetscInt i, j, nr = bA->nr, nc = bA->nc; 286 287 PetscFunctionBegin; 288 for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(x, bA->isglobal.row[i], &bx[i])); 289 for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(z, bA->isglobal.col[i], &bz[i])); 290 for (j = 0; j < nc; j++) { 291 if (y != z) { 292 Vec by; 293 PetscCall(VecGetSubVector(y, bA->isglobal.col[j], &by)); 294 PetscCall(VecCopy(by, bz[j])); 295 PetscCall(VecRestoreSubVector(y, bA->isglobal.col[j], &by)); 296 } 297 for (i = 0; i < nr; i++) { 298 if (!bA->m[i][j]) continue; 299 /* z[j] <- y[j] + (A[i][j])^T * x[i] */ 300 PetscCall(MatMultTransposeAdd(bA->m[i][j], bx[i], bz[j], bz[j])); 301 } 302 } 303 for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i])); 304 for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(z, bA->isglobal.col[i], &bz[i])); 305 PetscFunctionReturn(PETSC_SUCCESS); 306 } 307 308 static PetscErrorCode MatTranspose_Nest(Mat A, MatReuse reuse, Mat *B) 309 { 310 Mat_Nest *bA = (Mat_Nest *)A->data, *bC; 311 Mat C; 312 PetscInt i, j, nr = bA->nr, nc = bA->nc; 313 314 PetscFunctionBegin; 315 if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B)); 316 PetscCheck(reuse != MAT_INPLACE_MATRIX || nr == nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Square nested matrix only for in-place"); 317 318 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 319 Mat *subs; 320 IS *is_row, *is_col; 321 322 PetscCall(PetscCalloc1(nr * nc, &subs)); 323 PetscCall(PetscMalloc2(nr, &is_row, nc, &is_col)); 324 PetscCall(MatNestGetISs(A, is_row, is_col)); 325 if (reuse == MAT_INPLACE_MATRIX) { 326 for (i = 0; i < nr; i++) { 327 for (j = 0; j < nc; j++) subs[i + nr * j] = bA->m[i][j]; 328 } 329 } 330 331 PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nc, is_col, nr, is_row, subs, &C)); 332 PetscCall(PetscFree(subs)); 333 PetscCall(PetscFree2(is_row, is_col)); 334 } else { 335 C = *B; 336 } 337 338 bC = (Mat_Nest *)C->data; 339 for (i = 0; i < nr; i++) { 340 for (j = 0; j < nc; j++) { 341 if (bA->m[i][j]) { 342 PetscCall(MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]))); 343 } else { 344 bC->m[j][i] = NULL; 345 } 346 } 347 } 348 349 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 350 *B = C; 351 } else { 352 PetscCall(MatHeaderMerge(A, &C)); 353 } 354 PetscFunctionReturn(PETSC_SUCCESS); 355 } 356 357 static PetscErrorCode MatNestDestroyISList(PetscInt n, IS **list) 358 { 359 IS *lst = *list; 360 PetscInt i; 361 362 PetscFunctionBegin; 363 if (!lst) PetscFunctionReturn(PETSC_SUCCESS); 364 for (i = 0; i < n; i++) 365 if (lst[i]) PetscCall(ISDestroy(&lst[i])); 366 PetscCall(PetscFree(lst)); 367 *list = NULL; 368 PetscFunctionReturn(PETSC_SUCCESS); 369 } 370 371 static PetscErrorCode MatReset_Nest(Mat A) 372 { 373 Mat_Nest *vs = (Mat_Nest *)A->data; 374 PetscInt i, j; 375 376 PetscFunctionBegin; 377 /* release the matrices and the place holders */ 378 PetscCall(MatNestDestroyISList(vs->nr, &vs->isglobal.row)); 379 PetscCall(MatNestDestroyISList(vs->nc, &vs->isglobal.col)); 380 PetscCall(MatNestDestroyISList(vs->nr, &vs->islocal.row)); 381 PetscCall(MatNestDestroyISList(vs->nc, &vs->islocal.col)); 382 383 PetscCall(PetscFree(vs->row_len)); 384 PetscCall(PetscFree(vs->col_len)); 385 PetscCall(PetscFree(vs->nnzstate)); 386 387 PetscCall(PetscFree2(vs->left, vs->right)); 388 389 /* release the matrices and the place holders */ 390 if (vs->m) { 391 for (i = 0; i < vs->nr; i++) { 392 for (j = 0; j < vs->nc; j++) PetscCall(MatDestroy(&vs->m[i][j])); 393 PetscCall(PetscFree(vs->m[i])); 394 } 395 PetscCall(PetscFree(vs->m)); 396 } 397 398 /* restore defaults */ 399 vs->nr = 0; 400 vs->nc = 0; 401 vs->splitassembly = PETSC_FALSE; 402 PetscFunctionReturn(PETSC_SUCCESS); 403 } 404 405 static PetscErrorCode MatDestroy_Nest(Mat A) 406 { 407 PetscFunctionBegin; 408 PetscCall(MatReset_Nest(A)); 409 PetscCall(PetscFree(A->data)); 410 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", NULL)); 411 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", NULL)); 412 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", NULL)); 413 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", NULL)); 414 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", NULL)); 415 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", NULL)); 416 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", NULL)); 417 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", NULL)); 418 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", NULL)); 419 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", NULL)); 420 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", NULL)); 421 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", NULL)); 422 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", NULL)); 423 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", NULL)); 424 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", NULL)); 425 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", NULL)); 426 PetscFunctionReturn(PETSC_SUCCESS); 427 } 428 429 static PetscErrorCode MatMissingDiagonal_Nest(Mat mat, PetscBool *missing, PetscInt *dd) 430 { 431 Mat_Nest *vs = (Mat_Nest *)mat->data; 432 PetscInt i; 433 434 PetscFunctionBegin; 435 if (dd) *dd = 0; 436 if (!vs->nr) { 437 *missing = PETSC_TRUE; 438 PetscFunctionReturn(PETSC_SUCCESS); 439 } 440 *missing = PETSC_FALSE; 441 for (i = 0; i < vs->nr && !(*missing); i++) { 442 *missing = PETSC_TRUE; 443 if (vs->m[i][i]) { 444 PetscCall(MatMissingDiagonal(vs->m[i][i], missing, NULL)); 445 PetscCheck(!*missing || !dd, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "First missing entry not yet implemented"); 446 } 447 } 448 PetscFunctionReturn(PETSC_SUCCESS); 449 } 450 451 static PetscErrorCode MatAssemblyBegin_Nest(Mat A, MatAssemblyType type) 452 { 453 Mat_Nest *vs = (Mat_Nest *)A->data; 454 PetscInt i, j; 455 PetscBool nnzstate = PETSC_FALSE; 456 457 PetscFunctionBegin; 458 for (i = 0; i < vs->nr; i++) { 459 for (j = 0; j < vs->nc; j++) { 460 PetscObjectState subnnzstate = 0; 461 if (vs->m[i][j]) { 462 PetscCall(MatAssemblyBegin(vs->m[i][j], type)); 463 if (!vs->splitassembly) { 464 /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested 465 * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was 466 * already performing an assembly, but the result would by more complicated and appears to offer less 467 * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an 468 * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives". 469 */ 470 PetscCall(MatAssemblyEnd(vs->m[i][j], type)); 471 PetscCall(MatGetNonzeroState(vs->m[i][j], &subnnzstate)); 472 } 473 } 474 nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i * vs->nc + j] != subnnzstate); 475 vs->nnzstate[i * vs->nc + j] = subnnzstate; 476 } 477 } 478 if (nnzstate) A->nonzerostate++; 479 PetscFunctionReturn(PETSC_SUCCESS); 480 } 481 482 static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type) 483 { 484 Mat_Nest *vs = (Mat_Nest *)A->data; 485 PetscInt i, j; 486 487 PetscFunctionBegin; 488 for (i = 0; i < vs->nr; i++) { 489 for (j = 0; j < vs->nc; j++) { 490 if (vs->m[i][j]) { 491 if (vs->splitassembly) PetscCall(MatAssemblyEnd(vs->m[i][j], type)); 492 } 493 } 494 } 495 PetscFunctionReturn(PETSC_SUCCESS); 496 } 497 498 static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A, PetscInt row, Mat *B) 499 { 500 Mat_Nest *vs = (Mat_Nest *)A->data; 501 PetscInt j; 502 Mat sub; 503 504 PetscFunctionBegin; 505 sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */ 506 for (j = 0; !sub && j < vs->nc; j++) sub = vs->m[row][j]; 507 if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */ 508 *B = sub; 509 PetscFunctionReturn(PETSC_SUCCESS); 510 } 511 512 static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A, PetscInt col, Mat *B) 513 { 514 Mat_Nest *vs = (Mat_Nest *)A->data; 515 PetscInt i; 516 Mat sub; 517 518 PetscFunctionBegin; 519 sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */ 520 for (i = 0; !sub && i < vs->nr; i++) sub = vs->m[i][col]; 521 if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */ 522 *B = sub; 523 PetscFunctionReturn(PETSC_SUCCESS); 524 } 525 526 static PetscErrorCode MatNestFindISRange(Mat A, PetscInt n, const IS list[], IS is, PetscInt *begin, PetscInt *end) 527 { 528 PetscInt i, j, size, m; 529 PetscBool flg; 530 IS out, concatenate[2]; 531 532 PetscFunctionBegin; 533 PetscAssertPointer(list, 3); 534 PetscValidHeaderSpecific(is, IS_CLASSID, 4); 535 if (begin) { 536 PetscAssertPointer(begin, 5); 537 *begin = -1; 538 } 539 if (end) { 540 PetscAssertPointer(end, 6); 541 *end = -1; 542 } 543 for (i = 0; i < n; i++) { 544 if (!list[i]) continue; 545 PetscCall(ISEqualUnsorted(list[i], is, &flg)); 546 if (flg) { 547 if (begin) *begin = i; 548 if (end) *end = i + 1; 549 PetscFunctionReturn(PETSC_SUCCESS); 550 } 551 } 552 PetscCall(ISGetSize(is, &size)); 553 for (i = 0; i < n - 1; i++) { 554 if (!list[i]) continue; 555 m = 0; 556 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), 2, list + i, &out)); 557 PetscCall(ISGetSize(out, &m)); 558 for (j = i + 2; j < n && m < size; j++) { 559 if (list[j]) { 560 concatenate[0] = out; 561 concatenate[1] = list[j]; 562 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), 2, concatenate, &out)); 563 PetscCall(ISDestroy(concatenate)); 564 PetscCall(ISGetSize(out, &m)); 565 } 566 } 567 if (m == size) { 568 PetscCall(ISEqualUnsorted(out, is, &flg)); 569 if (flg) { 570 if (begin) *begin = i; 571 if (end) *end = j; 572 PetscCall(ISDestroy(&out)); 573 PetscFunctionReturn(PETSC_SUCCESS); 574 } 575 } 576 PetscCall(ISDestroy(&out)); 577 } 578 PetscFunctionReturn(PETSC_SUCCESS); 579 } 580 581 static PetscErrorCode MatNestFillEmptyMat_Private(Mat A, PetscInt i, PetscInt j, Mat *B) 582 { 583 Mat_Nest *vs = (Mat_Nest *)A->data; 584 PetscInt lr, lc; 585 586 PetscFunctionBegin; 587 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B)); 588 PetscCall(ISGetLocalSize(vs->isglobal.row[i], &lr)); 589 PetscCall(ISGetLocalSize(vs->isglobal.col[j], &lc)); 590 PetscCall(MatSetSizes(*B, lr, lc, PETSC_DECIDE, PETSC_DECIDE)); 591 PetscCall(MatSetType(*B, MATAIJ)); 592 PetscCall(MatSeqAIJSetPreallocation(*B, 0, NULL)); 593 PetscCall(MatMPIAIJSetPreallocation(*B, 0, NULL, 0, NULL)); 594 PetscCall(MatSetUp(*B)); 595 PetscCall(MatSetOption(*B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 596 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 597 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 598 PetscFunctionReturn(PETSC_SUCCESS); 599 } 600 601 static PetscErrorCode MatNestGetBlock_Private(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *B) 602 { 603 Mat_Nest *vs = (Mat_Nest *)A->data; 604 Mat *a; 605 PetscInt i, j, k, l, nr = rend - rbegin, nc = cend - cbegin; 606 char keyname[256]; 607 PetscBool *b; 608 PetscBool flg; 609 610 PetscFunctionBegin; 611 *B = NULL; 612 PetscCall(PetscSNPrintf(keyname, sizeof(keyname), "NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT, rbegin, rend, cbegin, cend)); 613 PetscCall(PetscObjectQuery((PetscObject)A, keyname, (PetscObject *)B)); 614 if (*B) PetscFunctionReturn(PETSC_SUCCESS); 615 616 PetscCall(PetscMalloc2(nr * nc, &a, nr * nc, &b)); 617 for (i = 0; i < nr; i++) { 618 for (j = 0; j < nc; j++) { 619 a[i * nc + j] = vs->m[rbegin + i][cbegin + j]; 620 b[i * nc + j] = PETSC_FALSE; 621 } 622 } 623 if (nc != vs->nc && nr != vs->nr) { 624 for (i = 0; i < nr; i++) { 625 for (j = 0; j < nc; j++) { 626 flg = PETSC_FALSE; 627 for (k = 0; (k < nr && !flg); k++) { 628 if (a[j + k * nc]) flg = PETSC_TRUE; 629 } 630 if (flg) { 631 flg = PETSC_FALSE; 632 for (l = 0; (l < nc && !flg); l++) { 633 if (a[i * nc + l]) flg = PETSC_TRUE; 634 } 635 } 636 if (!flg) { 637 b[i * nc + j] = PETSC_TRUE; 638 PetscCall(MatNestFillEmptyMat_Private(A, rbegin + i, cbegin + j, a + i * nc + j)); 639 } 640 } 641 } 642 } 643 PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nr, nr != vs->nr ? NULL : vs->isglobal.row, nc, nc != vs->nc ? NULL : vs->isglobal.col, a, B)); 644 for (i = 0; i < nr; i++) { 645 for (j = 0; j < nc; j++) { 646 if (b[i * nc + j]) PetscCall(MatDestroy(a + i * nc + j)); 647 } 648 } 649 PetscCall(PetscFree2(a, b)); 650 (*B)->assembled = A->assembled; 651 PetscCall(PetscObjectCompose((PetscObject)A, keyname, (PetscObject)*B)); 652 PetscCall(PetscObjectDereference((PetscObject)*B)); /* Leave the only remaining reference in the composition */ 653 PetscFunctionReturn(PETSC_SUCCESS); 654 } 655 656 static PetscErrorCode MatNestFindSubMat(Mat A, struct MatNestISPair *is, IS isrow, IS iscol, Mat *B) 657 { 658 Mat_Nest *vs = (Mat_Nest *)A->data; 659 PetscInt rbegin, rend, cbegin, cend; 660 661 PetscFunctionBegin; 662 PetscCall(MatNestFindISRange(A, vs->nr, is->row, isrow, &rbegin, &rend)); 663 PetscCall(MatNestFindISRange(A, vs->nc, is->col, iscol, &cbegin, &cend)); 664 if (rend == rbegin + 1 && cend == cbegin + 1) { 665 if (!vs->m[rbegin][cbegin]) PetscCall(MatNestFillEmptyMat_Private(A, rbegin, cbegin, vs->m[rbegin] + cbegin)); 666 *B = vs->m[rbegin][cbegin]; 667 } else if (rbegin != -1 && cbegin != -1) { 668 PetscCall(MatNestGetBlock_Private(A, rbegin, rend, cbegin, cend, B)); 669 } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Could not find index set"); 670 PetscFunctionReturn(PETSC_SUCCESS); 671 } 672 673 /* 674 TODO: This does not actually returns a submatrix we can modify 675 */ 676 static PetscErrorCode MatCreateSubMatrix_Nest(Mat A, IS isrow, IS iscol, MatReuse reuse, Mat *B) 677 { 678 Mat_Nest *vs = (Mat_Nest *)A->data; 679 Mat sub; 680 681 PetscFunctionBegin; 682 PetscCall(MatNestFindSubMat(A, &vs->isglobal, isrow, iscol, &sub)); 683 switch (reuse) { 684 case MAT_INITIAL_MATRIX: 685 if (sub) PetscCall(PetscObjectReference((PetscObject)sub)); 686 *B = sub; 687 break; 688 case MAT_REUSE_MATRIX: 689 PetscCheck(sub == *B, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Submatrix was not used before in this call"); 690 break; 691 case MAT_IGNORE_MATRIX: /* Nothing to do */ 692 break; 693 case MAT_INPLACE_MATRIX: /* Nothing to do */ 694 SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_INPLACE_MATRIX is not supported yet"); 695 } 696 PetscFunctionReturn(PETSC_SUCCESS); 697 } 698 699 static PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B) 700 { 701 Mat_Nest *vs = (Mat_Nest *)A->data; 702 Mat sub; 703 704 PetscFunctionBegin; 705 PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub)); 706 /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */ 707 if (sub) PetscCall(PetscObjectReference((PetscObject)sub)); 708 *B = sub; 709 PetscFunctionReturn(PETSC_SUCCESS); 710 } 711 712 static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B) 713 { 714 Mat_Nest *vs = (Mat_Nest *)A->data; 715 Mat sub; 716 717 PetscFunctionBegin; 718 PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub)); 719 PetscCheck(*B == sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has not been gotten"); 720 if (sub) { 721 PetscCheck(((PetscObject)sub)->refct > 1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has had reference count decremented too many times"); 722 PetscCall(MatDestroy(B)); 723 } 724 PetscFunctionReturn(PETSC_SUCCESS); 725 } 726 727 static PetscErrorCode MatGetDiagonal_Nest(Mat A, Vec v) 728 { 729 Mat_Nest *bA = (Mat_Nest *)A->data; 730 PetscInt i; 731 732 PetscFunctionBegin; 733 for (i = 0; i < bA->nr; i++) { 734 Vec bv; 735 PetscCall(VecGetSubVector(v, bA->isglobal.row[i], &bv)); 736 if (bA->m[i][i]) { 737 PetscCall(MatGetDiagonal(bA->m[i][i], bv)); 738 } else { 739 PetscCall(VecSet(bv, 0.0)); 740 } 741 PetscCall(VecRestoreSubVector(v, bA->isglobal.row[i], &bv)); 742 } 743 PetscFunctionReturn(PETSC_SUCCESS); 744 } 745 746 static PetscErrorCode MatDiagonalScale_Nest(Mat A, Vec l, Vec r) 747 { 748 Mat_Nest *bA = (Mat_Nest *)A->data; 749 Vec bl, *br; 750 PetscInt i, j; 751 752 PetscFunctionBegin; 753 PetscCall(PetscCalloc1(bA->nc, &br)); 754 if (r) { 755 for (j = 0; j < bA->nc; j++) PetscCall(VecGetSubVector(r, bA->isglobal.col[j], &br[j])); 756 } 757 bl = NULL; 758 for (i = 0; i < bA->nr; i++) { 759 if (l) PetscCall(VecGetSubVector(l, bA->isglobal.row[i], &bl)); 760 for (j = 0; j < bA->nc; j++) { 761 if (bA->m[i][j]) PetscCall(MatDiagonalScale(bA->m[i][j], bl, br[j])); 762 } 763 if (l) PetscCall(VecRestoreSubVector(l, bA->isglobal.row[i], &bl)); 764 } 765 if (r) { 766 for (j = 0; j < bA->nc; j++) PetscCall(VecRestoreSubVector(r, bA->isglobal.col[j], &br[j])); 767 } 768 PetscCall(PetscFree(br)); 769 PetscFunctionReturn(PETSC_SUCCESS); 770 } 771 772 static PetscErrorCode MatScale_Nest(Mat A, PetscScalar a) 773 { 774 Mat_Nest *bA = (Mat_Nest *)A->data; 775 PetscInt i, j; 776 777 PetscFunctionBegin; 778 for (i = 0; i < bA->nr; i++) { 779 for (j = 0; j < bA->nc; j++) { 780 if (bA->m[i][j]) PetscCall(MatScale(bA->m[i][j], a)); 781 } 782 } 783 PetscFunctionReturn(PETSC_SUCCESS); 784 } 785 786 static PetscErrorCode MatShift_Nest(Mat A, PetscScalar a) 787 { 788 Mat_Nest *bA = (Mat_Nest *)A->data; 789 PetscInt i; 790 PetscBool nnzstate = PETSC_FALSE; 791 792 PetscFunctionBegin; 793 for (i = 0; i < bA->nr; i++) { 794 PetscObjectState subnnzstate = 0; 795 PetscCheck(bA->m[i][i], PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for shifting an empty diagonal block, insert a matrix in block (%" PetscInt_FMT ",%" PetscInt_FMT ")", i, i); 796 PetscCall(MatShift(bA->m[i][i], a)); 797 PetscCall(MatGetNonzeroState(bA->m[i][i], &subnnzstate)); 798 nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate); 799 bA->nnzstate[i * bA->nc + i] = subnnzstate; 800 } 801 if (nnzstate) A->nonzerostate++; 802 PetscFunctionReturn(PETSC_SUCCESS); 803 } 804 805 static PetscErrorCode MatDiagonalSet_Nest(Mat A, Vec D, InsertMode is) 806 { 807 Mat_Nest *bA = (Mat_Nest *)A->data; 808 PetscInt i; 809 PetscBool nnzstate = PETSC_FALSE; 810 811 PetscFunctionBegin; 812 for (i = 0; i < bA->nr; i++) { 813 PetscObjectState subnnzstate = 0; 814 Vec bv; 815 PetscCall(VecGetSubVector(D, bA->isglobal.row[i], &bv)); 816 if (bA->m[i][i]) { 817 PetscCall(MatDiagonalSet(bA->m[i][i], bv, is)); 818 PetscCall(MatGetNonzeroState(bA->m[i][i], &subnnzstate)); 819 } 820 PetscCall(VecRestoreSubVector(D, bA->isglobal.row[i], &bv)); 821 nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate); 822 bA->nnzstate[i * bA->nc + i] = subnnzstate; 823 } 824 if (nnzstate) A->nonzerostate++; 825 PetscFunctionReturn(PETSC_SUCCESS); 826 } 827 828 static PetscErrorCode MatSetRandom_Nest(Mat A, PetscRandom rctx) 829 { 830 Mat_Nest *bA = (Mat_Nest *)A->data; 831 PetscInt i, j; 832 833 PetscFunctionBegin; 834 for (i = 0; i < bA->nr; i++) { 835 for (j = 0; j < bA->nc; j++) { 836 if (bA->m[i][j]) PetscCall(MatSetRandom(bA->m[i][j], rctx)); 837 } 838 } 839 PetscFunctionReturn(PETSC_SUCCESS); 840 } 841 842 static PetscErrorCode MatCreateVecs_Nest(Mat A, Vec *right, Vec *left) 843 { 844 Mat_Nest *bA = (Mat_Nest *)A->data; 845 Vec *L, *R; 846 MPI_Comm comm; 847 PetscInt i, j; 848 849 PetscFunctionBegin; 850 PetscCall(PetscObjectGetComm((PetscObject)A, &comm)); 851 if (right) { 852 /* allocate R */ 853 PetscCall(PetscMalloc1(bA->nc, &R)); 854 /* Create the right vectors */ 855 for (j = 0; j < bA->nc; j++) { 856 for (i = 0; i < bA->nr; i++) { 857 if (bA->m[i][j]) { 858 PetscCall(MatCreateVecs(bA->m[i][j], &R[j], NULL)); 859 break; 860 } 861 } 862 PetscCheck(i != bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column."); 863 } 864 PetscCall(VecCreateNest(comm, bA->nc, bA->isglobal.col, R, right)); 865 /* hand back control to the nest vector */ 866 for (j = 0; j < bA->nc; j++) PetscCall(VecDestroy(&R[j])); 867 PetscCall(PetscFree(R)); 868 } 869 870 if (left) { 871 /* allocate L */ 872 PetscCall(PetscMalloc1(bA->nr, &L)); 873 /* Create the left vectors */ 874 for (i = 0; i < bA->nr; i++) { 875 for (j = 0; j < bA->nc; j++) { 876 if (bA->m[i][j]) { 877 PetscCall(MatCreateVecs(bA->m[i][j], NULL, &L[i])); 878 break; 879 } 880 } 881 PetscCheck(j != bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row."); 882 } 883 884 PetscCall(VecCreateNest(comm, bA->nr, bA->isglobal.row, L, left)); 885 for (i = 0; i < bA->nr; i++) PetscCall(VecDestroy(&L[i])); 886 887 PetscCall(PetscFree(L)); 888 } 889 PetscFunctionReturn(PETSC_SUCCESS); 890 } 891 892 static PetscErrorCode MatView_Nest(Mat A, PetscViewer viewer) 893 { 894 Mat_Nest *bA = (Mat_Nest *)A->data; 895 PetscBool isascii, viewSub = PETSC_FALSE; 896 PetscInt i, j; 897 898 PetscFunctionBegin; 899 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 900 if (isascii) { 901 PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_view_nest_sub", &viewSub, NULL)); 902 PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix object:\n")); 903 PetscCall(PetscViewerASCIIPushTab(viewer)); 904 PetscCall(PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", bA->nr, bA->nc)); 905 906 PetscCall(PetscViewerASCIIPrintf(viewer, "MatNest structure:\n")); 907 for (i = 0; i < bA->nr; i++) { 908 for (j = 0; j < bA->nc; j++) { 909 MatType type; 910 char name[256] = "", prefix[256] = ""; 911 PetscInt NR, NC; 912 PetscBool isNest = PETSC_FALSE; 913 914 if (!bA->m[i][j]) { 915 PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL\n", i, j)); 916 continue; 917 } 918 PetscCall(MatGetSize(bA->m[i][j], &NR, &NC)); 919 PetscCall(MatGetType(bA->m[i][j], &type)); 920 if (((PetscObject)bA->m[i][j])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bA->m[i][j])->name)); 921 if (((PetscObject)bA->m[i][j])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bA->m[i][j])->prefix)); 922 PetscCall(PetscObjectTypeCompare((PetscObject)bA->m[i][j], MATNEST, &isNest)); 923 924 PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", i, j, name, prefix, type, NR, NC)); 925 926 if (isNest || viewSub) { 927 PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */ 928 PetscCall(MatView(bA->m[i][j], viewer)); 929 PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */ 930 } 931 } 932 } 933 PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */ 934 } 935 PetscFunctionReturn(PETSC_SUCCESS); 936 } 937 938 static PetscErrorCode MatZeroEntries_Nest(Mat A) 939 { 940 Mat_Nest *bA = (Mat_Nest *)A->data; 941 PetscInt i, j; 942 943 PetscFunctionBegin; 944 for (i = 0; i < bA->nr; i++) { 945 for (j = 0; j < bA->nc; j++) { 946 if (!bA->m[i][j]) continue; 947 PetscCall(MatZeroEntries(bA->m[i][j])); 948 } 949 } 950 PetscFunctionReturn(PETSC_SUCCESS); 951 } 952 953 static PetscErrorCode MatCopy_Nest(Mat A, Mat B, MatStructure str) 954 { 955 Mat_Nest *bA = (Mat_Nest *)A->data, *bB = (Mat_Nest *)B->data; 956 PetscInt i, j, nr = bA->nr, nc = bA->nc; 957 PetscBool nnzstate = PETSC_FALSE; 958 959 PetscFunctionBegin; 960 PetscCheck(nr == bB->nr && nc == bB->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Cannot copy a Mat_Nest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ") to a Mat_Nest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ")", bB->nr, bB->nc, nr, nc); 961 for (i = 0; i < nr; i++) { 962 for (j = 0; j < nc; j++) { 963 PetscObjectState subnnzstate = 0; 964 if (bA->m[i][j] && bB->m[i][j]) { 965 PetscCall(MatCopy(bA->m[i][j], bB->m[i][j], str)); 966 } else PetscCheck(!bA->m[i][j] && !bB->m[i][j], PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT, i, j); 967 PetscCall(MatGetNonzeroState(bB->m[i][j], &subnnzstate)); 968 nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i * nc + j] != subnnzstate); 969 bB->nnzstate[i * nc + j] = subnnzstate; 970 } 971 } 972 if (nnzstate) B->nonzerostate++; 973 PetscFunctionReturn(PETSC_SUCCESS); 974 } 975 976 static PetscErrorCode MatAXPY_Nest(Mat Y, PetscScalar a, Mat X, MatStructure str) 977 { 978 Mat_Nest *bY = (Mat_Nest *)Y->data, *bX = (Mat_Nest *)X->data; 979 PetscInt i, j, nr = bY->nr, nc = bY->nc; 980 PetscBool nnzstate = PETSC_FALSE; 981 982 PetscFunctionBegin; 983 PetscCheck(nr == bX->nr && nc == bX->nc, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_INCOMP, "Cannot AXPY a MatNest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ") with a MatNest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ")", bX->nr, bX->nc, nr, nc); 984 for (i = 0; i < nr; i++) { 985 for (j = 0; j < nc; j++) { 986 PetscObjectState subnnzstate = 0; 987 if (bY->m[i][j] && bX->m[i][j]) { 988 PetscCall(MatAXPY(bY->m[i][j], a, bX->m[i][j], str)); 989 } else if (bX->m[i][j]) { 990 Mat M; 991 992 PetscCheck(str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT ". Use DIFFERENT_NONZERO_PATTERN or UNKNOWN_NONZERO_PATTERN", i, j); 993 PetscCall(MatDuplicate(bX->m[i][j], MAT_COPY_VALUES, &M)); 994 PetscCall(MatNestSetSubMat(Y, i, j, M)); 995 PetscCall(MatDestroy(&M)); 996 } 997 if (bY->m[i][j]) PetscCall(MatGetNonzeroState(bY->m[i][j], &subnnzstate)); 998 nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i * nc + j] != subnnzstate); 999 bY->nnzstate[i * nc + j] = subnnzstate; 1000 } 1001 } 1002 if (nnzstate) Y->nonzerostate++; 1003 PetscFunctionReturn(PETSC_SUCCESS); 1004 } 1005 1006 static PetscErrorCode MatDuplicate_Nest(Mat A, MatDuplicateOption op, Mat *B) 1007 { 1008 Mat_Nest *bA = (Mat_Nest *)A->data; 1009 Mat *b; 1010 PetscInt i, j, nr = bA->nr, nc = bA->nc; 1011 1012 PetscFunctionBegin; 1013 PetscCall(PetscMalloc1(nr * nc, &b)); 1014 for (i = 0; i < nr; i++) { 1015 for (j = 0; j < nc; j++) { 1016 if (bA->m[i][j]) { 1017 PetscCall(MatDuplicate(bA->m[i][j], op, &b[i * nc + j])); 1018 } else { 1019 b[i * nc + j] = NULL; 1020 } 1021 } 1022 } 1023 PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nr, bA->isglobal.row, nc, bA->isglobal.col, b, B)); 1024 /* Give the new MatNest exclusive ownership */ 1025 for (i = 0; i < nr * nc; i++) PetscCall(MatDestroy(&b[i])); 1026 PetscCall(PetscFree(b)); 1027 1028 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 1029 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 1030 PetscFunctionReturn(PETSC_SUCCESS); 1031 } 1032 1033 /* nest api */ 1034 static PetscErrorCode MatNestGetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat *mat) 1035 { 1036 Mat_Nest *bA = (Mat_Nest *)A->data; 1037 1038 PetscFunctionBegin; 1039 PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1); 1040 PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1); 1041 *mat = bA->m[idxm][jdxm]; 1042 PetscFunctionReturn(PETSC_SUCCESS); 1043 } 1044 1045 /*@ 1046 MatNestGetSubMat - Returns a single, sub-matrix from a `MATNEST` 1047 1048 Not Collective 1049 1050 Input Parameters: 1051 + A - `MATNEST` matrix 1052 . idxm - index of the matrix within the nest matrix 1053 - jdxm - index of the matrix within the nest matrix 1054 1055 Output Parameter: 1056 . sub - matrix at index `idxm`, `jdxm` within the nest matrix 1057 1058 Level: developer 1059 1060 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestSetSubMat()`, 1061 `MatNestGetLocalISs()`, `MatNestGetISs()` 1062 @*/ 1063 PetscErrorCode MatNestGetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat *sub) 1064 { 1065 PetscFunctionBegin; 1066 PetscUseMethod(A, "MatNestGetSubMat_C", (Mat, PetscInt, PetscInt, Mat *), (A, idxm, jdxm, sub)); 1067 PetscFunctionReturn(PETSC_SUCCESS); 1068 } 1069 1070 static PetscErrorCode MatNestSetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat mat) 1071 { 1072 Mat_Nest *bA = (Mat_Nest *)A->data; 1073 PetscInt m, n, M, N, mi, ni, Mi, Ni; 1074 1075 PetscFunctionBegin; 1076 PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1); 1077 PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1); 1078 PetscCall(MatGetLocalSize(mat, &m, &n)); 1079 PetscCall(MatGetSize(mat, &M, &N)); 1080 PetscCall(ISGetLocalSize(bA->isglobal.row[idxm], &mi)); 1081 PetscCall(ISGetSize(bA->isglobal.row[idxm], &Mi)); 1082 PetscCall(ISGetLocalSize(bA->isglobal.col[jdxm], &ni)); 1083 PetscCall(ISGetSize(bA->isglobal.col[jdxm], &Ni)); 1084 PetscCheck(M == Mi && N == Ni, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_INCOMP, "Submatrix dimension (%" PetscInt_FMT ",%" PetscInt_FMT ") incompatible with nest block (%" PetscInt_FMT ",%" PetscInt_FMT ")", M, N, Mi, Ni); 1085 PetscCheck(m == mi && n == ni, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_INCOMP, "Submatrix local dimension (%" PetscInt_FMT ",%" PetscInt_FMT ") incompatible with nest block (%" PetscInt_FMT ",%" PetscInt_FMT ")", m, n, mi, ni); 1086 1087 /* do not increase object state */ 1088 if (mat == bA->m[idxm][jdxm]) PetscFunctionReturn(PETSC_SUCCESS); 1089 1090 PetscCall(PetscObjectReference((PetscObject)mat)); 1091 PetscCall(MatDestroy(&bA->m[idxm][jdxm])); 1092 bA->m[idxm][jdxm] = mat; 1093 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 1094 PetscCall(MatGetNonzeroState(mat, &bA->nnzstate[idxm * bA->nc + jdxm])); 1095 A->nonzerostate++; 1096 PetscFunctionReturn(PETSC_SUCCESS); 1097 } 1098 1099 /*@ 1100 MatNestSetSubMat - Set a single submatrix in the `MATNEST` 1101 1102 Logically Collective 1103 1104 Input Parameters: 1105 + A - `MATNEST` matrix 1106 . idxm - index of the matrix within the nest matrix 1107 . jdxm - index of the matrix within the nest matrix 1108 - sub - matrix at index `idxm`, `jdxm` within the nest matrix 1109 1110 Level: developer 1111 1112 Notes: 1113 The new submatrix must have the same size and communicator as that block of the nest. 1114 1115 This increments the reference count of the submatrix. 1116 1117 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestSetSubMats()`, `MatNestGetSubMats()`, `MatNestGetLocalISs()`, `MatCreateNest()`, 1118 `MatNestGetSubMat()`, `MatNestGetISs()`, `MatNestGetSize()` 1119 @*/ 1120 PetscErrorCode MatNestSetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat sub) 1121 { 1122 PetscFunctionBegin; 1123 PetscUseMethod(A, "MatNestSetSubMat_C", (Mat, PetscInt, PetscInt, Mat), (A, idxm, jdxm, sub)); 1124 PetscFunctionReturn(PETSC_SUCCESS); 1125 } 1126 1127 static PetscErrorCode MatNestGetSubMats_Nest(Mat A, PetscInt *M, PetscInt *N, Mat ***mat) 1128 { 1129 Mat_Nest *bA = (Mat_Nest *)A->data; 1130 1131 PetscFunctionBegin; 1132 if (M) *M = bA->nr; 1133 if (N) *N = bA->nc; 1134 if (mat) *mat = bA->m; 1135 PetscFunctionReturn(PETSC_SUCCESS); 1136 } 1137 1138 /*@C 1139 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a `MATNEST` matrix. 1140 1141 Not Collective 1142 1143 Input Parameter: 1144 . A - nest matrix 1145 1146 Output Parameters: 1147 + M - number of rows in the nest matrix 1148 . N - number of cols in the nest matrix 1149 - mat - array of matrices 1150 1151 Level: developer 1152 1153 Note: 1154 The user should not free the array `mat`. 1155 1156 Fortran Notes: 1157 This routine has a calling sequence 1158 $ call MatNestGetSubMats(A, M, N, mat, ierr) 1159 where the space allocated for the optional argument `mat` is assumed large enough (if provided). 1160 Matrices in `mat` are returned in row-major order, see `MatCreateNest()` for an example. 1161 1162 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatCreateNest()`, 1163 `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()` 1164 @*/ 1165 PetscErrorCode MatNestGetSubMats(Mat A, PetscInt *M, PetscInt *N, Mat ***mat) 1166 { 1167 PetscFunctionBegin; 1168 PetscUseMethod(A, "MatNestGetSubMats_C", (Mat, PetscInt *, PetscInt *, Mat ***), (A, M, N, mat)); 1169 PetscFunctionReturn(PETSC_SUCCESS); 1170 } 1171 1172 static PetscErrorCode MatNestGetSize_Nest(Mat A, PetscInt *M, PetscInt *N) 1173 { 1174 Mat_Nest *bA = (Mat_Nest *)A->data; 1175 1176 PetscFunctionBegin; 1177 if (M) *M = bA->nr; 1178 if (N) *N = bA->nc; 1179 PetscFunctionReturn(PETSC_SUCCESS); 1180 } 1181 1182 /*@ 1183 MatNestGetSize - Returns the size of the `MATNEST` matrix. 1184 1185 Not Collective 1186 1187 Input Parameter: 1188 . A - `MATNEST` matrix 1189 1190 Output Parameters: 1191 + M - number of rows in the nested mat 1192 - N - number of cols in the nested mat 1193 1194 Level: developer 1195 1196 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestGetLocalISs()`, 1197 `MatNestGetISs()` 1198 @*/ 1199 PetscErrorCode MatNestGetSize(Mat A, PetscInt *M, PetscInt *N) 1200 { 1201 PetscFunctionBegin; 1202 PetscUseMethod(A, "MatNestGetSize_C", (Mat, PetscInt *, PetscInt *), (A, M, N)); 1203 PetscFunctionReturn(PETSC_SUCCESS); 1204 } 1205 1206 static PetscErrorCode MatNestGetISs_Nest(Mat A, IS rows[], IS cols[]) 1207 { 1208 Mat_Nest *vs = (Mat_Nest *)A->data; 1209 PetscInt i; 1210 1211 PetscFunctionBegin; 1212 if (rows) 1213 for (i = 0; i < vs->nr; i++) rows[i] = vs->isglobal.row[i]; 1214 if (cols) 1215 for (i = 0; i < vs->nc; i++) cols[i] = vs->isglobal.col[i]; 1216 PetscFunctionReturn(PETSC_SUCCESS); 1217 } 1218 1219 /*@C 1220 MatNestGetISs - Returns the index sets partitioning the row and column spaces of a `MATNEST` 1221 1222 Not Collective 1223 1224 Input Parameter: 1225 . A - `MATNEST` matrix 1226 1227 Output Parameters: 1228 + rows - array of row index sets 1229 - cols - array of column index sets 1230 1231 Level: advanced 1232 1233 Note: 1234 The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s. 1235 1236 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`, 1237 `MatCreateNest()`, `MatNestSetSubMats()` 1238 @*/ 1239 PetscErrorCode MatNestGetISs(Mat A, IS rows[], IS cols[]) 1240 { 1241 PetscFunctionBegin; 1242 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1243 PetscUseMethod(A, "MatNestGetISs_C", (Mat, IS[], IS[]), (A, rows, cols)); 1244 PetscFunctionReturn(PETSC_SUCCESS); 1245 } 1246 1247 static PetscErrorCode MatNestGetLocalISs_Nest(Mat A, IS rows[], IS cols[]) 1248 { 1249 Mat_Nest *vs = (Mat_Nest *)A->data; 1250 PetscInt i; 1251 1252 PetscFunctionBegin; 1253 if (rows) 1254 for (i = 0; i < vs->nr; i++) rows[i] = vs->islocal.row[i]; 1255 if (cols) 1256 for (i = 0; i < vs->nc; i++) cols[i] = vs->islocal.col[i]; 1257 PetscFunctionReturn(PETSC_SUCCESS); 1258 } 1259 1260 /*@C 1261 MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces of a `MATNEST` 1262 1263 Not Collective 1264 1265 Input Parameter: 1266 . A - `MATNEST` matrix 1267 1268 Output Parameters: 1269 + rows - array of row index sets (or `NULL` to ignore) 1270 - cols - array of column index sets (or `NULL` to ignore) 1271 1272 Level: advanced 1273 1274 Note: 1275 The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s. 1276 1277 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`, 1278 `MatNestSetSubMats()`, `MatNestSetSubMat()` 1279 @*/ 1280 PetscErrorCode MatNestGetLocalISs(Mat A, IS rows[], IS cols[]) 1281 { 1282 PetscFunctionBegin; 1283 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1284 PetscUseMethod(A, "MatNestGetLocalISs_C", (Mat, IS[], IS[]), (A, rows, cols)); 1285 PetscFunctionReturn(PETSC_SUCCESS); 1286 } 1287 1288 static PetscErrorCode MatNestSetVecType_Nest(Mat A, VecType vtype) 1289 { 1290 PetscBool flg; 1291 1292 PetscFunctionBegin; 1293 PetscCall(PetscStrcmp(vtype, VECNEST, &flg)); 1294 /* In reality, this only distinguishes VECNEST and "other" */ 1295 if (flg) A->ops->getvecs = MatCreateVecs_Nest; 1296 else A->ops->getvecs = (PetscErrorCode(*)(Mat, Vec *, Vec *))0; 1297 PetscFunctionReturn(PETSC_SUCCESS); 1298 } 1299 1300 /*@C 1301 MatNestSetVecType - Sets the type of `Vec` returned by `MatCreateVecs()` 1302 1303 Not Collective 1304 1305 Input Parameters: 1306 + A - `MATNEST` matrix 1307 - vtype - `VecType` to use for creating vectors 1308 1309 Level: developer 1310 1311 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateVecs()`, `MatCreateNest()`, `VecType` 1312 @*/ 1313 PetscErrorCode MatNestSetVecType(Mat A, VecType vtype) 1314 { 1315 PetscFunctionBegin; 1316 PetscTryMethod(A, "MatNestSetVecType_C", (Mat, VecType), (A, vtype)); 1317 PetscFunctionReturn(PETSC_SUCCESS); 1318 } 1319 1320 static PetscErrorCode MatNestSetSubMats_Nest(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[]) 1321 { 1322 Mat_Nest *s = (Mat_Nest *)A->data; 1323 PetscInt i, j, m, n, M, N; 1324 PetscBool cong, isstd, sametype = PETSC_FALSE; 1325 VecType vtype, type; 1326 1327 PetscFunctionBegin; 1328 PetscCall(MatReset_Nest(A)); 1329 1330 s->nr = nr; 1331 s->nc = nc; 1332 1333 /* Create space for submatrices */ 1334 PetscCall(PetscMalloc1(nr, &s->m)); 1335 for (i = 0; i < nr; i++) PetscCall(PetscMalloc1(nc, &s->m[i])); 1336 for (i = 0; i < nr; i++) { 1337 for (j = 0; j < nc; j++) { 1338 s->m[i][j] = a[i * nc + j]; 1339 if (a[i * nc + j]) PetscCall(PetscObjectReference((PetscObject)a[i * nc + j])); 1340 } 1341 } 1342 PetscCall(MatGetVecType(A, &vtype)); 1343 PetscCall(PetscStrcmp(vtype, VECSTANDARD, &isstd)); 1344 if (isstd) { 1345 /* check if all blocks have the same vectype */ 1346 vtype = NULL; 1347 for (i = 0; i < nr; i++) { 1348 for (j = 0; j < nc; j++) { 1349 if (a[i * nc + j]) { 1350 if (!vtype) { /* first visited block */ 1351 PetscCall(MatGetVecType(a[i * nc + j], &vtype)); 1352 sametype = PETSC_TRUE; 1353 } else if (sametype) { 1354 PetscCall(MatGetVecType(a[i * nc + j], &type)); 1355 PetscCall(PetscStrcmp(vtype, type, &sametype)); 1356 } 1357 } 1358 } 1359 } 1360 if (sametype) { /* propagate vectype */ 1361 PetscCall(MatSetVecType(A, vtype)); 1362 } 1363 } 1364 1365 PetscCall(MatSetUp_NestIS_Private(A, nr, is_row, nc, is_col)); 1366 1367 PetscCall(PetscMalloc1(nr, &s->row_len)); 1368 PetscCall(PetscMalloc1(nc, &s->col_len)); 1369 for (i = 0; i < nr; i++) s->row_len[i] = -1; 1370 for (j = 0; j < nc; j++) s->col_len[j] = -1; 1371 1372 PetscCall(PetscCalloc1(nr * nc, &s->nnzstate)); 1373 for (i = 0; i < nr; i++) { 1374 for (j = 0; j < nc; j++) { 1375 if (s->m[i][j]) PetscCall(MatGetNonzeroState(s->m[i][j], &s->nnzstate[i * nc + j])); 1376 } 1377 } 1378 1379 PetscCall(MatNestGetSizes_Private(A, &m, &n, &M, &N)); 1380 1381 PetscCall(PetscLayoutSetSize(A->rmap, M)); 1382 PetscCall(PetscLayoutSetLocalSize(A->rmap, m)); 1383 PetscCall(PetscLayoutSetSize(A->cmap, N)); 1384 PetscCall(PetscLayoutSetLocalSize(A->cmap, n)); 1385 1386 PetscCall(PetscLayoutSetUp(A->rmap)); 1387 PetscCall(PetscLayoutSetUp(A->cmap)); 1388 1389 /* disable operations that are not supported for non-square matrices, 1390 or matrices for which is_row != is_col */ 1391 PetscCall(MatHasCongruentLayouts(A, &cong)); 1392 if (cong && nr != nc) cong = PETSC_FALSE; 1393 if (cong) { 1394 for (i = 0; cong && i < nr; i++) PetscCall(ISEqualUnsorted(s->isglobal.row[i], s->isglobal.col[i], &cong)); 1395 } 1396 if (!cong) { 1397 A->ops->missingdiagonal = NULL; 1398 A->ops->getdiagonal = NULL; 1399 A->ops->shift = NULL; 1400 A->ops->diagonalset = NULL; 1401 } 1402 1403 PetscCall(PetscCalloc2(nr, &s->left, nc, &s->right)); 1404 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 1405 A->nonzerostate++; 1406 PetscFunctionReturn(PETSC_SUCCESS); 1407 } 1408 1409 /*@ 1410 MatNestSetSubMats - Sets the nested submatrices in a `MATNEST` 1411 1412 Collective 1413 1414 Input Parameters: 1415 + A - `MATNEST` matrix 1416 . nr - number of nested row blocks 1417 . is_row - index sets for each nested row block, or `NULL` to make contiguous 1418 . nc - number of nested column blocks 1419 . is_col - index sets for each nested column block, or `NULL` to make contiguous 1420 - a - array of nr*nc submatrices, empty submatrices can be passed using `NULL` 1421 1422 Level: advanced 1423 1424 Notes: 1425 This always resets any submatrix information previously set 1426 1427 In both C and Fortran, `a` must be a row-major order array containing the matrices. See 1428 `MatCreateNest()` for an example. 1429 1430 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, `MatNestGetSubMats()` 1431 @*/ 1432 PetscErrorCode MatNestSetSubMats(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[]) 1433 { 1434 PetscInt i; 1435 1436 PetscFunctionBegin; 1437 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1438 PetscCheck(nr >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of rows cannot be negative"); 1439 if (nr && is_row) { 1440 PetscAssertPointer(is_row, 3); 1441 for (i = 0; i < nr; i++) PetscValidHeaderSpecific(is_row[i], IS_CLASSID, 3); 1442 } 1443 PetscCheck(nc >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of columns cannot be negative"); 1444 if (nc && is_col) { 1445 PetscAssertPointer(is_col, 5); 1446 for (i = 0; i < nc; i++) PetscValidHeaderSpecific(is_col[i], IS_CLASSID, 5); 1447 } 1448 if (nr * nc > 0) PetscAssertPointer(a, 6); 1449 PetscUseMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a)); 1450 PetscFunctionReturn(PETSC_SUCCESS); 1451 } 1452 1453 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A, PetscInt n, const IS islocal[], const IS isglobal[], PetscBool colflg, ISLocalToGlobalMapping *ltog) 1454 { 1455 PetscBool flg; 1456 PetscInt i, j, m, mi, *ix; 1457 1458 PetscFunctionBegin; 1459 *ltog = NULL; 1460 for (i = 0, m = 0, flg = PETSC_FALSE; i < n; i++) { 1461 if (islocal[i]) { 1462 PetscCall(ISGetLocalSize(islocal[i], &mi)); 1463 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1464 } else { 1465 PetscCall(ISGetLocalSize(isglobal[i], &mi)); 1466 } 1467 m += mi; 1468 } 1469 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 1470 1471 PetscCall(PetscMalloc1(m, &ix)); 1472 for (i = 0, m = 0; i < n; i++) { 1473 ISLocalToGlobalMapping smap = NULL; 1474 Mat sub = NULL; 1475 PetscSF sf; 1476 PetscLayout map; 1477 const PetscInt *ix2; 1478 1479 if (!colflg) { 1480 PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub)); 1481 } else { 1482 PetscCall(MatNestFindNonzeroSubMatCol(A, i, &sub)); 1483 } 1484 if (sub) { 1485 if (!colflg) { 1486 PetscCall(MatGetLocalToGlobalMapping(sub, &smap, NULL)); 1487 } else { 1488 PetscCall(MatGetLocalToGlobalMapping(sub, NULL, &smap)); 1489 } 1490 } 1491 /* 1492 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1493 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1494 */ 1495 PetscCall(ISGetIndices(isglobal[i], &ix2)); 1496 if (islocal[i]) { 1497 PetscInt *ilocal, *iremote; 1498 PetscInt mil, nleaves; 1499 1500 PetscCall(ISGetLocalSize(islocal[i], &mi)); 1501 PetscCheck(smap, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing local to global map"); 1502 for (j = 0; j < mi; j++) ix[m + j] = j; 1503 PetscCall(ISLocalToGlobalMappingApply(smap, mi, ix + m, ix + m)); 1504 1505 /* PetscSFSetGraphLayout does not like negative indices */ 1506 PetscCall(PetscMalloc2(mi, &ilocal, mi, &iremote)); 1507 for (j = 0, nleaves = 0; j < mi; j++) { 1508 if (ix[m + j] < 0) continue; 1509 ilocal[nleaves] = j; 1510 iremote[nleaves] = ix[m + j]; 1511 nleaves++; 1512 } 1513 PetscCall(ISGetLocalSize(isglobal[i], &mil)); 1514 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf)); 1515 PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)A), &map)); 1516 PetscCall(PetscLayoutSetLocalSize(map, mil)); 1517 PetscCall(PetscLayoutSetUp(map)); 1518 PetscCall(PetscSFSetGraphLayout(sf, map, nleaves, ilocal, PETSC_USE_POINTER, iremote)); 1519 PetscCall(PetscLayoutDestroy(&map)); 1520 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE)); 1521 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE)); 1522 PetscCall(PetscSFDestroy(&sf)); 1523 PetscCall(PetscFree2(ilocal, iremote)); 1524 } else { 1525 PetscCall(ISGetLocalSize(isglobal[i], &mi)); 1526 for (j = 0; j < mi; j++) ix[m + j] = ix2[i]; 1527 } 1528 PetscCall(ISRestoreIndices(isglobal[i], &ix2)); 1529 m += mi; 1530 } 1531 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, m, ix, PETSC_OWN_POINTER, ltog)); 1532 PetscFunctionReturn(PETSC_SUCCESS); 1533 } 1534 1535 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1536 /* 1537 nprocessors = NP 1538 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1539 proc 0: => (g_0,h_0,) 1540 proc 1: => (g_1,h_1,) 1541 ... 1542 proc nprocs-1: => (g_NP-1,h_NP-1,) 1543 1544 proc 0: proc 1: proc nprocs-1: 1545 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1546 1547 proc 0: 1548 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1549 proc 1: 1550 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1551 1552 proc NP-1: 1553 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1554 */ 1555 static PetscErrorCode MatSetUp_NestIS_Private(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[]) 1556 { 1557 Mat_Nest *vs = (Mat_Nest *)A->data; 1558 PetscInt i, j, offset, n, nsum, bs; 1559 Mat sub = NULL; 1560 1561 PetscFunctionBegin; 1562 PetscCall(PetscMalloc1(nr, &vs->isglobal.row)); 1563 PetscCall(PetscMalloc1(nc, &vs->isglobal.col)); 1564 if (is_row) { /* valid IS is passed in */ 1565 /* refs on is[] are incremented */ 1566 for (i = 0; i < vs->nr; i++) { 1567 PetscCall(PetscObjectReference((PetscObject)is_row[i])); 1568 1569 vs->isglobal.row[i] = is_row[i]; 1570 } 1571 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1572 nsum = 0; 1573 for (i = 0; i < vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1574 PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub)); 1575 PetscCheck(sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "No nonzero submatrix in row %" PetscInt_FMT, i); 1576 PetscCall(MatGetLocalSize(sub, &n, NULL)); 1577 PetscCheck(n >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Sizes have not yet been set for submatrix"); 1578 nsum += n; 1579 } 1580 PetscCallMPI(MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A))); 1581 offset -= nsum; 1582 for (i = 0; i < vs->nr; i++) { 1583 PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub)); 1584 PetscCall(MatGetLocalSize(sub, &n, NULL)); 1585 PetscCall(MatGetBlockSizes(sub, &bs, NULL)); 1586 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.row[i])); 1587 PetscCall(ISSetBlockSize(vs->isglobal.row[i], bs)); 1588 offset += n; 1589 } 1590 } 1591 1592 if (is_col) { /* valid IS is passed in */ 1593 /* refs on is[] are incremented */ 1594 for (j = 0; j < vs->nc; j++) { 1595 PetscCall(PetscObjectReference((PetscObject)is_col[j])); 1596 1597 vs->isglobal.col[j] = is_col[j]; 1598 } 1599 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1600 offset = A->cmap->rstart; 1601 nsum = 0; 1602 for (j = 0; j < vs->nc; j++) { 1603 PetscCall(MatNestFindNonzeroSubMatCol(A, j, &sub)); 1604 PetscCheck(sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "No nonzero submatrix in column %" PetscInt_FMT, i); 1605 PetscCall(MatGetLocalSize(sub, NULL, &n)); 1606 PetscCheck(n >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Sizes have not yet been set for submatrix"); 1607 nsum += n; 1608 } 1609 PetscCallMPI(MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A))); 1610 offset -= nsum; 1611 for (j = 0; j < vs->nc; j++) { 1612 PetscCall(MatNestFindNonzeroSubMatCol(A, j, &sub)); 1613 PetscCall(MatGetLocalSize(sub, NULL, &n)); 1614 PetscCall(MatGetBlockSizes(sub, NULL, &bs)); 1615 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.col[j])); 1616 PetscCall(ISSetBlockSize(vs->isglobal.col[j], bs)); 1617 offset += n; 1618 } 1619 } 1620 1621 /* Set up the local ISs */ 1622 PetscCall(PetscMalloc1(vs->nr, &vs->islocal.row)); 1623 PetscCall(PetscMalloc1(vs->nc, &vs->islocal.col)); 1624 for (i = 0, offset = 0; i < vs->nr; i++) { 1625 IS isloc; 1626 ISLocalToGlobalMapping rmap = NULL; 1627 PetscInt nlocal, bs; 1628 PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub)); 1629 if (sub) PetscCall(MatGetLocalToGlobalMapping(sub, &rmap, NULL)); 1630 if (rmap) { 1631 PetscCall(MatGetBlockSizes(sub, &bs, NULL)); 1632 PetscCall(ISLocalToGlobalMappingGetSize(rmap, &nlocal)); 1633 PetscCall(ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc)); 1634 PetscCall(ISSetBlockSize(isloc, bs)); 1635 } else { 1636 nlocal = 0; 1637 isloc = NULL; 1638 } 1639 vs->islocal.row[i] = isloc; 1640 offset += nlocal; 1641 } 1642 for (i = 0, offset = 0; i < vs->nc; i++) { 1643 IS isloc; 1644 ISLocalToGlobalMapping cmap = NULL; 1645 PetscInt nlocal, bs; 1646 PetscCall(MatNestFindNonzeroSubMatCol(A, i, &sub)); 1647 if (sub) PetscCall(MatGetLocalToGlobalMapping(sub, NULL, &cmap)); 1648 if (cmap) { 1649 PetscCall(MatGetBlockSizes(sub, NULL, &bs)); 1650 PetscCall(ISLocalToGlobalMappingGetSize(cmap, &nlocal)); 1651 PetscCall(ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc)); 1652 PetscCall(ISSetBlockSize(isloc, bs)); 1653 } else { 1654 nlocal = 0; 1655 isloc = NULL; 1656 } 1657 vs->islocal.col[i] = isloc; 1658 offset += nlocal; 1659 } 1660 1661 /* Set up the aggregate ISLocalToGlobalMapping */ 1662 { 1663 ISLocalToGlobalMapping rmap, cmap; 1664 PetscCall(MatNestCreateAggregateL2G_Private(A, vs->nr, vs->islocal.row, vs->isglobal.row, PETSC_FALSE, &rmap)); 1665 PetscCall(MatNestCreateAggregateL2G_Private(A, vs->nc, vs->islocal.col, vs->isglobal.col, PETSC_TRUE, &cmap)); 1666 if (rmap && cmap) PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap)); 1667 PetscCall(ISLocalToGlobalMappingDestroy(&rmap)); 1668 PetscCall(ISLocalToGlobalMappingDestroy(&cmap)); 1669 } 1670 1671 if (PetscDefined(USE_DEBUG)) { 1672 for (i = 0; i < vs->nr; i++) { 1673 for (j = 0; j < vs->nc; j++) { 1674 PetscInt m, n, M, N, mi, ni, Mi, Ni; 1675 Mat B = vs->m[i][j]; 1676 if (!B) continue; 1677 PetscCall(MatGetSize(B, &M, &N)); 1678 PetscCall(MatGetLocalSize(B, &m, &n)); 1679 PetscCall(ISGetSize(vs->isglobal.row[i], &Mi)); 1680 PetscCall(ISGetSize(vs->isglobal.col[j], &Ni)); 1681 PetscCall(ISGetLocalSize(vs->isglobal.row[i], &mi)); 1682 PetscCall(ISGetLocalSize(vs->isglobal.col[j], &ni)); 1683 PetscCheck(M == Mi && N == Ni, PetscObjectComm((PetscObject)sub), PETSC_ERR_ARG_INCOMP, "Global sizes (%" PetscInt_FMT ",%" PetscInt_FMT ") of nested submatrix (%" PetscInt_FMT ",%" PetscInt_FMT ") do not agree with space defined by index sets (%" PetscInt_FMT ",%" PetscInt_FMT ")", M, N, i, j, Mi, Ni); 1684 PetscCheck(m == mi && n == ni, PetscObjectComm((PetscObject)sub), PETSC_ERR_ARG_INCOMP, "Local sizes (%" PetscInt_FMT ",%" PetscInt_FMT ") of nested submatrix (%" PetscInt_FMT ",%" PetscInt_FMT ") do not agree with space defined by index sets (%" PetscInt_FMT ",%" PetscInt_FMT ")", m, n, i, j, mi, ni); 1685 } 1686 } 1687 } 1688 1689 /* Set A->assembled if all non-null blocks are currently assembled */ 1690 for (i = 0; i < vs->nr; i++) { 1691 for (j = 0; j < vs->nc; j++) { 1692 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(PETSC_SUCCESS); 1693 } 1694 } 1695 A->assembled = PETSC_TRUE; 1696 PetscFunctionReturn(PETSC_SUCCESS); 1697 } 1698 1699 /*@C 1700 MatCreateNest - Creates a new `MATNEST` matrix containing several nested submatrices, each stored separately 1701 1702 Collective 1703 1704 Input Parameters: 1705 + comm - Communicator for the new `MATNEST` 1706 . nr - number of nested row blocks 1707 . is_row - index sets for each nested row block, or `NULL` to make contiguous 1708 . nc - number of nested column blocks 1709 . is_col - index sets for each nested column block, or `NULL` to make contiguous 1710 - a - array of nr*nc submatrices, empty submatrices can be passed using `NULL` 1711 1712 Output Parameter: 1713 . B - new matrix 1714 1715 Note: 1716 In both C and Fortran, `a` must be a row-major order array holding references to the matrices. 1717 For instance, to represent the matrix 1718 $\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22}\end{bmatrix}$ 1719 one should use `Mat a[4]={A11,A12,A21,A22}`. 1720 1721 Level: advanced 1722 1723 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `VecCreateNest()`, `DMCreateMatrix()`, `MatNestSetSubMat()`, 1724 `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatNestGetSize()`, 1725 `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()` 1726 @*/ 1727 PetscErrorCode MatCreateNest(MPI_Comm comm, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[], Mat *B) 1728 { 1729 Mat A; 1730 1731 PetscFunctionBegin; 1732 *B = NULL; 1733 PetscCall(MatCreate(comm, &A)); 1734 PetscCall(MatSetType(A, MATNEST)); 1735 A->preallocated = PETSC_TRUE; 1736 PetscCall(MatNestSetSubMats(A, nr, is_row, nc, is_col, a)); 1737 *B = A; 1738 PetscFunctionReturn(PETSC_SUCCESS); 1739 } 1740 1741 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1742 { 1743 Mat_Nest *nest = (Mat_Nest *)A->data; 1744 Mat *trans; 1745 PetscScalar **avv; 1746 PetscScalar *vv; 1747 PetscInt **aii, **ajj; 1748 PetscInt *ii, *jj, *ci; 1749 PetscInt nr, nc, nnz, i, j; 1750 PetscBool done; 1751 1752 PetscFunctionBegin; 1753 PetscCall(MatGetSize(A, &nr, &nc)); 1754 if (reuse == MAT_REUSE_MATRIX) { 1755 PetscInt rnr; 1756 1757 PetscCall(MatGetRowIJ(*newmat, 0, PETSC_FALSE, PETSC_FALSE, &rnr, (const PetscInt **)&ii, (const PetscInt **)&jj, &done)); 1758 PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatGetRowIJ"); 1759 PetscCheck(rnr == nr, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of rows"); 1760 PetscCall(MatSeqAIJGetArray(*newmat, &vv)); 1761 } 1762 /* extract CSR for nested SeqAIJ matrices */ 1763 nnz = 0; 1764 PetscCall(PetscCalloc4(nest->nr * nest->nc, &aii, nest->nr * nest->nc, &ajj, nest->nr * nest->nc, &avv, nest->nr * nest->nc, &trans)); 1765 for (i = 0; i < nest->nr; ++i) { 1766 for (j = 0; j < nest->nc; ++j) { 1767 Mat B = nest->m[i][j]; 1768 if (B) { 1769 PetscScalar *naa; 1770 PetscInt *nii, *njj, nnr; 1771 PetscBool istrans; 1772 1773 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans)); 1774 if (istrans) { 1775 Mat Bt; 1776 1777 PetscCall(MatTransposeGetMat(B, &Bt)); 1778 PetscCall(MatTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j])); 1779 B = trans[i * nest->nc + j]; 1780 } else { 1781 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans)); 1782 if (istrans) { 1783 Mat Bt; 1784 1785 PetscCall(MatHermitianTransposeGetMat(B, &Bt)); 1786 PetscCall(MatHermitianTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j])); 1787 B = trans[i * nest->nc + j]; 1788 } 1789 } 1790 PetscCall(MatGetRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&nii, (const PetscInt **)&njj, &done)); 1791 PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatGetRowIJ"); 1792 PetscCall(MatSeqAIJGetArray(B, &naa)); 1793 nnz += nii[nnr]; 1794 1795 aii[i * nest->nc + j] = nii; 1796 ajj[i * nest->nc + j] = njj; 1797 avv[i * nest->nc + j] = naa; 1798 } 1799 } 1800 } 1801 if (reuse != MAT_REUSE_MATRIX) { 1802 PetscCall(PetscMalloc1(nr + 1, &ii)); 1803 PetscCall(PetscMalloc1(nnz, &jj)); 1804 PetscCall(PetscMalloc1(nnz, &vv)); 1805 } else { 1806 PetscCheck(nnz == ii[nr], PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of nonzeros"); 1807 } 1808 1809 /* new row pointer */ 1810 PetscCall(PetscArrayzero(ii, nr + 1)); 1811 for (i = 0; i < nest->nr; ++i) { 1812 PetscInt ncr, rst; 1813 1814 PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL)); 1815 PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr)); 1816 for (j = 0; j < nest->nc; ++j) { 1817 if (aii[i * nest->nc + j]) { 1818 PetscInt *nii = aii[i * nest->nc + j]; 1819 PetscInt ir; 1820 1821 for (ir = rst; ir < ncr + rst; ++ir) { 1822 ii[ir + 1] += nii[1] - nii[0]; 1823 nii++; 1824 } 1825 } 1826 } 1827 } 1828 for (i = 0; i < nr; i++) ii[i + 1] += ii[i]; 1829 1830 /* construct CSR for the new matrix */ 1831 PetscCall(PetscCalloc1(nr, &ci)); 1832 for (i = 0; i < nest->nr; ++i) { 1833 PetscInt ncr, rst; 1834 1835 PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL)); 1836 PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr)); 1837 for (j = 0; j < nest->nc; ++j) { 1838 if (aii[i * nest->nc + j]) { 1839 PetscScalar *nvv = avv[i * nest->nc + j]; 1840 PetscInt *nii = aii[i * nest->nc + j]; 1841 PetscInt *njj = ajj[i * nest->nc + j]; 1842 PetscInt ir, cst; 1843 1844 PetscCall(ISStrideGetInfo(nest->isglobal.col[j], &cst, NULL)); 1845 for (ir = rst; ir < ncr + rst; ++ir) { 1846 PetscInt ij, rsize = nii[1] - nii[0], ist = ii[ir] + ci[ir]; 1847 1848 for (ij = 0; ij < rsize; ij++) { 1849 jj[ist + ij] = *njj + cst; 1850 vv[ist + ij] = *nvv; 1851 njj++; 1852 nvv++; 1853 } 1854 ci[ir] += rsize; 1855 nii++; 1856 } 1857 } 1858 } 1859 } 1860 PetscCall(PetscFree(ci)); 1861 1862 /* restore info */ 1863 for (i = 0; i < nest->nr; ++i) { 1864 for (j = 0; j < nest->nc; ++j) { 1865 Mat B = nest->m[i][j]; 1866 if (B) { 1867 PetscInt nnr = 0, k = i * nest->nc + j; 1868 1869 B = (trans[k] ? trans[k] : B); 1870 PetscCall(MatRestoreRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&aii[k], (const PetscInt **)&ajj[k], &done)); 1871 PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatRestoreRowIJ"); 1872 PetscCall(MatSeqAIJRestoreArray(B, &avv[k])); 1873 PetscCall(MatDestroy(&trans[k])); 1874 } 1875 } 1876 } 1877 PetscCall(PetscFree4(aii, ajj, avv, trans)); 1878 1879 /* finalize newmat */ 1880 if (reuse == MAT_INITIAL_MATRIX) { 1881 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, newmat)); 1882 } else if (reuse == MAT_INPLACE_MATRIX) { 1883 Mat B; 1884 1885 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, &B)); 1886 PetscCall(MatHeaderReplace(A, &B)); 1887 } 1888 PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY)); 1889 PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY)); 1890 { 1891 Mat_SeqAIJ *a = (Mat_SeqAIJ *)((*newmat)->data); 1892 a->free_a = PETSC_TRUE; 1893 a->free_ij = PETSC_TRUE; 1894 } 1895 PetscFunctionReturn(PETSC_SUCCESS); 1896 } 1897 1898 PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y, PetscScalar a, Mat X) 1899 { 1900 Mat_Nest *nest = (Mat_Nest *)X->data; 1901 PetscInt i, j, k, rstart; 1902 PetscBool flg; 1903 1904 PetscFunctionBegin; 1905 /* Fill by row */ 1906 for (j = 0; j < nest->nc; ++j) { 1907 /* Using global column indices and ISAllGather() is not scalable. */ 1908 IS bNis; 1909 PetscInt bN; 1910 const PetscInt *bNindices; 1911 PetscCall(ISAllGather(nest->isglobal.col[j], &bNis)); 1912 PetscCall(ISGetSize(bNis, &bN)); 1913 PetscCall(ISGetIndices(bNis, &bNindices)); 1914 for (i = 0; i < nest->nr; ++i) { 1915 Mat B = nest->m[i][j], D = NULL; 1916 PetscInt bm, br; 1917 const PetscInt *bmindices; 1918 if (!B) continue; 1919 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "")); 1920 if (flg) { 1921 PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D)); 1922 PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D)); 1923 PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D)); 1924 B = D; 1925 } 1926 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, "")); 1927 if (flg) { 1928 if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D)); 1929 else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D)); 1930 B = D; 1931 } 1932 PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm)); 1933 PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices)); 1934 PetscCall(MatGetOwnershipRange(B, &rstart, NULL)); 1935 for (br = 0; br < bm; ++br) { 1936 PetscInt row = bmindices[br], brncols, *cols; 1937 const PetscInt *brcols; 1938 const PetscScalar *brcoldata; 1939 PetscScalar *vals = NULL; 1940 PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, &brcoldata)); 1941 PetscCall(PetscMalloc1(brncols, &cols)); 1942 for (k = 0; k < brncols; k++) cols[k] = bNindices[brcols[k]]; 1943 /* 1944 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 1945 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 1946 */ 1947 if (a != 1.0) { 1948 PetscCall(PetscMalloc1(brncols, &vals)); 1949 for (k = 0; k < brncols; k++) vals[k] = a * brcoldata[k]; 1950 PetscCall(MatSetValues(Y, 1, &row, brncols, cols, vals, ADD_VALUES)); 1951 PetscCall(PetscFree(vals)); 1952 } else { 1953 PetscCall(MatSetValues(Y, 1, &row, brncols, cols, brcoldata, ADD_VALUES)); 1954 } 1955 PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, &brcoldata)); 1956 PetscCall(PetscFree(cols)); 1957 } 1958 if (D) PetscCall(MatDestroy(&D)); 1959 PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices)); 1960 } 1961 PetscCall(ISRestoreIndices(bNis, &bNindices)); 1962 PetscCall(ISDestroy(&bNis)); 1963 } 1964 PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY)); 1965 PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY)); 1966 PetscFunctionReturn(PETSC_SUCCESS); 1967 } 1968 1969 static PetscErrorCode MatConvert_Nest_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1970 { 1971 Mat_Nest *nest = (Mat_Nest *)A->data; 1972 PetscInt m, n, M, N, i, j, k, *dnnz, *onnz = NULL, rstart, cstart, cend; 1973 PetscMPIInt size; 1974 Mat C; 1975 1976 PetscFunctionBegin; 1977 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 1978 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 1979 PetscInt nf; 1980 PetscBool fast; 1981 1982 PetscCall(PetscStrcmp(newtype, MATAIJ, &fast)); 1983 if (!fast) PetscCall(PetscStrcmp(newtype, MATSEQAIJ, &fast)); 1984 for (i = 0; i < nest->nr && fast; ++i) { 1985 for (j = 0; j < nest->nc && fast; ++j) { 1986 Mat B = nest->m[i][j]; 1987 if (B) { 1988 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &fast)); 1989 if (!fast) { 1990 PetscBool istrans; 1991 1992 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans)); 1993 if (istrans) { 1994 Mat Bt; 1995 1996 PetscCall(MatTransposeGetMat(B, &Bt)); 1997 PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast)); 1998 } else { 1999 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans)); 2000 if (istrans) { 2001 Mat Bt; 2002 2003 PetscCall(MatHermitianTransposeGetMat(B, &Bt)); 2004 PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast)); 2005 } 2006 } 2007 } 2008 } 2009 } 2010 } 2011 for (i = 0, nf = 0; i < nest->nr && fast; ++i) { 2012 PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i], ISSTRIDE, &fast)); 2013 if (fast) { 2014 PetscInt f, s; 2015 2016 PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &f, &s)); 2017 if (f != nf || s != 1) { 2018 fast = PETSC_FALSE; 2019 } else { 2020 PetscCall(ISGetSize(nest->isglobal.row[i], &f)); 2021 nf += f; 2022 } 2023 } 2024 } 2025 for (i = 0, nf = 0; i < nest->nc && fast; ++i) { 2026 PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i], ISSTRIDE, &fast)); 2027 if (fast) { 2028 PetscInt f, s; 2029 2030 PetscCall(ISStrideGetInfo(nest->isglobal.col[i], &f, &s)); 2031 if (f != nf || s != 1) { 2032 fast = PETSC_FALSE; 2033 } else { 2034 PetscCall(ISGetSize(nest->isglobal.col[i], &f)); 2035 nf += f; 2036 } 2037 } 2038 } 2039 if (fast) { 2040 PetscCall(MatConvert_Nest_SeqAIJ_fast(A, newtype, reuse, newmat)); 2041 PetscFunctionReturn(PETSC_SUCCESS); 2042 } 2043 } 2044 PetscCall(MatGetSize(A, &M, &N)); 2045 PetscCall(MatGetLocalSize(A, &m, &n)); 2046 PetscCall(MatGetOwnershipRangeColumn(A, &cstart, &cend)); 2047 if (reuse == MAT_REUSE_MATRIX) C = *newmat; 2048 else { 2049 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); 2050 PetscCall(MatSetType(C, newtype)); 2051 PetscCall(MatSetSizes(C, m, n, M, N)); 2052 } 2053 PetscCall(PetscMalloc1(2 * m, &dnnz)); 2054 if (m) { 2055 onnz = dnnz + m; 2056 for (k = 0; k < m; k++) { 2057 dnnz[k] = 0; 2058 onnz[k] = 0; 2059 } 2060 } 2061 for (j = 0; j < nest->nc; ++j) { 2062 IS bNis; 2063 PetscInt bN; 2064 const PetscInt *bNindices; 2065 PetscBool flg; 2066 /* Using global column indices and ISAllGather() is not scalable. */ 2067 PetscCall(ISAllGather(nest->isglobal.col[j], &bNis)); 2068 PetscCall(ISGetSize(bNis, &bN)); 2069 PetscCall(ISGetIndices(bNis, &bNindices)); 2070 for (i = 0; i < nest->nr; ++i) { 2071 PetscSF bmsf; 2072 PetscSFNode *iremote; 2073 Mat B = nest->m[i][j], D = NULL; 2074 PetscInt bm, *sub_dnnz, *sub_onnz, br; 2075 const PetscInt *bmindices; 2076 if (!B) continue; 2077 PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm)); 2078 PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices)); 2079 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf)); 2080 PetscCall(PetscMalloc1(bm, &iremote)); 2081 PetscCall(PetscMalloc1(bm, &sub_dnnz)); 2082 PetscCall(PetscMalloc1(bm, &sub_onnz)); 2083 for (k = 0; k < bm; ++k) { 2084 sub_dnnz[k] = 0; 2085 sub_onnz[k] = 0; 2086 } 2087 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "")); 2088 if (flg) { 2089 PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D)); 2090 PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D)); 2091 PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D)); 2092 B = D; 2093 } 2094 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, "")); 2095 if (flg) { 2096 if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D)); 2097 else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D)); 2098 B = D; 2099 } 2100 /* 2101 Locate the owners for all of the locally-owned global row indices for this row block. 2102 These determine the roots of PetscSF used to communicate preallocation data to row owners. 2103 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 2104 */ 2105 PetscCall(MatGetOwnershipRange(B, &rstart, NULL)); 2106 for (br = 0; br < bm; ++br) { 2107 PetscInt row = bmindices[br], brncols, col; 2108 const PetscInt *brcols; 2109 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 2110 PetscMPIInt rowowner = 0; 2111 PetscCall(PetscLayoutFindOwnerIndex(A->rmap, row, &rowowner, &rowrel)); 2112 /* how many roots */ 2113 iremote[br].rank = rowowner; 2114 iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 2115 /* get nonzero pattern */ 2116 PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, NULL)); 2117 for (k = 0; k < brncols; k++) { 2118 col = bNindices[brcols[k]]; 2119 if (col >= A->cmap->range[rowowner] && col < A->cmap->range[rowowner + 1]) { 2120 sub_dnnz[br]++; 2121 } else { 2122 sub_onnz[br]++; 2123 } 2124 } 2125 PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, NULL)); 2126 } 2127 if (D) PetscCall(MatDestroy(&D)); 2128 PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices)); 2129 /* bsf will have to take care of disposing of bedges. */ 2130 PetscCall(PetscSFSetGraph(bmsf, m, bm, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 2131 PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM)); 2132 PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM)); 2133 PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM)); 2134 PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM)); 2135 PetscCall(PetscFree(sub_dnnz)); 2136 PetscCall(PetscFree(sub_onnz)); 2137 PetscCall(PetscSFDestroy(&bmsf)); 2138 } 2139 PetscCall(ISRestoreIndices(bNis, &bNindices)); 2140 PetscCall(ISDestroy(&bNis)); 2141 } 2142 /* Resize preallocation if overestimated */ 2143 for (i = 0; i < m; i++) { 2144 dnnz[i] = PetscMin(dnnz[i], A->cmap->n); 2145 onnz[i] = PetscMin(onnz[i], A->cmap->N - A->cmap->n); 2146 } 2147 PetscCall(MatSeqAIJSetPreallocation(C, 0, dnnz)); 2148 PetscCall(MatMPIAIJSetPreallocation(C, 0, dnnz, 0, onnz)); 2149 PetscCall(PetscFree(dnnz)); 2150 PetscCall(MatAXPY_Dense_Nest(C, 1.0, A)); 2151 if (reuse == MAT_INPLACE_MATRIX) { 2152 PetscCall(MatHeaderReplace(A, &C)); 2153 } else *newmat = C; 2154 PetscFunctionReturn(PETSC_SUCCESS); 2155 } 2156 2157 static PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 2158 { 2159 Mat B; 2160 PetscInt m, n, M, N; 2161 2162 PetscFunctionBegin; 2163 PetscCall(MatGetSize(A, &M, &N)); 2164 PetscCall(MatGetLocalSize(A, &m, &n)); 2165 if (reuse == MAT_REUSE_MATRIX) { 2166 B = *newmat; 2167 PetscCall(MatZeroEntries(B)); 2168 } else { 2169 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, PETSC_DECIDE, M, N, NULL, &B)); 2170 } 2171 PetscCall(MatAXPY_Dense_Nest(B, 1.0, A)); 2172 if (reuse == MAT_INPLACE_MATRIX) { 2173 PetscCall(MatHeaderReplace(A, &B)); 2174 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 2175 PetscFunctionReturn(PETSC_SUCCESS); 2176 } 2177 2178 static PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has) 2179 { 2180 Mat_Nest *bA = (Mat_Nest *)mat->data; 2181 MatOperation opAdd; 2182 PetscInt i, j, nr = bA->nr, nc = bA->nc; 2183 PetscBool flg; 2184 PetscFunctionBegin; 2185 2186 *has = PETSC_FALSE; 2187 if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) { 2188 opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD); 2189 for (j = 0; j < nc; j++) { 2190 for (i = 0; i < nr; i++) { 2191 if (!bA->m[i][j]) continue; 2192 PetscCall(MatHasOperation(bA->m[i][j], opAdd, &flg)); 2193 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 2194 } 2195 } 2196 } 2197 if (((void **)mat->ops)[op]) *has = PETSC_TRUE; 2198 PetscFunctionReturn(PETSC_SUCCESS); 2199 } 2200 2201 /*MC 2202 MATNEST - "nest" - Matrix type consisting of nested submatrices, each stored separately. 2203 2204 Level: intermediate 2205 2206 Notes: 2207 This matrix type permits scalable use of `PCFIELDSPLIT` and avoids the large memory costs of extracting submatrices. 2208 It allows the use of symmetric and block formats for parts of multi-physics simulations. 2209 It is usually used with `DMCOMPOSITE` and `DMCreateMatrix()` 2210 2211 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero 2212 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes 2213 than the nest matrix. 2214 2215 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, 2216 `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`, 2217 `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()` 2218 M*/ 2219 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 2220 { 2221 Mat_Nest *s; 2222 2223 PetscFunctionBegin; 2224 PetscCall(PetscNew(&s)); 2225 A->data = (void *)s; 2226 2227 s->nr = -1; 2228 s->nc = -1; 2229 s->m = NULL; 2230 s->splitassembly = PETSC_FALSE; 2231 2232 PetscCall(PetscMemzero(A->ops, sizeof(*A->ops))); 2233 2234 A->ops->mult = MatMult_Nest; 2235 A->ops->multadd = MatMultAdd_Nest; 2236 A->ops->multtranspose = MatMultTranspose_Nest; 2237 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 2238 A->ops->transpose = MatTranspose_Nest; 2239 A->ops->assemblybegin = MatAssemblyBegin_Nest; 2240 A->ops->assemblyend = MatAssemblyEnd_Nest; 2241 A->ops->zeroentries = MatZeroEntries_Nest; 2242 A->ops->copy = MatCopy_Nest; 2243 A->ops->axpy = MatAXPY_Nest; 2244 A->ops->duplicate = MatDuplicate_Nest; 2245 A->ops->createsubmatrix = MatCreateSubMatrix_Nest; 2246 A->ops->destroy = MatDestroy_Nest; 2247 A->ops->view = MatView_Nest; 2248 A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 2249 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 2250 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 2251 A->ops->getdiagonal = MatGetDiagonal_Nest; 2252 A->ops->diagonalscale = MatDiagonalScale_Nest; 2253 A->ops->scale = MatScale_Nest; 2254 A->ops->shift = MatShift_Nest; 2255 A->ops->diagonalset = MatDiagonalSet_Nest; 2256 A->ops->setrandom = MatSetRandom_Nest; 2257 A->ops->hasoperation = MatHasOperation_Nest; 2258 A->ops->missingdiagonal = MatMissingDiagonal_Nest; 2259 2260 A->spptr = NULL; 2261 A->assembled = PETSC_FALSE; 2262 2263 /* expose Nest api's */ 2264 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", MatNestGetSubMat_Nest)); 2265 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", MatNestSetSubMat_Nest)); 2266 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", MatNestGetSubMats_Nest)); 2267 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", MatNestGetSize_Nest)); 2268 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", MatNestGetISs_Nest)); 2269 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", MatNestGetLocalISs_Nest)); 2270 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", MatNestSetVecType_Nest)); 2271 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", MatNestSetSubMats_Nest)); 2272 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ)); 2273 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ)); 2274 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", MatConvert_Nest_AIJ)); 2275 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", MatConvert_Nest_IS)); 2276 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", MatConvert_Nest_Dense)); 2277 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", MatConvert_Nest_Dense)); 2278 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", MatProductSetFromOptions_Nest_Dense)); 2279 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", MatProductSetFromOptions_Nest_Dense)); 2280 2281 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATNEST)); 2282 PetscFunctionReturn(PETSC_SUCCESS); 2283 } 2284