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