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