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