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 Notes: 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 @*/ 474 PetscErrorCode MatCreateComposite(MPI_Comm comm, PetscInt nmat, const Mat *mats, Mat *mat) { 475 PetscInt m, n, M, N, i; 476 477 PetscFunctionBegin; 478 PetscCheck(nmat >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must pass in at least one matrix"); 479 PetscValidPointer(mat, 4); 480 481 PetscCall(MatGetLocalSize(mats[0], PETSC_IGNORE, &n)); 482 PetscCall(MatGetLocalSize(mats[nmat - 1], &m, PETSC_IGNORE)); 483 PetscCall(MatGetSize(mats[0], PETSC_IGNORE, &N)); 484 PetscCall(MatGetSize(mats[nmat - 1], &M, PETSC_IGNORE)); 485 PetscCall(MatCreate(comm, mat)); 486 PetscCall(MatSetSizes(*mat, m, n, M, N)); 487 PetscCall(MatSetType(*mat, MATCOMPOSITE)); 488 for (i = 0; i < nmat; i++) PetscCall(MatCompositeAddMat(*mat, mats[i])); 489 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 490 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 491 PetscFunctionReturn(0); 492 } 493 494 static PetscErrorCode MatCompositeAddMat_Composite(Mat mat, Mat smat) { 495 Mat_Composite *shell = (Mat_Composite *)mat->data; 496 Mat_CompositeLink ilink, next = shell->head; 497 498 PetscFunctionBegin; 499 PetscCall(PetscNewLog(mat, &ilink)); 500 ilink->next = NULL; 501 PetscCall(PetscObjectReference((PetscObject)smat)); 502 ilink->mat = smat; 503 504 if (!next) shell->head = ilink; 505 else { 506 while (next->next) { next = next->next; } 507 next->next = ilink; 508 ilink->prev = next; 509 } 510 shell->tail = ilink; 511 shell->nmat += 1; 512 513 /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */ 514 if (shell->scalings) { 515 PetscCall(PetscRealloc(sizeof(PetscScalar) * shell->nmat, &shell->scalings)); 516 shell->scalings[shell->nmat - 1] = 1.0; 517 } 518 PetscFunctionReturn(0); 519 } 520 521 /*@ 522 MatCompositeAddMat - Add another matrix to a composite matrix. 523 524 Collective on Mat 525 526 Input Parameters: 527 + mat - the composite matrix 528 - smat - the partial matrix 529 530 Level: advanced 531 532 .seealso: `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE` 533 @*/ 534 PetscErrorCode MatCompositeAddMat(Mat mat, Mat smat) { 535 PetscFunctionBegin; 536 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 537 PetscValidHeaderSpecific(smat, MAT_CLASSID, 2); 538 PetscUseMethod(mat, "MatCompositeAddMat_C", (Mat, Mat), (mat, smat)); 539 PetscFunctionReturn(0); 540 } 541 542 static PetscErrorCode MatCompositeSetType_Composite(Mat mat, MatCompositeType type) { 543 Mat_Composite *b = (Mat_Composite *)mat->data; 544 545 PetscFunctionBegin; 546 b->type = type; 547 if (type == MAT_COMPOSITE_MULTIPLICATIVE) { 548 mat->ops->getdiagonal = NULL; 549 mat->ops->mult = MatMult_Composite_Multiplicative; 550 mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative; 551 b->merge_mvctx = PETSC_FALSE; 552 } else { 553 mat->ops->getdiagonal = MatGetDiagonal_Composite; 554 mat->ops->mult = MatMult_Composite; 555 mat->ops->multtranspose = MatMultTranspose_Composite; 556 } 557 PetscFunctionReturn(0); 558 } 559 560 /*@ 561 MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product. 562 563 Logically Collective on Mat 564 565 Input Parameters: 566 . mat - the composite matrix 567 568 Level: advanced 569 570 .seealso: `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeGetType()`, `MATCOMPOSITE` 571 572 @*/ 573 PetscErrorCode MatCompositeSetType(Mat mat, MatCompositeType type) { 574 PetscFunctionBegin; 575 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 576 PetscValidLogicalCollectiveEnum(mat, type, 2); 577 PetscUseMethod(mat, "MatCompositeSetType_C", (Mat, MatCompositeType), (mat, type)); 578 PetscFunctionReturn(0); 579 } 580 581 static PetscErrorCode MatCompositeGetType_Composite(Mat mat, MatCompositeType *type) { 582 Mat_Composite *b = (Mat_Composite *)mat->data; 583 584 PetscFunctionBegin; 585 *type = b->type; 586 PetscFunctionReturn(0); 587 } 588 589 /*@ 590 MatCompositeGetType - Returns type of composite. 591 592 Not Collective 593 594 Input Parameter: 595 . mat - the composite matrix 596 597 Output Parameter: 598 . type - type of composite 599 600 Level: advanced 601 602 .seealso: `MatCreateComposite()`, `MatCompositeSetType()`, `MATCOMPOSITE` 603 604 @*/ 605 PetscErrorCode MatCompositeGetType(Mat mat, MatCompositeType *type) { 606 PetscFunctionBegin; 607 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 608 PetscValidPointer(type, 2); 609 PetscUseMethod(mat, "MatCompositeGetType_C", (Mat, MatCompositeType *), (mat, type)); 610 PetscFunctionReturn(0); 611 } 612 613 static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat, MatStructure str) { 614 Mat_Composite *b = (Mat_Composite *)mat->data; 615 616 PetscFunctionBegin; 617 b->structure = str; 618 PetscFunctionReturn(0); 619 } 620 621 /*@ 622 MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix. 623 624 Not Collective 625 626 Input Parameters: 627 + mat - the composite matrix 628 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN (default) or SUBSET_NONZERO_PATTERN 629 630 Level: advanced 631 632 Notes: 633 Information about the matrices structure is used in MatCompositeMerge() for additive composite matrix. 634 635 .seealso: `MatAXPY()`, `MatCreateComposite()`, `MatCompositeMerge()` `MatCompositeGetMatStructure()`, `MATCOMPOSITE` 636 637 @*/ 638 PetscErrorCode MatCompositeSetMatStructure(Mat mat, MatStructure str) { 639 PetscFunctionBegin; 640 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 641 PetscUseMethod(mat, "MatCompositeSetMatStructure_C", (Mat, MatStructure), (mat, str)); 642 PetscFunctionReturn(0); 643 } 644 645 static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat, MatStructure *str) { 646 Mat_Composite *b = (Mat_Composite *)mat->data; 647 648 PetscFunctionBegin; 649 *str = b->structure; 650 PetscFunctionReturn(0); 651 } 652 653 /*@ 654 MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix. 655 656 Not Collective 657 658 Input Parameter: 659 . mat - the composite matrix 660 661 Output Parameter: 662 . str - structure of the matrices 663 664 Level: advanced 665 666 .seealso: `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MATCOMPOSITE` 667 668 @*/ 669 PetscErrorCode MatCompositeGetMatStructure(Mat mat, MatStructure *str) { 670 PetscFunctionBegin; 671 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 672 PetscValidPointer(str, 2); 673 PetscUseMethod(mat, "MatCompositeGetMatStructure_C", (Mat, MatStructure *), (mat, str)); 674 PetscFunctionReturn(0); 675 } 676 677 static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat, MatCompositeMergeType type) { 678 Mat_Composite *shell = (Mat_Composite *)mat->data; 679 680 PetscFunctionBegin; 681 shell->mergetype = type; 682 PetscFunctionReturn(0); 683 } 684 685 /*@ 686 MatCompositeSetMergeType - Sets order of MatCompositeMerge(). 687 688 Logically Collective on Mat 689 690 Input Parameters: 691 + mat - the composite matrix 692 - type - MAT_COMPOSITE_MERGE RIGHT (default) to start merge from right with the first added matrix (mat[0]), 693 MAT_COMPOSITE_MERGE_LEFT to start merge from left with the last added matrix (mat[nmat-1]) 694 695 Level: advanced 696 697 Notes: 698 The resulting matrix is the same regardles of the MergeType. Only the order of operation is changed. 699 If set to MAT_COMPOSITE_MERGE_RIGHT the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0]))) 700 otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0]. 701 702 .seealso: `MatCreateComposite()`, `MatCompositeMerge()`, `MATCOMPOSITE` 703 704 @*/ 705 PetscErrorCode MatCompositeSetMergeType(Mat mat, MatCompositeMergeType type) { 706 PetscFunctionBegin; 707 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 708 PetscValidLogicalCollectiveEnum(mat, type, 2); 709 PetscUseMethod(mat, "MatCompositeSetMergeType_C", (Mat, MatCompositeMergeType), (mat, type)); 710 PetscFunctionReturn(0); 711 } 712 713 static PetscErrorCode MatCompositeMerge_Composite(Mat mat) { 714 Mat_Composite *shell = (Mat_Composite *)mat->data; 715 Mat_CompositeLink next = shell->head, prev = shell->tail; 716 Mat tmat, newmat; 717 Vec left, right; 718 PetscScalar scale; 719 PetscInt i; 720 721 PetscFunctionBegin; 722 PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()"); 723 scale = shell->scale; 724 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 725 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 726 i = 0; 727 PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat)); 728 if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i++])); 729 while ((next = next->next)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i++] : 1.0), next->mat, shell->structure)); 730 } else { 731 i = shell->nmat - 1; 732 PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat)); 733 if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i--])); 734 while ((prev = prev->prev)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i--] : 1.0), prev->mat, shell->structure)); 735 } 736 } else { 737 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 738 PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat)); 739 while ((next = next->next)) { 740 PetscCall(MatMatMult(next->mat, tmat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat)); 741 PetscCall(MatDestroy(&tmat)); 742 tmat = newmat; 743 } 744 } else { 745 PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat)); 746 while ((prev = prev->prev)) { 747 PetscCall(MatMatMult(tmat, prev->mat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat)); 748 PetscCall(MatDestroy(&tmat)); 749 tmat = newmat; 750 } 751 } 752 if (shell->scalings) { 753 for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i]; 754 } 755 } 756 757 if ((left = shell->left)) PetscCall(PetscObjectReference((PetscObject)left)); 758 if ((right = shell->right)) PetscCall(PetscObjectReference((PetscObject)right)); 759 760 PetscCall(MatHeaderReplace(mat, &tmat)); 761 762 PetscCall(MatDiagonalScale(mat, left, right)); 763 PetscCall(MatScale(mat, scale)); 764 PetscCall(VecDestroy(&left)); 765 PetscCall(VecDestroy(&right)); 766 PetscFunctionReturn(0); 767 } 768 769 /*@ 770 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 771 by summing or computing the product of all the matrices inside the composite matrix. 772 773 Collective 774 775 Input Parameter: 776 . mat - the composite matrix 777 778 Options Database Keys: 779 + -mat_composite_merge - merge in MatAssemblyEnd() 780 - -mat_composite_merge_type - set merge direction 781 782 Level: advanced 783 784 Notes: 785 The MatType of the resulting matrix will be the same as the MatType of the FIRST 786 matrix in the composite matrix. 787 788 .seealso: `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MatCompositeSetMergeType()`, `MATCOMPOSITE` 789 790 @*/ 791 PetscErrorCode MatCompositeMerge(Mat mat) { 792 PetscFunctionBegin; 793 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 794 PetscUseMethod(mat, "MatCompositeMerge_C", (Mat), (mat)); 795 PetscFunctionReturn(0); 796 } 797 798 static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat, PetscInt *nmat) { 799 Mat_Composite *shell = (Mat_Composite *)mat->data; 800 801 PetscFunctionBegin; 802 *nmat = shell->nmat; 803 PetscFunctionReturn(0); 804 } 805 806 /*@ 807 MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix. 808 809 Not Collective 810 811 Input Parameter: 812 . mat - the composite matrix 813 814 Output Parameter: 815 . nmat - number of matrices in the composite matrix 816 817 Level: advanced 818 819 .seealso: `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE` 820 821 @*/ 822 PetscErrorCode MatCompositeGetNumberMat(Mat mat, PetscInt *nmat) { 823 PetscFunctionBegin; 824 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 825 PetscValidIntPointer(nmat, 2); 826 PetscUseMethod(mat, "MatCompositeGetNumberMat_C", (Mat, PetscInt *), (mat, nmat)); 827 PetscFunctionReturn(0); 828 } 829 830 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat, PetscInt i, Mat *Ai) { 831 Mat_Composite *shell = (Mat_Composite *)mat->data; 832 Mat_CompositeLink ilink; 833 PetscInt k; 834 835 PetscFunctionBegin; 836 PetscCheck(i < shell->nmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_OUTOFRANGE, "index out of range: %" PetscInt_FMT " >= %" PetscInt_FMT, i, shell->nmat); 837 ilink = shell->head; 838 for (k = 0; k < i; k++) { ilink = ilink->next; } 839 *Ai = ilink->mat; 840 PetscFunctionReturn(0); 841 } 842 843 /*@ 844 MatCompositeGetMat - Returns the ith matrix from the composite matrix. 845 846 Logically Collective on Mat 847 848 Input Parameters: 849 + mat - the composite matrix 850 - i - the number of requested matrix 851 852 Output Parameter: 853 . Ai - ith matrix in composite 854 855 Level: advanced 856 857 .seealso: `MatCreateComposite()`, `MatCompositeGetNumberMat()`, `MatCompositeAddMat()`, `MATCOMPOSITE` 858 859 @*/ 860 PetscErrorCode MatCompositeGetMat(Mat mat, PetscInt i, Mat *Ai) { 861 PetscFunctionBegin; 862 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 863 PetscValidLogicalCollectiveInt(mat, i, 2); 864 PetscValidPointer(Ai, 3); 865 PetscUseMethod(mat, "MatCompositeGetMat_C", (Mat, PetscInt, Mat *), (mat, i, Ai)); 866 PetscFunctionReturn(0); 867 } 868 869 PetscErrorCode MatCompositeSetScalings_Composite(Mat mat, const PetscScalar *scalings) { 870 Mat_Composite *shell = (Mat_Composite *)mat->data; 871 PetscInt nmat; 872 873 PetscFunctionBegin; 874 PetscCall(MatCompositeGetNumberMat(mat, &nmat)); 875 if (!shell->scalings) PetscCall(PetscMalloc1(nmat, &shell->scalings)); 876 PetscCall(PetscArraycpy(shell->scalings, scalings, nmat)); 877 PetscFunctionReturn(0); 878 } 879 880 /*@ 881 MatCompositeSetScalings - Sets separate scaling factors for component matrices. 882 883 Logically Collective on Mat 884 885 Input Parameters: 886 + mat - the composite matrix 887 - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat) 888 889 Level: advanced 890 891 .seealso: `MatScale()`, `MatDiagonalScale()`, `MATCOMPOSITE` 892 893 @*/ 894 PetscErrorCode MatCompositeSetScalings(Mat mat, const PetscScalar *scalings) { 895 PetscFunctionBegin; 896 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 897 PetscValidScalarPointer(scalings, 2); 898 PetscValidLogicalCollectiveScalar(mat, *scalings, 2); 899 PetscUseMethod(mat, "MatCompositeSetScalings_C", (Mat, const PetscScalar *), (mat, scalings)); 900 PetscFunctionReturn(0); 901 } 902 903 static struct _MatOps MatOps_Values = {NULL, 904 NULL, 905 NULL, 906 MatMult_Composite, 907 MatMultAdd_Composite, 908 /* 5*/ MatMultTranspose_Composite, 909 MatMultTransposeAdd_Composite, 910 NULL, 911 NULL, 912 NULL, 913 /* 10*/ NULL, 914 NULL, 915 NULL, 916 NULL, 917 NULL, 918 /* 15*/ NULL, 919 NULL, 920 MatGetDiagonal_Composite, 921 MatDiagonalScale_Composite, 922 NULL, 923 /* 20*/ NULL, 924 MatAssemblyEnd_Composite, 925 NULL, 926 NULL, 927 /* 24*/ NULL, 928 NULL, 929 NULL, 930 NULL, 931 NULL, 932 /* 29*/ NULL, 933 NULL, 934 NULL, 935 NULL, 936 NULL, 937 /* 34*/ NULL, 938 NULL, 939 NULL, 940 NULL, 941 NULL, 942 /* 39*/ NULL, 943 NULL, 944 NULL, 945 NULL, 946 NULL, 947 /* 44*/ NULL, 948 MatScale_Composite, 949 MatShift_Basic, 950 NULL, 951 NULL, 952 /* 49*/ NULL, 953 NULL, 954 NULL, 955 NULL, 956 NULL, 957 /* 54*/ NULL, 958 NULL, 959 NULL, 960 NULL, 961 NULL, 962 /* 59*/ NULL, 963 MatDestroy_Composite, 964 NULL, 965 NULL, 966 NULL, 967 /* 64*/ NULL, 968 NULL, 969 NULL, 970 NULL, 971 NULL, 972 /* 69*/ NULL, 973 NULL, 974 NULL, 975 NULL, 976 NULL, 977 /* 74*/ NULL, 978 NULL, 979 MatSetFromOptions_Composite, 980 NULL, 981 NULL, 982 /* 79*/ NULL, 983 NULL, 984 NULL, 985 NULL, 986 NULL, 987 /* 84*/ NULL, 988 NULL, 989 NULL, 990 NULL, 991 NULL, 992 /* 89*/ NULL, 993 NULL, 994 NULL, 995 NULL, 996 NULL, 997 /* 94*/ NULL, 998 NULL, 999 NULL, 1000 NULL, 1001 NULL, 1002 /*99*/ NULL, 1003 NULL, 1004 NULL, 1005 NULL, 1006 NULL, 1007 /*104*/ NULL, 1008 NULL, 1009 NULL, 1010 NULL, 1011 NULL, 1012 /*109*/ NULL, 1013 NULL, 1014 NULL, 1015 NULL, 1016 NULL, 1017 /*114*/ NULL, 1018 NULL, 1019 NULL, 1020 NULL, 1021 NULL, 1022 /*119*/ NULL, 1023 NULL, 1024 NULL, 1025 NULL, 1026 NULL, 1027 /*124*/ NULL, 1028 NULL, 1029 NULL, 1030 NULL, 1031 NULL, 1032 /*129*/ NULL, 1033 NULL, 1034 NULL, 1035 NULL, 1036 NULL, 1037 /*134*/ NULL, 1038 NULL, 1039 NULL, 1040 NULL, 1041 NULL, 1042 /*139*/ NULL, 1043 NULL, 1044 NULL, 1045 NULL, 1046 NULL, 1047 /*144*/ NULL, 1048 NULL, 1049 NULL, 1050 NULL, 1051 NULL, 1052 NULL, 1053 /*150*/ NULL}; 1054 1055 /*MC 1056 MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices. 1057 The matrices need to have a correct size and parallel layout for the sum or product to be valid. 1058 1059 Notes: 1060 to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE); 1061 1062 Level: advanced 1063 1064 .seealso: `MatCreateComposite()`, `MatCompositeSetScalings()`, `MatCompositeAddMat()`, `MatSetType()`, `MatCompositeSetType()`, `MatCompositeGetType()`, `MatCompositeSetMatStructure()`, `MatCompositeGetMatStructure()`, `MatCompositeMerge()`, `MatCompositeSetMergeType()`, `MatCompositeGetNumberMat()`, `MatCompositeGetMat()` 1065 M*/ 1066 1067 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A) { 1068 Mat_Composite *b; 1069 1070 PetscFunctionBegin; 1071 PetscCall(PetscNewLog(A, &b)); 1072 A->data = (void *)b; 1073 PetscCall(PetscMemcpy(A->ops, &MatOps_Values, sizeof(struct _MatOps))); 1074 1075 PetscCall(PetscLayoutSetUp(A->rmap)); 1076 PetscCall(PetscLayoutSetUp(A->cmap)); 1077 1078 A->assembled = PETSC_TRUE; 1079 A->preallocated = PETSC_TRUE; 1080 b->type = MAT_COMPOSITE_ADDITIVE; 1081 b->scale = 1.0; 1082 b->nmat = 0; 1083 b->merge = PETSC_FALSE; 1084 b->mergetype = MAT_COMPOSITE_MERGE_RIGHT; 1085 b->structure = DIFFERENT_NONZERO_PATTERN; 1086 b->merge_mvctx = PETSC_TRUE; 1087 1088 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeAddMat_C", MatCompositeAddMat_Composite)); 1089 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetType_C", MatCompositeSetType_Composite)); 1090 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetType_C", MatCompositeGetType_Composite)); 1091 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMergeType_C", MatCompositeSetMergeType_Composite)); 1092 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMatStructure_C", MatCompositeSetMatStructure_Composite)); 1093 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMatStructure_C", MatCompositeGetMatStructure_Composite)); 1094 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeMerge_C", MatCompositeMerge_Composite)); 1095 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetNumberMat_C", MatCompositeGetNumberMat_Composite)); 1096 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMat_C", MatCompositeGetMat_Composite)); 1097 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetScalings_C", MatCompositeSetScalings_Composite)); 1098 1099 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCOMPOSITE)); 1100 PetscFunctionReturn(0); 1101 } 1102