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 PetscViewerFormat format; 932 933 PetscCall(PetscViewerGetFormat(viewer, &format)); 934 if (format == PETSC_VIEWER_ASCII_MATLAB) { 935 Mat T; 936 937 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &T)); 938 PetscCall(MatView(T, viewer)); 939 PetscCall(MatDestroy(&T)); 940 PetscFunctionReturn(PETSC_SUCCESS); 941 } 942 PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_view_nest_sub", &viewSub, NULL)); 943 PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix object:\n")); 944 PetscCall(PetscViewerASCIIPushTab(viewer)); 945 PetscCall(PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", bA->nr, bA->nc)); 946 947 PetscCall(PetscViewerASCIIPrintf(viewer, "MatNest structure:\n")); 948 for (i = 0; i < bA->nr; i++) { 949 for (j = 0; j < bA->nc; j++) { 950 MatType type; 951 char name[256] = "", prefix[256] = ""; 952 PetscInt NR, NC; 953 PetscBool isNest = PETSC_FALSE; 954 955 if (!bA->m[i][j]) { 956 PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL\n", i, j)); 957 continue; 958 } 959 PetscCall(MatGetSize(bA->m[i][j], &NR, &NC)); 960 PetscCall(MatGetType(bA->m[i][j], &type)); 961 if (((PetscObject)bA->m[i][j])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bA->m[i][j])->name)); 962 if (((PetscObject)bA->m[i][j])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bA->m[i][j])->prefix)); 963 PetscCall(PetscObjectTypeCompare((PetscObject)bA->m[i][j], MATNEST, &isNest)); 964 965 PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", i, j, name, prefix, type, NR, NC)); 966 967 if (isNest || viewSub) { 968 PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */ 969 PetscCall(MatView(bA->m[i][j], viewer)); 970 PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */ 971 } 972 } 973 } 974 PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */ 975 } 976 PetscFunctionReturn(PETSC_SUCCESS); 977 } 978 979 static PetscErrorCode MatZeroEntries_Nest(Mat A) 980 { 981 Mat_Nest *bA = (Mat_Nest *)A->data; 982 PetscInt i, j; 983 984 PetscFunctionBegin; 985 for (i = 0; i < bA->nr; i++) { 986 for (j = 0; j < bA->nc; j++) { 987 if (!bA->m[i][j]) continue; 988 PetscCall(MatZeroEntries(bA->m[i][j])); 989 } 990 } 991 PetscFunctionReturn(PETSC_SUCCESS); 992 } 993 994 static PetscErrorCode MatCopy_Nest(Mat A, Mat B, MatStructure str) 995 { 996 Mat_Nest *bA = (Mat_Nest *)A->data, *bB = (Mat_Nest *)B->data; 997 PetscInt i, j, nr = bA->nr, nc = bA->nc; 998 PetscBool nnzstate = PETSC_FALSE; 999 1000 PetscFunctionBegin; 1001 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); 1002 for (i = 0; i < nr; i++) { 1003 for (j = 0; j < nc; j++) { 1004 PetscObjectState subnnzstate = 0; 1005 if (bA->m[i][j] && bB->m[i][j]) { 1006 PetscCall(MatCopy(bA->m[i][j], bB->m[i][j], str)); 1007 PetscCall(MatGetNonzeroState(bB->m[i][j], &subnnzstate)); 1008 nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i * nc + j] != subnnzstate); 1009 bB->nnzstate[i * nc + j] = subnnzstate; 1010 } else if (bA->m[i][j]) { // bB->m[i][j] is NULL 1011 Mat M; 1012 1013 PetscCheck(str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT ". Use DIFFERENT_NONZERO_PATTERN or UNKNOWN_NONZERO_PATTERN", i, j); 1014 PetscCall(MatDuplicate(bA->m[i][j], MAT_COPY_VALUES, &M)); 1015 PetscCall(MatNestSetSubMat(B, i, j, M)); 1016 PetscCall(MatDestroy(&M)); 1017 } else if (bB->m[i][j]) { // bA->m[i][j] is NULL 1018 PetscCheck(str == DIFFERENT_NONZERO_PATTERN || str == SUBSET_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT ". Use DIFFERENT_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN or UNKNOWN_NONZERO_PATTERN", i, j); 1019 PetscCall(MatNestSetSubMat(B, i, j, NULL)); 1020 } 1021 } 1022 } 1023 if (nnzstate) B->nonzerostate++; 1024 PetscFunctionReturn(PETSC_SUCCESS); 1025 } 1026 1027 static PetscErrorCode MatAXPY_Nest(Mat Y, PetscScalar a, Mat X, MatStructure str) 1028 { 1029 Mat_Nest *bY = (Mat_Nest *)Y->data, *bX = (Mat_Nest *)X->data; 1030 PetscInt i, j, nr = bY->nr, nc = bY->nc; 1031 PetscBool nnzstate = PETSC_FALSE; 1032 1033 PetscFunctionBegin; 1034 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); 1035 for (i = 0; i < nr; i++) { 1036 for (j = 0; j < nc; j++) { 1037 PetscObjectState subnnzstate = 0; 1038 if (bY->m[i][j] && bX->m[i][j]) { 1039 PetscCall(MatAXPY(bY->m[i][j], a, bX->m[i][j], str)); 1040 } else if (bX->m[i][j]) { 1041 Mat M; 1042 1043 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); 1044 PetscCall(MatDuplicate(bX->m[i][j], MAT_COPY_VALUES, &M)); 1045 PetscCall(MatNestSetSubMat(Y, i, j, M)); 1046 PetscCall(MatDestroy(&M)); 1047 } 1048 if (bY->m[i][j]) PetscCall(MatGetNonzeroState(bY->m[i][j], &subnnzstate)); 1049 nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i * nc + j] != subnnzstate); 1050 bY->nnzstate[i * nc + j] = subnnzstate; 1051 } 1052 } 1053 if (nnzstate) Y->nonzerostate++; 1054 PetscFunctionReturn(PETSC_SUCCESS); 1055 } 1056 1057 static PetscErrorCode MatDuplicate_Nest(Mat A, MatDuplicateOption op, Mat *B) 1058 { 1059 Mat_Nest *bA = (Mat_Nest *)A->data; 1060 Mat *b; 1061 PetscInt i, j, nr = bA->nr, nc = bA->nc; 1062 1063 PetscFunctionBegin; 1064 PetscCall(PetscMalloc1(nr * nc, &b)); 1065 for (i = 0; i < nr; i++) { 1066 for (j = 0; j < nc; j++) { 1067 if (bA->m[i][j]) { 1068 PetscCall(MatDuplicate(bA->m[i][j], op, &b[i * nc + j])); 1069 } else { 1070 b[i * nc + j] = NULL; 1071 } 1072 } 1073 } 1074 PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nr, bA->isglobal.row, nc, bA->isglobal.col, b, B)); 1075 /* Give the new MatNest exclusive ownership */ 1076 for (i = 0; i < nr * nc; i++) PetscCall(MatDestroy(&b[i])); 1077 PetscCall(PetscFree(b)); 1078 1079 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY)); 1080 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY)); 1081 PetscFunctionReturn(PETSC_SUCCESS); 1082 } 1083 1084 /* nest api */ 1085 static PetscErrorCode MatNestGetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat *mat) 1086 { 1087 Mat_Nest *bA = (Mat_Nest *)A->data; 1088 1089 PetscFunctionBegin; 1090 PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1); 1091 PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1); 1092 *mat = bA->m[idxm][jdxm]; 1093 PetscFunctionReturn(PETSC_SUCCESS); 1094 } 1095 1096 /*@ 1097 MatNestGetSubMat - Returns a single, sub-matrix from a `MATNEST` 1098 1099 Not Collective 1100 1101 Input Parameters: 1102 + A - `MATNEST` matrix 1103 . idxm - index of the matrix within the nest matrix 1104 - jdxm - index of the matrix within the nest matrix 1105 1106 Output Parameter: 1107 . sub - matrix at index `idxm`, `jdxm` within the nest matrix 1108 1109 Level: developer 1110 1111 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestSetSubMat()`, 1112 `MatNestGetLocalISs()`, `MatNestGetISs()` 1113 @*/ 1114 PetscErrorCode MatNestGetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat *sub) 1115 { 1116 PetscFunctionBegin; 1117 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1118 PetscValidLogicalCollectiveInt(A, idxm, 2); 1119 PetscValidLogicalCollectiveInt(A, jdxm, 3); 1120 PetscAssertPointer(sub, 4); 1121 PetscUseMethod(A, "MatNestGetSubMat_C", (Mat, PetscInt, PetscInt, Mat *), (A, idxm, jdxm, sub)); 1122 PetscFunctionReturn(PETSC_SUCCESS); 1123 } 1124 1125 static PetscErrorCode MatNestSetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat mat) 1126 { 1127 Mat_Nest *bA = (Mat_Nest *)A->data; 1128 PetscInt m, n, M, N, mi, ni, Mi, Ni; 1129 1130 PetscFunctionBegin; 1131 PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1); 1132 PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1); 1133 if (mat) { 1134 PetscCall(MatGetLocalSize(mat, &m, &n)); 1135 PetscCall(MatGetSize(mat, &M, &N)); 1136 PetscCall(ISGetLocalSize(bA->isglobal.row[idxm], &mi)); 1137 PetscCall(ISGetSize(bA->isglobal.row[idxm], &Mi)); 1138 PetscCall(ISGetLocalSize(bA->isglobal.col[jdxm], &ni)); 1139 PetscCall(ISGetSize(bA->isglobal.col[jdxm], &Ni)); 1140 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); 1141 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); 1142 } 1143 1144 /* do not increase object state */ 1145 if (mat == bA->m[idxm][jdxm]) PetscFunctionReturn(PETSC_SUCCESS); 1146 1147 PetscCall(PetscObjectReference((PetscObject)mat)); 1148 PetscCall(MatDestroy(&bA->m[idxm][jdxm])); 1149 bA->m[idxm][jdxm] = mat; 1150 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 1151 if (mat) PetscCall(MatGetNonzeroState(mat, &bA->nnzstate[idxm * bA->nc + jdxm])); 1152 else bA->nnzstate[idxm * bA->nc + jdxm] = 0; 1153 A->nonzerostate++; 1154 PetscFunctionReturn(PETSC_SUCCESS); 1155 } 1156 1157 /*@ 1158 MatNestSetSubMat - Set a single submatrix in the `MATNEST` 1159 1160 Logically Collective 1161 1162 Input Parameters: 1163 + A - `MATNEST` matrix 1164 . idxm - index of the matrix within the nest matrix 1165 . jdxm - index of the matrix within the nest matrix 1166 - sub - matrix at index `idxm`, `jdxm` within the nest matrix 1167 1168 Level: developer 1169 1170 Notes: 1171 The new submatrix must have the same size and communicator as that block of the nest. 1172 1173 This increments the reference count of the submatrix. 1174 1175 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestSetSubMats()`, `MatNestGetSubMats()`, `MatNestGetLocalISs()`, `MatCreateNest()`, 1176 `MatNestGetSubMat()`, `MatNestGetISs()`, `MatNestGetSize()` 1177 @*/ 1178 PetscErrorCode MatNestSetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat sub) 1179 { 1180 PetscFunctionBegin; 1181 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1182 PetscValidLogicalCollectiveInt(A, idxm, 2); 1183 PetscValidLogicalCollectiveInt(A, jdxm, 3); 1184 if (sub) PetscValidHeaderSpecific(sub, MAT_CLASSID, 4); 1185 PetscTryMethod(A, "MatNestSetSubMat_C", (Mat, PetscInt, PetscInt, Mat), (A, idxm, jdxm, sub)); 1186 PetscFunctionReturn(PETSC_SUCCESS); 1187 } 1188 1189 static PetscErrorCode MatNestGetSubMats_Nest(Mat A, PetscInt *M, PetscInt *N, Mat ***mat) 1190 { 1191 Mat_Nest *bA = (Mat_Nest *)A->data; 1192 1193 PetscFunctionBegin; 1194 if (M) *M = bA->nr; 1195 if (N) *N = bA->nc; 1196 if (mat) *mat = bA->m; 1197 PetscFunctionReturn(PETSC_SUCCESS); 1198 } 1199 1200 /*@C 1201 MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a `MATNEST` matrix. 1202 1203 Not Collective 1204 1205 Input Parameter: 1206 . A - nest matrix 1207 1208 Output Parameters: 1209 + M - number of submatrix rows in the nest matrix 1210 . N - number of submatrix columns in the nest matrix 1211 - mat - array of matrices 1212 1213 Level: developer 1214 1215 Note: 1216 The user should not free the array `mat`. 1217 1218 Fortran Notes: 1219 This routine has a calling sequence `call MatNestGetSubMats(A, M, N, mat, ierr)` 1220 where the space allocated for the optional argument `mat` is assumed large enough (if provided). 1221 Matrices in `mat` are returned in row-major order, see `MatCreateNest()` for an example. 1222 1223 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatCreateNest()`, 1224 `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()` 1225 @*/ 1226 PetscErrorCode MatNestGetSubMats(Mat A, PetscInt *M, PetscInt *N, Mat ***mat) 1227 { 1228 PetscFunctionBegin; 1229 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1230 PetscUseMethod(A, "MatNestGetSubMats_C", (Mat, PetscInt *, PetscInt *, Mat ***), (A, M, N, mat)); 1231 PetscFunctionReturn(PETSC_SUCCESS); 1232 } 1233 1234 static PetscErrorCode MatNestGetSize_Nest(Mat A, PetscInt *M, PetscInt *N) 1235 { 1236 Mat_Nest *bA = (Mat_Nest *)A->data; 1237 1238 PetscFunctionBegin; 1239 if (M) *M = bA->nr; 1240 if (N) *N = bA->nc; 1241 PetscFunctionReturn(PETSC_SUCCESS); 1242 } 1243 1244 /*@ 1245 MatNestGetSize - Returns the size of the `MATNEST` matrix. 1246 1247 Not Collective 1248 1249 Input Parameter: 1250 . A - `MATNEST` matrix 1251 1252 Output Parameters: 1253 + M - number of rows in the nested mat 1254 - N - number of cols in the nested mat 1255 1256 Level: developer 1257 1258 Note: 1259 `size` refers to the number of submatrices in the row and column directions of the nested matrix 1260 1261 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestGetLocalISs()`, 1262 `MatNestGetISs()` 1263 @*/ 1264 PetscErrorCode MatNestGetSize(Mat A, PetscInt *M, PetscInt *N) 1265 { 1266 PetscFunctionBegin; 1267 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1268 PetscUseMethod(A, "MatNestGetSize_C", (Mat, PetscInt *, PetscInt *), (A, M, N)); 1269 PetscFunctionReturn(PETSC_SUCCESS); 1270 } 1271 1272 static PetscErrorCode MatNestGetISs_Nest(Mat A, IS rows[], IS cols[]) 1273 { 1274 Mat_Nest *vs = (Mat_Nest *)A->data; 1275 PetscInt i; 1276 1277 PetscFunctionBegin; 1278 if (rows) 1279 for (i = 0; i < vs->nr; i++) rows[i] = vs->isglobal.row[i]; 1280 if (cols) 1281 for (i = 0; i < vs->nc; i++) cols[i] = vs->isglobal.col[i]; 1282 PetscFunctionReturn(PETSC_SUCCESS); 1283 } 1284 1285 /*@ 1286 MatNestGetISs - Returns the index sets partitioning the row and column spaces of a `MATNEST` 1287 1288 Not Collective 1289 1290 Input Parameter: 1291 . A - `MATNEST` matrix 1292 1293 Output Parameters: 1294 + rows - array of row index sets (pass `NULL` to ignore) 1295 - cols - array of column index sets (pass `NULL` to ignore) 1296 1297 Level: advanced 1298 1299 Note: 1300 The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s. 1301 1302 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`, 1303 `MatCreateNest()`, `MatNestSetSubMats()` 1304 @*/ 1305 PetscErrorCode MatNestGetISs(Mat A, IS rows[], IS cols[]) 1306 { 1307 PetscFunctionBegin; 1308 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1309 PetscUseMethod(A, "MatNestGetISs_C", (Mat, IS[], IS[]), (A, rows, cols)); 1310 PetscFunctionReturn(PETSC_SUCCESS); 1311 } 1312 1313 static PetscErrorCode MatNestGetLocalISs_Nest(Mat A, IS rows[], IS cols[]) 1314 { 1315 Mat_Nest *vs = (Mat_Nest *)A->data; 1316 PetscInt i; 1317 1318 PetscFunctionBegin; 1319 if (rows) 1320 for (i = 0; i < vs->nr; i++) rows[i] = vs->islocal.row[i]; 1321 if (cols) 1322 for (i = 0; i < vs->nc; i++) cols[i] = vs->islocal.col[i]; 1323 PetscFunctionReturn(PETSC_SUCCESS); 1324 } 1325 1326 /*@ 1327 MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces of a `MATNEST` 1328 1329 Not Collective 1330 1331 Input Parameter: 1332 . A - `MATNEST` matrix 1333 1334 Output Parameters: 1335 + rows - array of row index sets (pass `NULL` to ignore) 1336 - cols - array of column index sets (pass `NULL` to ignore) 1337 1338 Level: advanced 1339 1340 Note: 1341 The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s. 1342 1343 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`, 1344 `MatNestSetSubMats()`, `MatNestSetSubMat()` 1345 @*/ 1346 PetscErrorCode MatNestGetLocalISs(Mat A, IS rows[], IS cols[]) 1347 { 1348 PetscFunctionBegin; 1349 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1350 PetscUseMethod(A, "MatNestGetLocalISs_C", (Mat, IS[], IS[]), (A, rows, cols)); 1351 PetscFunctionReturn(PETSC_SUCCESS); 1352 } 1353 1354 static PetscErrorCode MatNestSetVecType_Nest(Mat A, VecType vtype) 1355 { 1356 PetscBool flg; 1357 1358 PetscFunctionBegin; 1359 PetscCall(PetscStrcmp(vtype, VECNEST, &flg)); 1360 /* In reality, this only distinguishes VECNEST and "other" */ 1361 if (flg) A->ops->getvecs = MatCreateVecs_Nest; 1362 else A->ops->getvecs = NULL; 1363 PetscFunctionReturn(PETSC_SUCCESS); 1364 } 1365 1366 /*@ 1367 MatNestSetVecType - Sets the type of `Vec` returned by `MatCreateVecs()` 1368 1369 Not Collective 1370 1371 Input Parameters: 1372 + A - `MATNEST` matrix 1373 - vtype - `VecType` to use for creating vectors 1374 1375 Level: developer 1376 1377 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateVecs()`, `MatCreateNest()`, `VecType` 1378 @*/ 1379 PetscErrorCode MatNestSetVecType(Mat A, VecType vtype) 1380 { 1381 PetscFunctionBegin; 1382 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1383 PetscTryMethod(A, "MatNestSetVecType_C", (Mat, VecType), (A, vtype)); 1384 PetscFunctionReturn(PETSC_SUCCESS); 1385 } 1386 1387 static PetscErrorCode MatNestSetSubMats_Nest(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[]) 1388 { 1389 Mat_Nest *s = (Mat_Nest *)A->data; 1390 PetscInt i, j, m, n, M, N; 1391 PetscBool cong, isstd, sametype = PETSC_FALSE; 1392 VecType vtype, type; 1393 1394 PetscFunctionBegin; 1395 PetscCall(MatReset_Nest(A)); 1396 1397 s->nr = nr; 1398 s->nc = nc; 1399 1400 /* Create space for submatrices */ 1401 PetscCall(PetscMalloc1(nr, &s->m)); 1402 PetscCall(PetscMalloc1(nr * nc, &s->m[0])); 1403 for (i = 0; i < nr; i++) { 1404 s->m[i] = s->m[0] + i * nc; 1405 for (j = 0; j < nc; j++) { 1406 s->m[i][j] = a ? a[i * nc + j] : NULL; 1407 PetscCall(PetscObjectReference((PetscObject)s->m[i][j])); 1408 } 1409 } 1410 PetscCall(MatGetVecType(A, &vtype)); 1411 PetscCall(PetscStrcmp(vtype, VECSTANDARD, &isstd)); 1412 if (isstd) { 1413 /* check if all blocks have the same vectype */ 1414 vtype = NULL; 1415 for (i = 0; i < nr; i++) { 1416 for (j = 0; j < nc; j++) { 1417 if (s->m[i][j]) { 1418 if (!vtype) { /* first visited block */ 1419 PetscCall(MatGetVecType(s->m[i][j], &vtype)); 1420 sametype = PETSC_TRUE; 1421 } else if (sametype) { 1422 PetscCall(MatGetVecType(s->m[i][j], &type)); 1423 PetscCall(PetscStrcmp(vtype, type, &sametype)); 1424 } 1425 } 1426 } 1427 } 1428 if (sametype) { /* propagate vectype */ 1429 PetscCall(MatSetVecType(A, vtype)); 1430 } 1431 } 1432 1433 PetscCall(MatSetUp_NestIS_Private(A, nr, is_row, nc, is_col)); 1434 1435 PetscCall(PetscMalloc1(nr, &s->row_len)); 1436 PetscCall(PetscMalloc1(nc, &s->col_len)); 1437 for (i = 0; i < nr; i++) s->row_len[i] = -1; 1438 for (j = 0; j < nc; j++) s->col_len[j] = -1; 1439 1440 PetscCall(PetscCalloc1(nr * nc, &s->nnzstate)); 1441 for (i = 0; i < nr; i++) { 1442 for (j = 0; j < nc; j++) { 1443 if (s->m[i][j]) PetscCall(MatGetNonzeroState(s->m[i][j], &s->nnzstate[i * nc + j])); 1444 } 1445 } 1446 1447 PetscCall(MatNestGetSizes_Private(A, &m, &n, &M, &N)); 1448 1449 PetscCall(PetscLayoutSetSize(A->rmap, M)); 1450 PetscCall(PetscLayoutSetLocalSize(A->rmap, m)); 1451 PetscCall(PetscLayoutSetSize(A->cmap, N)); 1452 PetscCall(PetscLayoutSetLocalSize(A->cmap, n)); 1453 1454 PetscCall(PetscLayoutSetUp(A->rmap)); 1455 PetscCall(PetscLayoutSetUp(A->cmap)); 1456 1457 /* disable operations that are not supported for non-square matrices, 1458 or matrices for which is_row != is_col */ 1459 PetscCall(MatHasCongruentLayouts(A, &cong)); 1460 if (cong && nr != nc) cong = PETSC_FALSE; 1461 if (cong) { 1462 for (i = 0; cong && i < nr; i++) PetscCall(ISEqualUnsorted(s->isglobal.row[i], s->isglobal.col[i], &cong)); 1463 } 1464 if (!cong) { 1465 A->ops->missingdiagonal = NULL; 1466 A->ops->getdiagonal = NULL; 1467 A->ops->shift = NULL; 1468 A->ops->diagonalset = NULL; 1469 } 1470 1471 PetscCall(PetscCalloc2(nr, &s->left, nc, &s->right)); 1472 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 1473 A->nonzerostate++; 1474 PetscFunctionReturn(PETSC_SUCCESS); 1475 } 1476 1477 /*@ 1478 MatNestSetSubMats - Sets the nested submatrices in a `MATNEST` 1479 1480 Collective 1481 1482 Input Parameters: 1483 + A - `MATNEST` matrix 1484 . nr - number of nested row blocks 1485 . is_row - index sets for each nested row block, or `NULL` to make contiguous 1486 . nc - number of nested column blocks 1487 . is_col - index sets for each nested column block, or `NULL` to make contiguous 1488 - a - array of nr*nc submatrices, or `NULL` 1489 1490 Level: advanced 1491 1492 Notes: 1493 This always resets any block matrix information previously set. 1494 1495 Pass `NULL` in the corresponding entry of `a` for an empty block. 1496 1497 In both C and Fortran, `a` must be a row-major order array containing the matrices. See 1498 `MatCreateNest()` for an example. 1499 1500 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, `MatNestGetSubMats()` 1501 @*/ 1502 PetscErrorCode MatNestSetSubMats(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[]) 1503 { 1504 PetscFunctionBegin; 1505 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1506 PetscValidLogicalCollectiveInt(A, nr, 2); 1507 PetscCheck(nr >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of rows cannot be negative"); 1508 if (nr && is_row) { 1509 PetscAssertPointer(is_row, 3); 1510 for (PetscInt i = 0; i < nr; i++) PetscValidHeaderSpecific(is_row[i], IS_CLASSID, 3); 1511 } 1512 PetscValidLogicalCollectiveInt(A, nc, 4); 1513 PetscCheck(nc >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of columns cannot be negative"); 1514 if (nc && is_col) { 1515 PetscAssertPointer(is_col, 5); 1516 for (PetscInt i = 0; i < nc; i++) PetscValidHeaderSpecific(is_col[i], IS_CLASSID, 5); 1517 } 1518 PetscTryMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a)); 1519 PetscFunctionReturn(PETSC_SUCCESS); 1520 } 1521 1522 static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A, PetscInt n, const IS islocal[], const IS isglobal[], PetscBool colflg, ISLocalToGlobalMapping *ltog) 1523 { 1524 PetscBool flg; 1525 PetscInt i, j, m, mi, *ix; 1526 1527 PetscFunctionBegin; 1528 *ltog = NULL; 1529 for (i = 0, m = 0, flg = PETSC_FALSE; i < n; i++) { 1530 if (islocal[i]) { 1531 PetscCall(ISGetLocalSize(islocal[i], &mi)); 1532 flg = PETSC_TRUE; /* We found a non-trivial entry */ 1533 } else { 1534 PetscCall(ISGetLocalSize(isglobal[i], &mi)); 1535 } 1536 m += mi; 1537 } 1538 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 1539 1540 PetscCall(PetscMalloc1(m, &ix)); 1541 for (i = 0, m = 0; i < n; i++) { 1542 ISLocalToGlobalMapping smap = NULL; 1543 Mat sub = NULL; 1544 PetscSF sf; 1545 PetscLayout map; 1546 const PetscInt *ix2; 1547 1548 if (!colflg) { 1549 PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub)); 1550 } else { 1551 PetscCall(MatNestFindNonzeroSubMatCol(A, i, &sub)); 1552 } 1553 if (sub) { 1554 if (!colflg) { 1555 PetscCall(MatGetLocalToGlobalMapping(sub, &smap, NULL)); 1556 } else { 1557 PetscCall(MatGetLocalToGlobalMapping(sub, NULL, &smap)); 1558 } 1559 } 1560 /* 1561 Now we need to extract the monolithic global indices that correspond to the given split global indices. 1562 In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces. 1563 */ 1564 PetscCall(ISGetIndices(isglobal[i], &ix2)); 1565 if (islocal[i]) { 1566 PetscInt *ilocal, *iremote; 1567 PetscInt mil, nleaves; 1568 1569 PetscCall(ISGetLocalSize(islocal[i], &mi)); 1570 PetscCheck(smap, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing local to global map"); 1571 for (j = 0; j < mi; j++) ix[m + j] = j; 1572 PetscCall(ISLocalToGlobalMappingApply(smap, mi, ix + m, ix + m)); 1573 1574 /* PetscSFSetGraphLayout does not like negative indices */ 1575 PetscCall(PetscMalloc2(mi, &ilocal, mi, &iremote)); 1576 for (j = 0, nleaves = 0; j < mi; j++) { 1577 if (ix[m + j] < 0) continue; 1578 ilocal[nleaves] = j; 1579 iremote[nleaves] = ix[m + j]; 1580 nleaves++; 1581 } 1582 PetscCall(ISGetLocalSize(isglobal[i], &mil)); 1583 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &sf)); 1584 PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)A), &map)); 1585 PetscCall(PetscLayoutSetLocalSize(map, mil)); 1586 PetscCall(PetscLayoutSetUp(map)); 1587 PetscCall(PetscSFSetGraphLayout(sf, map, nleaves, ilocal, PETSC_USE_POINTER, iremote)); 1588 PetscCall(PetscLayoutDestroy(&map)); 1589 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE)); 1590 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE)); 1591 PetscCall(PetscSFDestroy(&sf)); 1592 PetscCall(PetscFree2(ilocal, iremote)); 1593 } else { 1594 PetscCall(ISGetLocalSize(isglobal[i], &mi)); 1595 for (j = 0; j < mi; j++) ix[m + j] = ix2[i]; 1596 } 1597 PetscCall(ISRestoreIndices(isglobal[i], &ix2)); 1598 m += mi; 1599 } 1600 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, m, ix, PETSC_OWN_POINTER, ltog)); 1601 PetscFunctionReturn(PETSC_SUCCESS); 1602 } 1603 1604 /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */ 1605 /* 1606 nprocessors = NP 1607 Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1)) 1608 proc 0: => (g_0,h_0,) 1609 proc 1: => (g_1,h_1,) 1610 ... 1611 proc nprocs-1: => (g_NP-1,h_NP-1,) 1612 1613 proc 0: proc 1: proc nprocs-1: 1614 is[0] = (0,1,2,...,nlocal(g_0)-1) (0,1,...,nlocal(g_1)-1) (0,1,...,nlocal(g_NP-1)) 1615 1616 proc 0: 1617 is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1) 1618 proc 1: 1619 is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1) 1620 1621 proc NP-1: 1622 is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1) 1623 */ 1624 static PetscErrorCode MatSetUp_NestIS_Private(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[]) 1625 { 1626 Mat_Nest *vs = (Mat_Nest *)A->data; 1627 PetscInt i, j, offset, n, nsum, bs; 1628 Mat sub = NULL; 1629 1630 PetscFunctionBegin; 1631 PetscCall(PetscMalloc1(nr, &vs->isglobal.row)); 1632 PetscCall(PetscMalloc1(nc, &vs->isglobal.col)); 1633 if (is_row) { /* valid IS is passed in */ 1634 /* refs on is[] are incremented */ 1635 for (i = 0; i < vs->nr; i++) { 1636 PetscCall(PetscObjectReference((PetscObject)is_row[i])); 1637 vs->isglobal.row[i] = is_row[i]; 1638 } 1639 } else { /* Create the ISs by inspecting sizes of a submatrix in each row */ 1640 nsum = 0; 1641 for (i = 0; i < vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */ 1642 PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub)); 1643 PetscCheck(sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "No nonzero submatrix in row %" PetscInt_FMT, i); 1644 PetscCall(MatGetLocalSize(sub, &n, NULL)); 1645 PetscCheck(n >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Sizes have not yet been set for submatrix"); 1646 nsum += n; 1647 } 1648 PetscCallMPI(MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A))); 1649 offset -= nsum; 1650 for (i = 0; i < vs->nr; i++) { 1651 PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub)); 1652 PetscCall(MatGetLocalSize(sub, &n, NULL)); 1653 PetscCall(MatGetBlockSizes(sub, &bs, NULL)); 1654 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.row[i])); 1655 PetscCall(ISSetBlockSize(vs->isglobal.row[i], bs)); 1656 offset += n; 1657 } 1658 } 1659 1660 if (is_col) { /* valid IS is passed in */ 1661 /* refs on is[] are incremented */ 1662 for (j = 0; j < vs->nc; j++) { 1663 PetscCall(PetscObjectReference((PetscObject)is_col[j])); 1664 vs->isglobal.col[j] = is_col[j]; 1665 } 1666 } else { /* Create the ISs by inspecting sizes of a submatrix in each column */ 1667 offset = A->cmap->rstart; 1668 nsum = 0; 1669 for (j = 0; j < vs->nc; j++) { 1670 PetscCall(MatNestFindNonzeroSubMatCol(A, j, &sub)); 1671 PetscCheck(sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "No nonzero submatrix in column %" PetscInt_FMT, i); 1672 PetscCall(MatGetLocalSize(sub, NULL, &n)); 1673 PetscCheck(n >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Sizes have not yet been set for submatrix"); 1674 nsum += n; 1675 } 1676 PetscCallMPI(MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A))); 1677 offset -= nsum; 1678 for (j = 0; j < vs->nc; j++) { 1679 PetscCall(MatNestFindNonzeroSubMatCol(A, j, &sub)); 1680 PetscCall(MatGetLocalSize(sub, NULL, &n)); 1681 PetscCall(MatGetBlockSizes(sub, NULL, &bs)); 1682 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.col[j])); 1683 PetscCall(ISSetBlockSize(vs->isglobal.col[j], bs)); 1684 offset += n; 1685 } 1686 } 1687 1688 /* Set up the local ISs */ 1689 PetscCall(PetscMalloc1(vs->nr, &vs->islocal.row)); 1690 PetscCall(PetscMalloc1(vs->nc, &vs->islocal.col)); 1691 for (i = 0, offset = 0; i < vs->nr; i++) { 1692 IS isloc; 1693 ISLocalToGlobalMapping rmap = NULL; 1694 PetscInt nlocal, bs; 1695 PetscCall(MatNestFindNonzeroSubMatRow(A, i, &sub)); 1696 if (sub) PetscCall(MatGetLocalToGlobalMapping(sub, &rmap, NULL)); 1697 if (rmap) { 1698 PetscCall(MatGetBlockSizes(sub, &bs, NULL)); 1699 PetscCall(ISLocalToGlobalMappingGetSize(rmap, &nlocal)); 1700 PetscCall(ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc)); 1701 PetscCall(ISSetBlockSize(isloc, bs)); 1702 } else { 1703 nlocal = 0; 1704 isloc = NULL; 1705 } 1706 vs->islocal.row[i] = isloc; 1707 offset += nlocal; 1708 } 1709 for (i = 0, offset = 0; i < vs->nc; i++) { 1710 IS isloc; 1711 ISLocalToGlobalMapping cmap = NULL; 1712 PetscInt nlocal, bs; 1713 PetscCall(MatNestFindNonzeroSubMatCol(A, i, &sub)); 1714 if (sub) PetscCall(MatGetLocalToGlobalMapping(sub, NULL, &cmap)); 1715 if (cmap) { 1716 PetscCall(MatGetBlockSizes(sub, NULL, &bs)); 1717 PetscCall(ISLocalToGlobalMappingGetSize(cmap, &nlocal)); 1718 PetscCall(ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc)); 1719 PetscCall(ISSetBlockSize(isloc, bs)); 1720 } else { 1721 nlocal = 0; 1722 isloc = NULL; 1723 } 1724 vs->islocal.col[i] = isloc; 1725 offset += nlocal; 1726 } 1727 1728 /* Set up the aggregate ISLocalToGlobalMapping */ 1729 { 1730 ISLocalToGlobalMapping rmap, cmap; 1731 PetscCall(MatNestCreateAggregateL2G_Private(A, vs->nr, vs->islocal.row, vs->isglobal.row, PETSC_FALSE, &rmap)); 1732 PetscCall(MatNestCreateAggregateL2G_Private(A, vs->nc, vs->islocal.col, vs->isglobal.col, PETSC_TRUE, &cmap)); 1733 if (rmap && cmap) PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap)); 1734 PetscCall(ISLocalToGlobalMappingDestroy(&rmap)); 1735 PetscCall(ISLocalToGlobalMappingDestroy(&cmap)); 1736 } 1737 1738 if (PetscDefined(USE_DEBUG)) { 1739 for (i = 0; i < vs->nr; i++) { 1740 for (j = 0; j < vs->nc; j++) { 1741 PetscInt m, n, M, N, mi, ni, Mi, Ni; 1742 Mat B = vs->m[i][j]; 1743 if (!B) continue; 1744 PetscCall(MatGetSize(B, &M, &N)); 1745 PetscCall(MatGetLocalSize(B, &m, &n)); 1746 PetscCall(ISGetSize(vs->isglobal.row[i], &Mi)); 1747 PetscCall(ISGetSize(vs->isglobal.col[j], &Ni)); 1748 PetscCall(ISGetLocalSize(vs->isglobal.row[i], &mi)); 1749 PetscCall(ISGetLocalSize(vs->isglobal.col[j], &ni)); 1750 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); 1751 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); 1752 } 1753 } 1754 } 1755 1756 /* Set A->assembled if all non-null blocks are currently assembled */ 1757 for (i = 0; i < vs->nr; i++) { 1758 for (j = 0; j < vs->nc; j++) { 1759 if (vs->m[i][j] && !vs->m[i][j]->assembled) PetscFunctionReturn(PETSC_SUCCESS); 1760 } 1761 } 1762 A->assembled = PETSC_TRUE; 1763 PetscFunctionReturn(PETSC_SUCCESS); 1764 } 1765 1766 /*@C 1767 MatCreateNest - Creates a new `MATNEST` matrix containing several nested submatrices, each stored separately 1768 1769 Collective 1770 1771 Input Parameters: 1772 + comm - Communicator for the new `MATNEST` 1773 . nr - number of nested row blocks 1774 . is_row - index sets for each nested row block, or `NULL` to make contiguous 1775 . nc - number of nested column blocks 1776 . is_col - index sets for each nested column block, or `NULL` to make contiguous 1777 - a - array of nr*nc submatrices, empty submatrices can be passed using `NULL` 1778 1779 Output Parameter: 1780 . B - new matrix 1781 1782 Note: 1783 In both C and Fortran, `a` must be a row-major order array holding references to the matrices. 1784 For instance, to represent the matrix 1785 $\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22}\end{bmatrix}$ 1786 one should use `Mat a[4]={A11,A12,A21,A22}`. 1787 1788 Level: advanced 1789 1790 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `VecCreateNest()`, `DMCreateMatrix()`, `MatNestSetSubMat()`, 1791 `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatNestGetSize()`, 1792 `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()` 1793 @*/ 1794 PetscErrorCode MatCreateNest(MPI_Comm comm, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[], Mat *B) PeNSS 1795 { 1796 PetscFunctionBegin; 1797 PetscCall(MatCreate(comm, B)); 1798 PetscCall(MatSetType(*B, MATNEST)); 1799 (*B)->preallocated = PETSC_TRUE; 1800 PetscCall(MatNestSetSubMats(*B, nr, is_row, nc, is_col, a)); 1801 PetscFunctionReturn(PETSC_SUCCESS); 1802 } 1803 1804 static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 1805 { 1806 Mat_Nest *nest = (Mat_Nest *)A->data; 1807 Mat *trans; 1808 PetscScalar **avv; 1809 PetscScalar *vv; 1810 PetscInt **aii, **ajj; 1811 PetscInt *ii, *jj, *ci; 1812 PetscInt nr, nc, nnz, i, j; 1813 PetscBool done; 1814 1815 PetscFunctionBegin; 1816 PetscCall(MatGetSize(A, &nr, &nc)); 1817 if (reuse == MAT_REUSE_MATRIX) { 1818 PetscInt rnr; 1819 1820 PetscCall(MatGetRowIJ(*newmat, 0, PETSC_FALSE, PETSC_FALSE, &rnr, (const PetscInt **)&ii, (const PetscInt **)&jj, &done)); 1821 PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "MatGetRowIJ"); 1822 PetscCheck(rnr == nr, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of rows"); 1823 PetscCall(MatSeqAIJGetArray(*newmat, &vv)); 1824 } 1825 /* extract CSR for nested SeqAIJ matrices */ 1826 nnz = 0; 1827 PetscCall(PetscCalloc4(nest->nr * nest->nc, &aii, nest->nr * nest->nc, &ajj, nest->nr * nest->nc, &avv, nest->nr * nest->nc, &trans)); 1828 for (i = 0; i < nest->nr; ++i) { 1829 for (j = 0; j < nest->nc; ++j) { 1830 Mat B = nest->m[i][j]; 1831 if (B) { 1832 PetscScalar *naa; 1833 PetscInt *nii, *njj, nnr; 1834 PetscBool istrans; 1835 1836 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans)); 1837 if (istrans) { 1838 Mat Bt; 1839 1840 PetscCall(MatTransposeGetMat(B, &Bt)); 1841 PetscCall(MatTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j])); 1842 B = trans[i * nest->nc + j]; 1843 } else { 1844 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans)); 1845 if (istrans) { 1846 Mat Bt; 1847 1848 PetscCall(MatHermitianTransposeGetMat(B, &Bt)); 1849 PetscCall(MatHermitianTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j])); 1850 B = trans[i * nest->nc + j]; 1851 } 1852 } 1853 PetscCall(MatGetRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&nii, (const PetscInt **)&njj, &done)); 1854 PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatGetRowIJ"); 1855 PetscCall(MatSeqAIJGetArray(B, &naa)); 1856 nnz += nii[nnr]; 1857 1858 aii[i * nest->nc + j] = nii; 1859 ajj[i * nest->nc + j] = njj; 1860 avv[i * nest->nc + j] = naa; 1861 } 1862 } 1863 } 1864 if (reuse != MAT_REUSE_MATRIX) { 1865 PetscCall(PetscMalloc1(nr + 1, &ii)); 1866 PetscCall(PetscMalloc1(nnz, &jj)); 1867 PetscCall(PetscMalloc1(nnz, &vv)); 1868 } else { 1869 PetscCheck(nnz == ii[nr], PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of nonzeros"); 1870 } 1871 1872 /* new row pointer */ 1873 PetscCall(PetscArrayzero(ii, nr + 1)); 1874 for (i = 0; i < nest->nr; ++i) { 1875 PetscInt ncr, rst; 1876 1877 PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL)); 1878 PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr)); 1879 for (j = 0; j < nest->nc; ++j) { 1880 if (aii[i * nest->nc + j]) { 1881 PetscInt *nii = aii[i * nest->nc + j]; 1882 PetscInt ir; 1883 1884 for (ir = rst; ir < ncr + rst; ++ir) { 1885 ii[ir + 1] += nii[1] - nii[0]; 1886 nii++; 1887 } 1888 } 1889 } 1890 } 1891 for (i = 0; i < nr; i++) ii[i + 1] += ii[i]; 1892 1893 /* construct CSR for the new matrix */ 1894 PetscCall(PetscCalloc1(nr, &ci)); 1895 for (i = 0; i < nest->nr; ++i) { 1896 PetscInt ncr, rst; 1897 1898 PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL)); 1899 PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr)); 1900 for (j = 0; j < nest->nc; ++j) { 1901 if (aii[i * nest->nc + j]) { 1902 PetscScalar *nvv = avv[i * nest->nc + j], vscale = 1.0, vshift = 0.0; 1903 PetscInt *nii = aii[i * nest->nc + j]; 1904 PetscInt *njj = ajj[i * nest->nc + j]; 1905 PetscInt ir, cst; 1906 1907 if (trans[i * nest->nc + j]) { 1908 vscale = ((Mat_Shell *)nest->m[i][j]->data)->vscale; 1909 vshift = ((Mat_Shell *)nest->m[i][j]->data)->vshift; 1910 } 1911 PetscCall(ISStrideGetInfo(nest->isglobal.col[j], &cst, NULL)); 1912 for (ir = rst; ir < ncr + rst; ++ir) { 1913 PetscInt ij, rsize = nii[1] - nii[0], ist = ii[ir] + ci[ir]; 1914 1915 for (ij = 0; ij < rsize; ij++) { 1916 jj[ist + ij] = *njj + cst; 1917 vv[ist + ij] = vscale * *nvv; 1918 if (PetscUnlikely(vshift != 0.0 && *njj == ir - rst)) vv[ist + ij] += vshift; 1919 njj++; 1920 nvv++; 1921 } 1922 ci[ir] += rsize; 1923 nii++; 1924 } 1925 } 1926 } 1927 } 1928 PetscCall(PetscFree(ci)); 1929 1930 /* restore info */ 1931 for (i = 0; i < nest->nr; ++i) { 1932 for (j = 0; j < nest->nc; ++j) { 1933 Mat B = nest->m[i][j]; 1934 if (B) { 1935 PetscInt nnr = 0, k = i * nest->nc + j; 1936 1937 B = (trans[k] ? trans[k] : B); 1938 PetscCall(MatRestoreRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&aii[k], (const PetscInt **)&ajj[k], &done)); 1939 PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatRestoreRowIJ"); 1940 PetscCall(MatSeqAIJRestoreArray(B, &avv[k])); 1941 PetscCall(MatDestroy(&trans[k])); 1942 } 1943 } 1944 } 1945 PetscCall(PetscFree4(aii, ajj, avv, trans)); 1946 1947 /* finalize newmat */ 1948 if (reuse == MAT_INITIAL_MATRIX) { 1949 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, newmat)); 1950 } else if (reuse == MAT_INPLACE_MATRIX) { 1951 Mat B; 1952 1953 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, &B)); 1954 PetscCall(MatHeaderReplace(A, &B)); 1955 } 1956 PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY)); 1957 PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY)); 1958 { 1959 Mat_SeqAIJ *a = (Mat_SeqAIJ *)((*newmat)->data); 1960 a->free_a = PETSC_TRUE; 1961 a->free_ij = PETSC_TRUE; 1962 } 1963 PetscFunctionReturn(PETSC_SUCCESS); 1964 } 1965 1966 PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y, PetscScalar a, Mat X) 1967 { 1968 Mat_Nest *nest = (Mat_Nest *)X->data; 1969 PetscInt i, j, k, rstart; 1970 PetscBool flg; 1971 1972 PetscFunctionBegin; 1973 /* Fill by row */ 1974 for (j = 0; j < nest->nc; ++j) { 1975 /* Using global column indices and ISAllGather() is not scalable. */ 1976 IS bNis; 1977 PetscInt bN; 1978 const PetscInt *bNindices; 1979 PetscCall(ISAllGather(nest->isglobal.col[j], &bNis)); 1980 PetscCall(ISGetSize(bNis, &bN)); 1981 PetscCall(ISGetIndices(bNis, &bNindices)); 1982 for (i = 0; i < nest->nr; ++i) { 1983 Mat B = nest->m[i][j], D = NULL; 1984 PetscInt bm, br; 1985 const PetscInt *bmindices; 1986 if (!B) continue; 1987 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "")); 1988 if (flg) { 1989 PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D)); 1990 PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D)); 1991 PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D)); 1992 B = D; 1993 } 1994 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, "")); 1995 if (flg) { 1996 if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D)); 1997 else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D)); 1998 B = D; 1999 } 2000 PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm)); 2001 PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices)); 2002 PetscCall(MatGetOwnershipRange(B, &rstart, NULL)); 2003 for (br = 0; br < bm; ++br) { 2004 PetscInt row = bmindices[br], brncols, *cols; 2005 const PetscInt *brcols; 2006 const PetscScalar *brcoldata; 2007 PetscScalar *vals = NULL; 2008 PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, &brcoldata)); 2009 PetscCall(PetscMalloc1(brncols, &cols)); 2010 for (k = 0; k < brncols; k++) cols[k] = bNindices[brcols[k]]; 2011 /* 2012 Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match. 2013 Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES. 2014 */ 2015 if (a != 1.0) { 2016 PetscCall(PetscMalloc1(brncols, &vals)); 2017 for (k = 0; k < brncols; k++) vals[k] = a * brcoldata[k]; 2018 PetscCall(MatSetValues(Y, 1, &row, brncols, cols, vals, ADD_VALUES)); 2019 PetscCall(PetscFree(vals)); 2020 } else { 2021 PetscCall(MatSetValues(Y, 1, &row, brncols, cols, brcoldata, ADD_VALUES)); 2022 } 2023 PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, &brcoldata)); 2024 PetscCall(PetscFree(cols)); 2025 } 2026 if (D) PetscCall(MatDestroy(&D)); 2027 PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices)); 2028 } 2029 PetscCall(ISRestoreIndices(bNis, &bNindices)); 2030 PetscCall(ISDestroy(&bNis)); 2031 } 2032 PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY)); 2033 PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY)); 2034 PetscFunctionReturn(PETSC_SUCCESS); 2035 } 2036 2037 static PetscErrorCode MatConvert_Nest_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 2038 { 2039 Mat_Nest *nest = (Mat_Nest *)A->data; 2040 PetscInt m, n, M, N, i, j, k, *dnnz, *onnz = NULL, rstart, cstart, cend; 2041 PetscMPIInt size; 2042 Mat C; 2043 2044 PetscFunctionBegin; 2045 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 2046 if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */ 2047 PetscInt nf; 2048 PetscBool fast; 2049 2050 PetscCall(PetscStrcmp(newtype, MATAIJ, &fast)); 2051 if (!fast) PetscCall(PetscStrcmp(newtype, MATSEQAIJ, &fast)); 2052 for (i = 0; i < nest->nr && fast; ++i) { 2053 for (j = 0; j < nest->nc && fast; ++j) { 2054 Mat B = nest->m[i][j]; 2055 if (B) { 2056 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &fast)); 2057 if (!fast) { 2058 PetscBool istrans; 2059 2060 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans)); 2061 if (istrans) { 2062 Mat Bt; 2063 2064 PetscCall(MatTransposeGetMat(B, &Bt)); 2065 PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast)); 2066 } else { 2067 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans)); 2068 if (istrans) { 2069 Mat Bt; 2070 2071 PetscCall(MatHermitianTransposeGetMat(B, &Bt)); 2072 PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast)); 2073 } 2074 } 2075 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); 2076 } 2077 } 2078 } 2079 } 2080 for (i = 0, nf = 0; i < nest->nr && fast; ++i) { 2081 PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i], ISSTRIDE, &fast)); 2082 if (fast) { 2083 PetscInt f, s; 2084 2085 PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &f, &s)); 2086 if (f != nf || s != 1) { 2087 fast = PETSC_FALSE; 2088 } else { 2089 PetscCall(ISGetSize(nest->isglobal.row[i], &f)); 2090 nf += f; 2091 } 2092 } 2093 } 2094 for (i = 0, nf = 0; i < nest->nc && fast; ++i) { 2095 PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i], ISSTRIDE, &fast)); 2096 if (fast) { 2097 PetscInt f, s; 2098 2099 PetscCall(ISStrideGetInfo(nest->isglobal.col[i], &f, &s)); 2100 if (f != nf || s != 1) { 2101 fast = PETSC_FALSE; 2102 } else { 2103 PetscCall(ISGetSize(nest->isglobal.col[i], &f)); 2104 nf += f; 2105 } 2106 } 2107 } 2108 if (fast) { 2109 PetscCall(MatConvert_Nest_SeqAIJ_fast(A, newtype, reuse, newmat)); 2110 PetscFunctionReturn(PETSC_SUCCESS); 2111 } 2112 } 2113 PetscCall(MatGetSize(A, &M, &N)); 2114 PetscCall(MatGetLocalSize(A, &m, &n)); 2115 PetscCall(MatGetOwnershipRangeColumn(A, &cstart, &cend)); 2116 if (reuse == MAT_REUSE_MATRIX) C = *newmat; 2117 else { 2118 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C)); 2119 PetscCall(MatSetType(C, newtype)); 2120 PetscCall(MatSetSizes(C, m, n, M, N)); 2121 } 2122 PetscCall(PetscMalloc1(2 * m, &dnnz)); 2123 if (m) { 2124 onnz = dnnz + m; 2125 for (k = 0; k < m; k++) { 2126 dnnz[k] = 0; 2127 onnz[k] = 0; 2128 } 2129 } 2130 for (j = 0; j < nest->nc; ++j) { 2131 IS bNis; 2132 PetscInt bN; 2133 const PetscInt *bNindices; 2134 PetscBool flg; 2135 /* Using global column indices and ISAllGather() is not scalable. */ 2136 PetscCall(ISAllGather(nest->isglobal.col[j], &bNis)); 2137 PetscCall(ISGetSize(bNis, &bN)); 2138 PetscCall(ISGetIndices(bNis, &bNindices)); 2139 for (i = 0; i < nest->nr; ++i) { 2140 PetscSF bmsf; 2141 PetscSFNode *iremote; 2142 Mat B = nest->m[i][j], D = NULL; 2143 PetscInt bm, *sub_dnnz, *sub_onnz, br; 2144 const PetscInt *bmindices; 2145 if (!B) continue; 2146 PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm)); 2147 PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices)); 2148 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf)); 2149 PetscCall(PetscMalloc1(bm, &iremote)); 2150 PetscCall(PetscMalloc1(bm, &sub_dnnz)); 2151 PetscCall(PetscMalloc1(bm, &sub_onnz)); 2152 for (k = 0; k < bm; ++k) { 2153 sub_dnnz[k] = 0; 2154 sub_onnz[k] = 0; 2155 } 2156 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "")); 2157 if (flg) { 2158 PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D)); 2159 PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D)); 2160 PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D)); 2161 B = D; 2162 } 2163 PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, "")); 2164 if (flg) { 2165 if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D)); 2166 else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D)); 2167 B = D; 2168 } 2169 /* 2170 Locate the owners for all of the locally-owned global row indices for this row block. 2171 These determine the roots of PetscSF used to communicate preallocation data to row owners. 2172 The roots correspond to the dnnz and onnz entries; thus, there are two roots per row. 2173 */ 2174 PetscCall(MatGetOwnershipRange(B, &rstart, NULL)); 2175 for (br = 0; br < bm; ++br) { 2176 PetscInt row = bmindices[br], brncols, col; 2177 const PetscInt *brcols; 2178 PetscInt rowrel = 0; /* row's relative index on its owner rank */ 2179 PetscMPIInt rowowner = 0; 2180 PetscCall(PetscLayoutFindOwnerIndex(A->rmap, row, &rowowner, &rowrel)); 2181 /* how many roots */ 2182 iremote[br].rank = rowowner; 2183 iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */ 2184 /* get nonzero pattern */ 2185 PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, NULL)); 2186 for (k = 0; k < brncols; k++) { 2187 col = bNindices[brcols[k]]; 2188 if (col >= A->cmap->range[rowowner] && col < A->cmap->range[rowowner + 1]) { 2189 sub_dnnz[br]++; 2190 } else { 2191 sub_onnz[br]++; 2192 } 2193 } 2194 PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, NULL)); 2195 } 2196 if (D) PetscCall(MatDestroy(&D)); 2197 PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices)); 2198 /* bsf will have to take care of disposing of bedges. */ 2199 PetscCall(PetscSFSetGraph(bmsf, m, bm, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER)); 2200 PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM)); 2201 PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM)); 2202 PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM)); 2203 PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM)); 2204 PetscCall(PetscFree(sub_dnnz)); 2205 PetscCall(PetscFree(sub_onnz)); 2206 PetscCall(PetscSFDestroy(&bmsf)); 2207 } 2208 PetscCall(ISRestoreIndices(bNis, &bNindices)); 2209 PetscCall(ISDestroy(&bNis)); 2210 } 2211 /* Resize preallocation if overestimated */ 2212 for (i = 0; i < m; i++) { 2213 dnnz[i] = PetscMin(dnnz[i], A->cmap->n); 2214 onnz[i] = PetscMin(onnz[i], A->cmap->N - A->cmap->n); 2215 } 2216 PetscCall(MatSeqAIJSetPreallocation(C, 0, dnnz)); 2217 PetscCall(MatMPIAIJSetPreallocation(C, 0, dnnz, 0, onnz)); 2218 PetscCall(PetscFree(dnnz)); 2219 PetscCall(MatAXPY_Dense_Nest(C, 1.0, A)); 2220 if (reuse == MAT_INPLACE_MATRIX) { 2221 PetscCall(MatHeaderReplace(A, &C)); 2222 } else *newmat = C; 2223 PetscFunctionReturn(PETSC_SUCCESS); 2224 } 2225 2226 static PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 2227 { 2228 Mat B; 2229 PetscInt m, n, M, N; 2230 2231 PetscFunctionBegin; 2232 PetscCall(MatGetSize(A, &M, &N)); 2233 PetscCall(MatGetLocalSize(A, &m, &n)); 2234 if (reuse == MAT_REUSE_MATRIX) { 2235 B = *newmat; 2236 PetscCall(MatZeroEntries(B)); 2237 } else { 2238 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, PETSC_DECIDE, M, N, NULL, &B)); 2239 } 2240 PetscCall(MatAXPY_Dense_Nest(B, 1.0, A)); 2241 if (reuse == MAT_INPLACE_MATRIX) { 2242 PetscCall(MatHeaderReplace(A, &B)); 2243 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 2244 PetscFunctionReturn(PETSC_SUCCESS); 2245 } 2246 2247 static PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has) 2248 { 2249 Mat_Nest *bA = (Mat_Nest *)mat->data; 2250 MatOperation opAdd; 2251 PetscInt i, j, nr = bA->nr, nc = bA->nc; 2252 PetscBool flg; 2253 2254 PetscFunctionBegin; 2255 *has = PETSC_FALSE; 2256 if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) { 2257 opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD); 2258 for (j = 0; j < nc; j++) { 2259 for (i = 0; i < nr; i++) { 2260 if (!bA->m[i][j]) continue; 2261 PetscCall(MatHasOperation(bA->m[i][j], opAdd, &flg)); 2262 if (!flg) PetscFunctionReturn(PETSC_SUCCESS); 2263 } 2264 } 2265 } 2266 if (((void **)mat->ops)[op]) *has = PETSC_TRUE; 2267 PetscFunctionReturn(PETSC_SUCCESS); 2268 } 2269 2270 /*MC 2271 MATNEST - "nest" - Matrix type consisting of nested submatrices, each stored separately. 2272 2273 Level: intermediate 2274 2275 Notes: 2276 This matrix type permits scalable use of `PCFIELDSPLIT` and avoids the large memory costs of extracting submatrices. 2277 It allows the use of symmetric and block formats for parts of multi-physics simulations. 2278 It is usually used with `DMCOMPOSITE` and `DMCreateMatrix()` 2279 2280 Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero 2281 rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes 2282 than the nest matrix. 2283 2284 .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, 2285 `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`, 2286 `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()` 2287 M*/ 2288 PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A) 2289 { 2290 Mat_Nest *s; 2291 2292 PetscFunctionBegin; 2293 PetscCall(PetscNew(&s)); 2294 A->data = (void *)s; 2295 2296 s->nr = -1; 2297 s->nc = -1; 2298 s->m = NULL; 2299 s->splitassembly = PETSC_FALSE; 2300 2301 PetscCall(PetscMemzero(A->ops, sizeof(*A->ops))); 2302 2303 A->ops->mult = MatMult_Nest; 2304 A->ops->multadd = MatMultAdd_Nest; 2305 A->ops->multtranspose = MatMultTranspose_Nest; 2306 A->ops->multtransposeadd = MatMultTransposeAdd_Nest; 2307 A->ops->transpose = MatTranspose_Nest; 2308 A->ops->multhermitiantranspose = MatMultHermitianTranspose_Nest; 2309 A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_Nest; 2310 A->ops->assemblybegin = MatAssemblyBegin_Nest; 2311 A->ops->assemblyend = MatAssemblyEnd_Nest; 2312 A->ops->zeroentries = MatZeroEntries_Nest; 2313 A->ops->copy = MatCopy_Nest; 2314 A->ops->axpy = MatAXPY_Nest; 2315 A->ops->duplicate = MatDuplicate_Nest; 2316 A->ops->createsubmatrix = MatCreateSubMatrix_Nest; 2317 A->ops->destroy = MatDestroy_Nest; 2318 A->ops->view = MatView_Nest; 2319 A->ops->getvecs = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */ 2320 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_Nest; 2321 A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest; 2322 A->ops->getdiagonal = MatGetDiagonal_Nest; 2323 A->ops->diagonalscale = MatDiagonalScale_Nest; 2324 A->ops->scale = MatScale_Nest; 2325 A->ops->shift = MatShift_Nest; 2326 A->ops->diagonalset = MatDiagonalSet_Nest; 2327 A->ops->setrandom = MatSetRandom_Nest; 2328 A->ops->hasoperation = MatHasOperation_Nest; 2329 A->ops->missingdiagonal = MatMissingDiagonal_Nest; 2330 2331 A->spptr = NULL; 2332 A->assembled = PETSC_FALSE; 2333 2334 /* expose Nest api's */ 2335 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", MatNestGetSubMat_Nest)); 2336 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", MatNestSetSubMat_Nest)); 2337 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", MatNestGetSubMats_Nest)); 2338 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", MatNestGetSize_Nest)); 2339 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", MatNestGetISs_Nest)); 2340 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", MatNestGetLocalISs_Nest)); 2341 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", MatNestSetVecType_Nest)); 2342 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", MatNestSetSubMats_Nest)); 2343 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ)); 2344 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ)); 2345 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", MatConvert_Nest_AIJ)); 2346 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", MatConvert_Nest_IS)); 2347 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", MatConvert_Nest_Dense)); 2348 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", MatConvert_Nest_Dense)); 2349 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", MatProductSetFromOptions_Nest_Dense)); 2350 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", MatProductSetFromOptions_Nest_Dense)); 2351 2352 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATNEST)); 2353 PetscFunctionReturn(PETSC_SUCCESS); 2354 } 2355