1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 const char *const MatCompositeMergeTypes[] = {"left", "right", "MatCompositeMergeType", "MAT_COMPOSITE_", NULL}; 5 6 typedef struct _Mat_CompositeLink *Mat_CompositeLink; 7 struct _Mat_CompositeLink { 8 Mat mat; 9 Vec work; 10 Mat_CompositeLink next, prev; 11 }; 12 13 typedef struct { 14 MatCompositeType type; 15 Mat_CompositeLink head, tail; 16 Vec work; 17 PetscScalar scale; /* scale factor supplied with MatScale() */ 18 Vec left, right; /* left and right diagonal scaling provided with MatDiagonalScale() */ 19 Vec leftwork, rightwork, leftwork2, rightwork2; /* Two pairs of working vectors */ 20 PetscInt nmat; 21 PetscBool merge; 22 MatCompositeMergeType mergetype; 23 MatStructure structure; 24 25 PetscScalar *scalings; 26 PetscBool merge_mvctx; /* Whether need to merge mvctx of component matrices */ 27 Vec *lvecs; /* [nmat] Basically, they are Mvctx->lvec of each component matrix */ 28 PetscScalar *larray; /* [len] Data arrays of lvecs[] are stored consecutively in larray */ 29 PetscInt len; /* Length of larray[] */ 30 Vec gvec; /* Union of lvecs[] without duplicated entries */ 31 PetscInt *location; /* A map that maps entries in garray[] to larray[] */ 32 VecScatter Mvctx; 33 } Mat_Composite; 34 35 PetscErrorCode MatDestroy_Composite(Mat mat) { 36 Mat_Composite *shell = (Mat_Composite *)mat->data; 37 Mat_CompositeLink next = shell->head, oldnext; 38 PetscInt i; 39 40 PetscFunctionBegin; 41 while (next) { 42 PetscCall(MatDestroy(&next->mat)); 43 if (next->work && (!next->next || next->work != next->next->work)) PetscCall(VecDestroy(&next->work)); 44 oldnext = next; 45 next = next->next; 46 PetscCall(PetscFree(oldnext)); 47 } 48 PetscCall(VecDestroy(&shell->work)); 49 PetscCall(VecDestroy(&shell->left)); 50 PetscCall(VecDestroy(&shell->right)); 51 PetscCall(VecDestroy(&shell->leftwork)); 52 PetscCall(VecDestroy(&shell->rightwork)); 53 PetscCall(VecDestroy(&shell->leftwork2)); 54 PetscCall(VecDestroy(&shell->rightwork2)); 55 56 if (shell->Mvctx) { 57 for (i = 0; i < shell->nmat; i++) PetscCall(VecDestroy(&shell->lvecs[i])); 58 PetscCall(PetscFree3(shell->location, shell->larray, shell->lvecs)); 59 PetscCall(PetscFree(shell->larray)); 60 PetscCall(VecDestroy(&shell->gvec)); 61 PetscCall(VecScatterDestroy(&shell->Mvctx)); 62 } 63 64 PetscCall(PetscFree(shell->scalings)); 65 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeAddMat_C", NULL)); 66 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetType_C", NULL)); 67 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetType_C", NULL)); 68 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMergeType_C", NULL)); 69 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMatStructure_C", NULL)); 70 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMatStructure_C", NULL)); 71 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeMerge_C", NULL)); 72 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetNumberMat_C", NULL)); 73 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMat_C", NULL)); 74 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetScalings_C", NULL)); 75 PetscCall(PetscFree(mat->data)); 76 PetscFunctionReturn(0); 77 } 78 79 PetscErrorCode MatMult_Composite_Multiplicative(Mat A, Vec x, Vec y) { 80 Mat_Composite *shell = (Mat_Composite *)A->data; 81 Mat_CompositeLink next = shell->head; 82 Vec in, out; 83 PetscScalar scale; 84 PetscInt i; 85 86 PetscFunctionBegin; 87 PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()"); 88 in = x; 89 if (shell->right) { 90 if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork)); 91 PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in)); 92 in = shell->rightwork; 93 } 94 while (next->next) { 95 if (!next->work) { /* should reuse previous work if the same size */ 96 PetscCall(MatCreateVecs(next->mat, NULL, &next->work)); 97 } 98 out = next->work; 99 PetscCall(MatMult(next->mat, in, out)); 100 in = out; 101 next = next->next; 102 } 103 PetscCall(MatMult(next->mat, in, y)); 104 if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y)); 105 scale = shell->scale; 106 if (shell->scalings) { 107 for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i]; 108 } 109 PetscCall(VecScale(y, scale)); 110 PetscFunctionReturn(0); 111 } 112 113 PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A, Vec x, Vec y) { 114 Mat_Composite *shell = (Mat_Composite *)A->data; 115 Mat_CompositeLink tail = shell->tail; 116 Vec in, out; 117 PetscScalar scale; 118 PetscInt i; 119 120 PetscFunctionBegin; 121 PetscCheck(tail, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()"); 122 in = x; 123 if (shell->left) { 124 if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork)); 125 PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in)); 126 in = shell->leftwork; 127 } 128 while (tail->prev) { 129 if (!tail->prev->work) { /* should reuse previous work if the same size */ 130 PetscCall(MatCreateVecs(tail->mat, NULL, &tail->prev->work)); 131 } 132 out = tail->prev->work; 133 PetscCall(MatMultTranspose(tail->mat, in, out)); 134 in = out; 135 tail = tail->prev; 136 } 137 PetscCall(MatMultTranspose(tail->mat, in, y)); 138 if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y)); 139 140 scale = shell->scale; 141 if (shell->scalings) { 142 for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i]; 143 } 144 PetscCall(VecScale(y, scale)); 145 PetscFunctionReturn(0); 146 } 147 148 PetscErrorCode MatMult_Composite(Mat mat, Vec x, Vec y) { 149 Mat_Composite *shell = (Mat_Composite *)mat->data; 150 Mat_CompositeLink cur = shell->head; 151 Vec in, y2, xin; 152 Mat A, B; 153 PetscInt i, j, k, n, nuniq, lo, hi, mid, *gindices, *buf, *tmp, tot; 154 const PetscScalar *vals; 155 const PetscInt *garray; 156 IS ix, iy; 157 PetscBool match; 158 159 PetscFunctionBegin; 160 PetscCheck(cur, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()"); 161 in = x; 162 if (shell->right) { 163 if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork)); 164 PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in)); 165 in = shell->rightwork; 166 } 167 168 /* Try to merge Mvctx when instructed but not yet done. We did not do it in MatAssemblyEnd() since at that time 169 we did not know whether mat is ADDITIVE or MULTIPLICATIVE. Only now we are assured mat is ADDITIVE and 170 it is legal to merge Mvctx, because all component matrices have the same size. 171 */ 172 if (shell->merge_mvctx && !shell->Mvctx) { 173 /* Currently only implemented for MATMPIAIJ */ 174 for (cur = shell->head; cur; cur = cur->next) { 175 PetscCall(PetscObjectTypeCompare((PetscObject)cur->mat, MATMPIAIJ, &match)); 176 if (!match) { 177 shell->merge_mvctx = PETSC_FALSE; 178 goto skip_merge_mvctx; 179 } 180 } 181 182 /* Go through matrices first time to count total number of nonzero off-diag columns (may have dups) */ 183 tot = 0; 184 for (cur = shell->head; cur; cur = cur->next) { 185 PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, NULL)); 186 PetscCall(MatGetLocalSize(B, NULL, &n)); 187 tot += n; 188 } 189 PetscCall(PetscMalloc3(tot, &shell->location, tot, &shell->larray, shell->nmat, &shell->lvecs)); 190 shell->len = tot; 191 192 /* Go through matrices second time to sort off-diag columns and remove dups */ 193 PetscCall(PetscMalloc1(tot, &gindices)); /* No Malloc2() since we will give one to petsc and free the other */ 194 PetscCall(PetscMalloc1(tot, &buf)); 195 nuniq = 0; /* Number of unique nonzero columns */ 196 for (cur = shell->head; cur; cur = cur->next) { 197 PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray)); 198 PetscCall(MatGetLocalSize(B, NULL, &n)); 199 /* Merge pre-sorted garray[0,n) and gindices[0,nuniq) to buf[] */ 200 i = j = k = 0; 201 while (i < n && j < nuniq) { 202 if (garray[i] < gindices[j]) buf[k++] = garray[i++]; 203 else if (garray[i] > gindices[j]) buf[k++] = gindices[j++]; 204 else { 205 buf[k++] = garray[i++]; 206 j++; 207 } 208 } 209 /* Copy leftover in garray[] or gindices[] */ 210 if (i < n) { 211 PetscCall(PetscArraycpy(buf + k, garray + i, n - i)); 212 nuniq = k + n - i; 213 } else if (j < nuniq) { 214 PetscCall(PetscArraycpy(buf + k, gindices + j, nuniq - j)); 215 nuniq = k + nuniq - j; 216 } else nuniq = k; 217 /* Swap gindices and buf to merge garray of the next matrix */ 218 tmp = gindices; 219 gindices = buf; 220 buf = tmp; 221 } 222 PetscCall(PetscFree(buf)); 223 224 /* Go through matrices third time to build a map from gindices[] to garray[] */ 225 tot = 0; 226 for (cur = shell->head, j = 0; cur; cur = cur->next, j++) { /* j-th matrix */ 227 PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray)); 228 PetscCall(MatGetLocalSize(B, NULL, &n)); 229 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n, NULL, &shell->lvecs[j])); 230 /* This is an optimized PetscFindInt(garray[i],nuniq,gindices,&shell->location[tot+i]), using the fact that garray[] is also sorted */ 231 lo = 0; 232 for (i = 0; i < n; i++) { 233 hi = nuniq; 234 while (hi - lo > 1) { 235 mid = lo + (hi - lo) / 2; 236 if (garray[i] < gindices[mid]) hi = mid; 237 else lo = mid; 238 } 239 shell->location[tot + i] = lo; /* gindices[lo] = garray[i] */ 240 lo++; /* Since garray[i+1] > garray[i], we can safely advance lo */ 241 } 242 tot += n; 243 } 244 245 /* Build merged Mvctx */ 246 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nuniq, gindices, PETSC_OWN_POINTER, &ix)); 247 PetscCall(ISCreateStride(PETSC_COMM_SELF, nuniq, 0, 1, &iy)); 248 PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat), 1, mat->cmap->n, mat->cmap->N, NULL, &xin)); 249 PetscCall(VecCreateSeq(PETSC_COMM_SELF, nuniq, &shell->gvec)); 250 PetscCall(VecScatterCreate(xin, ix, shell->gvec, iy, &shell->Mvctx)); 251 PetscCall(VecDestroy(&xin)); 252 PetscCall(ISDestroy(&ix)); 253 PetscCall(ISDestroy(&iy)); 254 } 255 256 skip_merge_mvctx: 257 PetscCall(VecSet(y, 0)); 258 if (!shell->leftwork2) PetscCall(VecDuplicate(y, &shell->leftwork2)); 259 y2 = shell->leftwork2; 260 261 if (shell->Mvctx) { /* Have a merged Mvctx */ 262 /* Suppose we want to compute y = sMx, where s is the scaling factor and A, B are matrix M's diagonal/off-diagonal part. We could do 263 in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an oppertunity to 264 overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach. 265 */ 266 PetscCall(VecScatterBegin(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD)); 267 PetscCall(VecScatterEnd(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD)); 268 269 PetscCall(VecGetArrayRead(shell->gvec, &vals)); 270 for (i = 0; i < shell->len; i++) shell->larray[i] = vals[shell->location[i]]; 271 PetscCall(VecRestoreArrayRead(shell->gvec, &vals)); 272 273 for (cur = shell->head, tot = i = 0; cur; cur = cur->next, i++) { /* i-th matrix */ 274 PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, &A, &B, NULL)); 275 PetscUseTypeMethod(A, mult, in, y2); 276 PetscCall(MatGetLocalSize(B, NULL, &n)); 277 PetscCall(VecPlaceArray(shell->lvecs[i], &shell->larray[tot])); 278 PetscCall((*B->ops->multadd)(B, shell->lvecs[i], y2, y2)); 279 PetscCall(VecResetArray(shell->lvecs[i])); 280 PetscCall(VecAXPY(y, (shell->scalings ? shell->scalings[i] : 1.0), y2)); 281 tot += n; 282 } 283 } else { 284 if (shell->scalings) { 285 for (cur = shell->head, i = 0; cur; cur = cur->next, i++) { 286 PetscCall(MatMult(cur->mat, in, y2)); 287 PetscCall(VecAXPY(y, shell->scalings[i], y2)); 288 } 289 } else { 290 for (cur = shell->head; cur; cur = cur->next) PetscCall(MatMultAdd(cur->mat, in, y, y)); 291 } 292 } 293 294 if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y)); 295 PetscCall(VecScale(y, shell->scale)); 296 PetscFunctionReturn(0); 297 } 298 299 PetscErrorCode MatMultTranspose_Composite(Mat A, Vec x, Vec y) { 300 Mat_Composite *shell = (Mat_Composite *)A->data; 301 Mat_CompositeLink next = shell->head; 302 Vec in, y2 = NULL; 303 PetscInt i; 304 305 PetscFunctionBegin; 306 PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()"); 307 in = x; 308 if (shell->left) { 309 if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork)); 310 PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in)); 311 in = shell->leftwork; 312 } 313 314 PetscCall(MatMultTranspose(next->mat, in, y)); 315 if (shell->scalings) { 316 PetscCall(VecScale(y, shell->scalings[0])); 317 if (!shell->rightwork2) PetscCall(VecDuplicate(y, &shell->rightwork2)); 318 y2 = shell->rightwork2; 319 } 320 i = 1; 321 while ((next = next->next)) { 322 if (!shell->scalings) PetscCall(MatMultTransposeAdd(next->mat, in, y, y)); 323 else { 324 PetscCall(MatMultTranspose(next->mat, in, y2)); 325 PetscCall(VecAXPY(y, shell->scalings[i++], y2)); 326 } 327 } 328 if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y)); 329 PetscCall(VecScale(y, shell->scale)); 330 PetscFunctionReturn(0); 331 } 332 333 PetscErrorCode MatMultAdd_Composite(Mat A, Vec x, Vec y, Vec z) { 334 Mat_Composite *shell = (Mat_Composite *)A->data; 335 336 PetscFunctionBegin; 337 if (y != z) { 338 PetscCall(MatMult(A, x, z)); 339 PetscCall(VecAXPY(z, 1.0, y)); 340 } else { 341 if (!shell->leftwork) PetscCall(VecDuplicate(z, &shell->leftwork)); 342 PetscCall(MatMult(A, x, shell->leftwork)); 343 PetscCall(VecCopy(y, z)); 344 PetscCall(VecAXPY(z, 1.0, shell->leftwork)); 345 } 346 PetscFunctionReturn(0); 347 } 348 349 PetscErrorCode MatMultTransposeAdd_Composite(Mat A, Vec x, Vec y, Vec z) { 350 Mat_Composite *shell = (Mat_Composite *)A->data; 351 352 PetscFunctionBegin; 353 if (y != z) { 354 PetscCall(MatMultTranspose(A, x, z)); 355 PetscCall(VecAXPY(z, 1.0, y)); 356 } else { 357 if (!shell->rightwork) PetscCall(VecDuplicate(z, &shell->rightwork)); 358 PetscCall(MatMultTranspose(A, x, shell->rightwork)); 359 PetscCall(VecCopy(y, z)); 360 PetscCall(VecAXPY(z, 1.0, shell->rightwork)); 361 } 362 PetscFunctionReturn(0); 363 } 364 365 PetscErrorCode MatGetDiagonal_Composite(Mat A, Vec v) { 366 Mat_Composite *shell = (Mat_Composite *)A->data; 367 Mat_CompositeLink next = shell->head; 368 PetscInt i; 369 370 PetscFunctionBegin; 371 PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()"); 372 PetscCheck(!shell->right && !shell->left, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get diagonal if left or right scaling"); 373 374 PetscCall(MatGetDiagonal(next->mat, v)); 375 if (shell->scalings) PetscCall(VecScale(v, shell->scalings[0])); 376 377 if (next->next && !shell->work) PetscCall(VecDuplicate(v, &shell->work)); 378 i = 1; 379 while ((next = next->next)) { 380 PetscCall(MatGetDiagonal(next->mat, shell->work)); 381 PetscCall(VecAXPY(v, (shell->scalings ? shell->scalings[i++] : 1.0), shell->work)); 382 } 383 PetscCall(VecScale(v, shell->scale)); 384 PetscFunctionReturn(0); 385 } 386 387 PetscErrorCode MatAssemblyEnd_Composite(Mat Y, MatAssemblyType t) { 388 Mat_Composite *shell = (Mat_Composite *)Y->data; 389 390 PetscFunctionBegin; 391 if (shell->merge) PetscCall(MatCompositeMerge(Y)); 392 PetscFunctionReturn(0); 393 } 394 395 PetscErrorCode MatScale_Composite(Mat inA, PetscScalar alpha) { 396 Mat_Composite *a = (Mat_Composite *)inA->data; 397 398 PetscFunctionBegin; 399 a->scale *= alpha; 400 PetscFunctionReturn(0); 401 } 402 403 PetscErrorCode MatDiagonalScale_Composite(Mat inA, Vec left, Vec right) { 404 Mat_Composite *a = (Mat_Composite *)inA->data; 405 406 PetscFunctionBegin; 407 if (left) { 408 if (!a->left) { 409 PetscCall(VecDuplicate(left, &a->left)); 410 PetscCall(VecCopy(left, a->left)); 411 } else { 412 PetscCall(VecPointwiseMult(a->left, left, a->left)); 413 } 414 } 415 if (right) { 416 if (!a->right) { 417 PetscCall(VecDuplicate(right, &a->right)); 418 PetscCall(VecCopy(right, a->right)); 419 } else { 420 PetscCall(VecPointwiseMult(a->right, right, a->right)); 421 } 422 } 423 PetscFunctionReturn(0); 424 } 425 426 PetscErrorCode MatSetFromOptions_Composite(Mat A, PetscOptionItems *PetscOptionsObject) { 427 Mat_Composite *a = (Mat_Composite *)A->data; 428 429 PetscFunctionBegin; 430 PetscOptionsHeadBegin(PetscOptionsObject, "MATCOMPOSITE options"); 431 PetscCall(PetscOptionsBool("-mat_composite_merge", "Merge at MatAssemblyEnd", "MatCompositeMerge", a->merge, &a->merge, NULL)); 432 PetscCall(PetscOptionsEnum("-mat_composite_merge_type", "Set composite merge direction", "MatCompositeSetMergeType", MatCompositeMergeTypes, (PetscEnum)a->mergetype, (PetscEnum *)&a->mergetype, NULL)); 433 PetscCall(PetscOptionsBool("-mat_composite_merge_mvctx", "Merge MatMult() vecscat contexts", "MatCreateComposite", a->merge_mvctx, &a->merge_mvctx, NULL)); 434 PetscOptionsHeadEnd(); 435 PetscFunctionReturn(0); 436 } 437 438 /*@ 439 MatCreateComposite - Creates a matrix as the sum or product of one or more matrices 440 441 Collective 442 443 Input Parameters: 444 + comm - MPI communicator 445 . nmat - number of matrices to put in 446 - mats - the matrices 447 448 Output Parameter: 449 . mat - the matrix 450 451 Options Database Keys: 452 + -mat_composite_merge - merge in `MatAssemblyEnd()` 453 . -mat_composite_merge_mvctx - merge Mvctx of component matrices to optimize communication in `MatMult()` for ADDITIVE matrices 454 - -mat_composite_merge_type - set merge direction 455 456 Level: advanced 457 458 Note: 459 Alternative construction 460 $ MatCreate(comm,&mat); 461 $ MatSetSizes(mat,m,n,M,N); 462 $ MatSetType(mat,MATCOMPOSITE); 463 $ MatCompositeAddMat(mat,mats[0]); 464 $ .... 465 $ MatCompositeAddMat(mat,mats[nmat-1]); 466 $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 467 $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); 468 469 For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0] 470 471 .seealso: `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCompositeGetMat()`, `MatCompositeMerge()`, `MatCompositeSetType()`, `MATCOMPOSITE` 472 @*/ 473 PetscErrorCode MatCreateComposite(MPI_Comm comm, PetscInt nmat, const Mat *mats, Mat *mat) { 474 PetscInt m, n, M, N, i; 475 476 PetscFunctionBegin; 477 PetscCheck(nmat >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must pass in at least one matrix"); 478 PetscValidPointer(mat, 4); 479 480 PetscCall(MatGetLocalSize(mats[0], PETSC_IGNORE, &n)); 481 PetscCall(MatGetLocalSize(mats[nmat - 1], &m, PETSC_IGNORE)); 482 PetscCall(MatGetSize(mats[0], PETSC_IGNORE, &N)); 483 PetscCall(MatGetSize(mats[nmat - 1], &M, PETSC_IGNORE)); 484 PetscCall(MatCreate(comm, mat)); 485 PetscCall(MatSetSizes(*mat, m, n, M, N)); 486 PetscCall(MatSetType(*mat, MATCOMPOSITE)); 487 for (i = 0; i < nmat; i++) PetscCall(MatCompositeAddMat(*mat, mats[i])); 488 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 489 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 490 PetscFunctionReturn(0); 491 } 492 493 static PetscErrorCode MatCompositeAddMat_Composite(Mat mat, Mat smat) { 494 Mat_Composite *shell = (Mat_Composite *)mat->data; 495 Mat_CompositeLink ilink, next = shell->head; 496 497 PetscFunctionBegin; 498 PetscCall(PetscNew(&ilink)); 499 ilink->next = NULL; 500 PetscCall(PetscObjectReference((PetscObject)smat)); 501 ilink->mat = smat; 502 503 if (!next) shell->head = ilink; 504 else { 505 while (next->next) next = next->next; 506 next->next = ilink; 507 ilink->prev = next; 508 } 509 shell->tail = ilink; 510 shell->nmat += 1; 511 512 /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */ 513 if (shell->scalings) { 514 PetscCall(PetscRealloc(sizeof(PetscScalar) * shell->nmat, &shell->scalings)); 515 shell->scalings[shell->nmat - 1] = 1.0; 516 } 517 PetscFunctionReturn(0); 518 } 519 520 /*@ 521 MatCompositeAddMat - Add another matrix to a composite matrix. 522 523 Collective on mat 524 525 Input Parameters: 526 + mat - the composite matrix 527 - smat - the partial matrix 528 529 Level: advanced 530 531 .seealso: `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE` 532 @*/ 533 PetscErrorCode MatCompositeAddMat(Mat mat, Mat smat) { 534 PetscFunctionBegin; 535 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 536 PetscValidHeaderSpecific(smat, MAT_CLASSID, 2); 537 PetscUseMethod(mat, "MatCompositeAddMat_C", (Mat, Mat), (mat, smat)); 538 PetscFunctionReturn(0); 539 } 540 541 static PetscErrorCode MatCompositeSetType_Composite(Mat mat, MatCompositeType type) { 542 Mat_Composite *b = (Mat_Composite *)mat->data; 543 544 PetscFunctionBegin; 545 b->type = type; 546 if (type == MAT_COMPOSITE_MULTIPLICATIVE) { 547 mat->ops->getdiagonal = NULL; 548 mat->ops->mult = MatMult_Composite_Multiplicative; 549 mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative; 550 b->merge_mvctx = PETSC_FALSE; 551 } else { 552 mat->ops->getdiagonal = MatGetDiagonal_Composite; 553 mat->ops->mult = MatMult_Composite; 554 mat->ops->multtranspose = MatMultTranspose_Composite; 555 } 556 PetscFunctionReturn(0); 557 } 558 559 /*@ 560 MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product. 561 562 Logically Collective on mat 563 564 Input Parameters: 565 . mat - the composite matrix 566 567 Level: advanced 568 569 .seealso: `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeGetType()`, `MATCOMPOSITE` 570 @*/ 571 PetscErrorCode MatCompositeSetType(Mat mat, MatCompositeType type) { 572 PetscFunctionBegin; 573 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 574 PetscValidLogicalCollectiveEnum(mat, type, 2); 575 PetscUseMethod(mat, "MatCompositeSetType_C", (Mat, MatCompositeType), (mat, type)); 576 PetscFunctionReturn(0); 577 } 578 579 static PetscErrorCode MatCompositeGetType_Composite(Mat mat, MatCompositeType *type) { 580 Mat_Composite *b = (Mat_Composite *)mat->data; 581 582 PetscFunctionBegin; 583 *type = b->type; 584 PetscFunctionReturn(0); 585 } 586 587 /*@ 588 MatCompositeGetType - Returns type of composite. 589 590 Not Collective 591 592 Input Parameter: 593 . mat - the composite matrix 594 595 Output Parameter: 596 . type - type of composite 597 598 Level: advanced 599 600 .seealso: `MatCreateComposite()`, `MatCompositeSetType()`, `MATCOMPOSITE`, `MatCompositeType` 601 @*/ 602 PetscErrorCode MatCompositeGetType(Mat mat, MatCompositeType *type) { 603 PetscFunctionBegin; 604 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 605 PetscValidPointer(type, 2); 606 PetscUseMethod(mat, "MatCompositeGetType_C", (Mat, MatCompositeType *), (mat, type)); 607 PetscFunctionReturn(0); 608 } 609 610 static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat, MatStructure str) { 611 Mat_Composite *b = (Mat_Composite *)mat->data; 612 613 PetscFunctionBegin; 614 b->structure = str; 615 PetscFunctionReturn(0); 616 } 617 618 /*@ 619 MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix. 620 621 Not Collective 622 623 Input Parameters: 624 + mat - the composite matrix 625 - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN` (default) or `SUBSET_NONZERO_PATTERN` 626 627 Level: advanced 628 629 Note: 630 Information about the matrices structure is used in `MatCompositeMerge()` for additive composite matrix. 631 632 .seealso: `MatAXPY()`, `MatCreateComposite()`, `MatCompositeMerge()` `MatCompositeGetMatStructure()`, `MATCOMPOSITE` 633 @*/ 634 PetscErrorCode MatCompositeSetMatStructure(Mat mat, MatStructure str) { 635 PetscFunctionBegin; 636 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 637 PetscUseMethod(mat, "MatCompositeSetMatStructure_C", (Mat, MatStructure), (mat, str)); 638 PetscFunctionReturn(0); 639 } 640 641 static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat, MatStructure *str) { 642 Mat_Composite *b = (Mat_Composite *)mat->data; 643 644 PetscFunctionBegin; 645 *str = b->structure; 646 PetscFunctionReturn(0); 647 } 648 649 /*@ 650 MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix. 651 652 Not Collective 653 654 Input Parameter: 655 . mat - the composite matrix 656 657 Output Parameter: 658 . str - structure of the matrices 659 660 Level: advanced 661 662 .seealso: `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MATCOMPOSITE` 663 @*/ 664 PetscErrorCode MatCompositeGetMatStructure(Mat mat, MatStructure *str) { 665 PetscFunctionBegin; 666 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 667 PetscValidPointer(str, 2); 668 PetscUseMethod(mat, "MatCompositeGetMatStructure_C", (Mat, MatStructure *), (mat, str)); 669 PetscFunctionReturn(0); 670 } 671 672 static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat, MatCompositeMergeType type) { 673 Mat_Composite *shell = (Mat_Composite *)mat->data; 674 675 PetscFunctionBegin; 676 shell->mergetype = type; 677 PetscFunctionReturn(0); 678 } 679 680 /*@ 681 MatCompositeSetMergeType - Sets order of `MatCompositeMerge()`. 682 683 Logically Collective on mat 684 685 Input Parameters: 686 + mat - the composite matrix 687 - type - `MAT_COMPOSITE_MERGE RIGHT` (default) to start merge from right with the first added matrix (mat[0]), 688 `MAT_COMPOSITE_MERGE_LEFT` to start merge from left with the last added matrix (mat[nmat-1]) 689 690 Level: advanced 691 692 Note: 693 The resulting matrix is the same regardles of the `MatCompositeMergeType`. Only the order of operation is changed. 694 If set to `MAT_COMPOSITE_MERGE_RIGHT` the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0]))) 695 otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0]. 696 697 .seealso: `MatCreateComposite()`, `MatCompositeMerge()`, `MATCOMPOSITE` 698 @*/ 699 PetscErrorCode MatCompositeSetMergeType(Mat mat, MatCompositeMergeType type) { 700 PetscFunctionBegin; 701 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 702 PetscValidLogicalCollectiveEnum(mat, type, 2); 703 PetscUseMethod(mat, "MatCompositeSetMergeType_C", (Mat, MatCompositeMergeType), (mat, type)); 704 PetscFunctionReturn(0); 705 } 706 707 static PetscErrorCode MatCompositeMerge_Composite(Mat mat) { 708 Mat_Composite *shell = (Mat_Composite *)mat->data; 709 Mat_CompositeLink next = shell->head, prev = shell->tail; 710 Mat tmat, newmat; 711 Vec left, right; 712 PetscScalar scale; 713 PetscInt i; 714 715 PetscFunctionBegin; 716 PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()"); 717 scale = shell->scale; 718 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 719 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 720 i = 0; 721 PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat)); 722 if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i++])); 723 while ((next = next->next)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i++] : 1.0), next->mat, shell->structure)); 724 } else { 725 i = shell->nmat - 1; 726 PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat)); 727 if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i--])); 728 while ((prev = prev->prev)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i--] : 1.0), prev->mat, shell->structure)); 729 } 730 } else { 731 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 732 PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat)); 733 while ((next = next->next)) { 734 PetscCall(MatMatMult(next->mat, tmat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat)); 735 PetscCall(MatDestroy(&tmat)); 736 tmat = newmat; 737 } 738 } else { 739 PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat)); 740 while ((prev = prev->prev)) { 741 PetscCall(MatMatMult(tmat, prev->mat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat)); 742 PetscCall(MatDestroy(&tmat)); 743 tmat = newmat; 744 } 745 } 746 if (shell->scalings) { 747 for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i]; 748 } 749 } 750 751 if ((left = shell->left)) PetscCall(PetscObjectReference((PetscObject)left)); 752 if ((right = shell->right)) PetscCall(PetscObjectReference((PetscObject)right)); 753 754 PetscCall(MatHeaderReplace(mat, &tmat)); 755 756 PetscCall(MatDiagonalScale(mat, left, right)); 757 PetscCall(MatScale(mat, scale)); 758 PetscCall(VecDestroy(&left)); 759 PetscCall(VecDestroy(&right)); 760 PetscFunctionReturn(0); 761 } 762 763 /*@ 764 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 765 by summing or computing the product of all the matrices inside the composite matrix. 766 767 Collective 768 769 Input Parameter: 770 . mat - the composite matrix 771 772 Options Database Keys: 773 + -mat_composite_merge - merge in `MatAssemblyEnd()` 774 - -mat_composite_merge_type - set merge direction 775 776 Level: advanced 777 778 Note: 779 The `MatType` of the resulting matrix will be the same as the `MatType` of the FIRST matrix in the composite matrix. 780 781 .seealso: `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MatCompositeSetMergeType()`, `MATCOMPOSITE` 782 @*/ 783 PetscErrorCode MatCompositeMerge(Mat mat) { 784 PetscFunctionBegin; 785 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 786 PetscUseMethod(mat, "MatCompositeMerge_C", (Mat), (mat)); 787 PetscFunctionReturn(0); 788 } 789 790 static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat, PetscInt *nmat) { 791 Mat_Composite *shell = (Mat_Composite *)mat->data; 792 793 PetscFunctionBegin; 794 *nmat = shell->nmat; 795 PetscFunctionReturn(0); 796 } 797 798 /*@ 799 MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix. 800 801 Not Collective 802 803 Input Parameter: 804 . mat - the composite matrix 805 806 Output Parameter: 807 . nmat - number of matrices in the composite matrix 808 809 Level: advanced 810 811 .seealso: `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE` 812 @*/ 813 PetscErrorCode MatCompositeGetNumberMat(Mat mat, PetscInt *nmat) { 814 PetscFunctionBegin; 815 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 816 PetscValidIntPointer(nmat, 2); 817 PetscUseMethod(mat, "MatCompositeGetNumberMat_C", (Mat, PetscInt *), (mat, nmat)); 818 PetscFunctionReturn(0); 819 } 820 821 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat, PetscInt i, Mat *Ai) { 822 Mat_Composite *shell = (Mat_Composite *)mat->data; 823 Mat_CompositeLink ilink; 824 PetscInt k; 825 826 PetscFunctionBegin; 827 PetscCheck(i < shell->nmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_OUTOFRANGE, "index out of range: %" PetscInt_FMT " >= %" PetscInt_FMT, i, shell->nmat); 828 ilink = shell->head; 829 for (k = 0; k < i; k++) ilink = ilink->next; 830 *Ai = ilink->mat; 831 PetscFunctionReturn(0); 832 } 833 834 /*@ 835 MatCompositeGetMat - Returns the ith matrix from the composite matrix. 836 837 Logically Collective on mat 838 839 Input Parameters: 840 + mat - the composite matrix 841 - i - the number of requested matrix 842 843 Output Parameter: 844 . Ai - ith matrix in composite 845 846 Level: advanced 847 848 .seealso: `MatCreateComposite()`, `MatCompositeGetNumberMat()`, `MatCompositeAddMat()`, `MATCOMPOSITE` 849 @*/ 850 PetscErrorCode MatCompositeGetMat(Mat mat, PetscInt i, Mat *Ai) { 851 PetscFunctionBegin; 852 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 853 PetscValidLogicalCollectiveInt(mat, i, 2); 854 PetscValidPointer(Ai, 3); 855 PetscUseMethod(mat, "MatCompositeGetMat_C", (Mat, PetscInt, Mat *), (mat, i, Ai)); 856 PetscFunctionReturn(0); 857 } 858 859 PetscErrorCode MatCompositeSetScalings_Composite(Mat mat, const PetscScalar *scalings) { 860 Mat_Composite *shell = (Mat_Composite *)mat->data; 861 PetscInt nmat; 862 863 PetscFunctionBegin; 864 PetscCall(MatCompositeGetNumberMat(mat, &nmat)); 865 if (!shell->scalings) PetscCall(PetscMalloc1(nmat, &shell->scalings)); 866 PetscCall(PetscArraycpy(shell->scalings, scalings, nmat)); 867 PetscFunctionReturn(0); 868 } 869 870 /*@ 871 MatCompositeSetScalings - Sets separate scaling factors for component matrices. 872 873 Logically Collective on mat 874 875 Input Parameters: 876 + mat - the composite matrix 877 - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat) 878 879 Level: advanced 880 881 .seealso: `MatScale()`, `MatDiagonalScale()`, `MATCOMPOSITE` 882 @*/ 883 PetscErrorCode MatCompositeSetScalings(Mat mat, const PetscScalar *scalings) { 884 PetscFunctionBegin; 885 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 886 PetscValidScalarPointer(scalings, 2); 887 PetscValidLogicalCollectiveScalar(mat, *scalings, 2); 888 PetscUseMethod(mat, "MatCompositeSetScalings_C", (Mat, const PetscScalar *), (mat, scalings)); 889 PetscFunctionReturn(0); 890 } 891 892 static struct _MatOps MatOps_Values = {NULL, 893 NULL, 894 NULL, 895 MatMult_Composite, 896 MatMultAdd_Composite, 897 /* 5*/ MatMultTranspose_Composite, 898 MatMultTransposeAdd_Composite, 899 NULL, 900 NULL, 901 NULL, 902 /* 10*/ NULL, 903 NULL, 904 NULL, 905 NULL, 906 NULL, 907 /* 15*/ NULL, 908 NULL, 909 MatGetDiagonal_Composite, 910 MatDiagonalScale_Composite, 911 NULL, 912 /* 20*/ NULL, 913 MatAssemblyEnd_Composite, 914 NULL, 915 NULL, 916 /* 24*/ NULL, 917 NULL, 918 NULL, 919 NULL, 920 NULL, 921 /* 29*/ NULL, 922 NULL, 923 NULL, 924 NULL, 925 NULL, 926 /* 34*/ NULL, 927 NULL, 928 NULL, 929 NULL, 930 NULL, 931 /* 39*/ NULL, 932 NULL, 933 NULL, 934 NULL, 935 NULL, 936 /* 44*/ NULL, 937 MatScale_Composite, 938 MatShift_Basic, 939 NULL, 940 NULL, 941 /* 49*/ NULL, 942 NULL, 943 NULL, 944 NULL, 945 NULL, 946 /* 54*/ NULL, 947 NULL, 948 NULL, 949 NULL, 950 NULL, 951 /* 59*/ NULL, 952 MatDestroy_Composite, 953 NULL, 954 NULL, 955 NULL, 956 /* 64*/ NULL, 957 NULL, 958 NULL, 959 NULL, 960 NULL, 961 /* 69*/ NULL, 962 NULL, 963 NULL, 964 NULL, 965 NULL, 966 /* 74*/ NULL, 967 NULL, 968 MatSetFromOptions_Composite, 969 NULL, 970 NULL, 971 /* 79*/ NULL, 972 NULL, 973 NULL, 974 NULL, 975 NULL, 976 /* 84*/ NULL, 977 NULL, 978 NULL, 979 NULL, 980 NULL, 981 /* 89*/ NULL, 982 NULL, 983 NULL, 984 NULL, 985 NULL, 986 /* 94*/ NULL, 987 NULL, 988 NULL, 989 NULL, 990 NULL, 991 /*99*/ NULL, 992 NULL, 993 NULL, 994 NULL, 995 NULL, 996 /*104*/ NULL, 997 NULL, 998 NULL, 999 NULL, 1000 NULL, 1001 /*109*/ NULL, 1002 NULL, 1003 NULL, 1004 NULL, 1005 NULL, 1006 /*114*/ NULL, 1007 NULL, 1008 NULL, 1009 NULL, 1010 NULL, 1011 /*119*/ NULL, 1012 NULL, 1013 NULL, 1014 NULL, 1015 NULL, 1016 /*124*/ NULL, 1017 NULL, 1018 NULL, 1019 NULL, 1020 NULL, 1021 /*129*/ NULL, 1022 NULL, 1023 NULL, 1024 NULL, 1025 NULL, 1026 /*134*/ NULL, 1027 NULL, 1028 NULL, 1029 NULL, 1030 NULL, 1031 /*139*/ NULL, 1032 NULL, 1033 NULL, 1034 NULL, 1035 NULL, 1036 /*144*/ NULL, 1037 NULL, 1038 NULL, 1039 NULL, 1040 NULL, 1041 NULL, 1042 /*150*/ NULL}; 1043 1044 /*MC 1045 MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices. 1046 The matrices need to have a correct size and parallel layout for the sum or product to be valid. 1047 1048 Note: 1049 To use the product of the matrices call `MatCompositeSetType`(mat,`MAT_COMPOSITE_MULTIPLICATIVE`); 1050 1051 Level: advanced 1052 1053 .seealso: `MatCreateComposite()`, `MatCompositeSetScalings()`, `MatCompositeAddMat()`, `MatSetType()`, `MatCompositeSetType()`, `MatCompositeGetType()`, 1054 `MatCompositeSetMatStructure()`, `MatCompositeGetMatStructure()`, `MatCompositeMerge()`, `MatCompositeSetMergeType()`, `MatCompositeGetNumberMat()`, `MatCompositeGetMat()` 1055 M*/ 1056 1057 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A) { 1058 Mat_Composite *b; 1059 1060 PetscFunctionBegin; 1061 PetscCall(PetscNew(&b)); 1062 A->data = (void *)b; 1063 PetscCall(PetscMemcpy(A->ops, &MatOps_Values, sizeof(struct _MatOps))); 1064 1065 PetscCall(PetscLayoutSetUp(A->rmap)); 1066 PetscCall(PetscLayoutSetUp(A->cmap)); 1067 1068 A->assembled = PETSC_TRUE; 1069 A->preallocated = PETSC_TRUE; 1070 b->type = MAT_COMPOSITE_ADDITIVE; 1071 b->scale = 1.0; 1072 b->nmat = 0; 1073 b->merge = PETSC_FALSE; 1074 b->mergetype = MAT_COMPOSITE_MERGE_RIGHT; 1075 b->structure = DIFFERENT_NONZERO_PATTERN; 1076 b->merge_mvctx = PETSC_TRUE; 1077 1078 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeAddMat_C", MatCompositeAddMat_Composite)); 1079 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetType_C", MatCompositeSetType_Composite)); 1080 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetType_C", MatCompositeGetType_Composite)); 1081 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMergeType_C", MatCompositeSetMergeType_Composite)); 1082 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMatStructure_C", MatCompositeSetMatStructure_Composite)); 1083 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMatStructure_C", MatCompositeGetMatStructure_Composite)); 1084 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeMerge_C", MatCompositeMerge_Composite)); 1085 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetNumberMat_C", MatCompositeGetNumberMat_Composite)); 1086 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMat_C", MatCompositeGetMat_Composite)); 1087 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetScalings_C", MatCompositeSetScalings_Composite)); 1088 1089 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCOMPOSITE)); 1090 PetscFunctionReturn(0); 1091 } 1092