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