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