1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct _Mat_CompositeLink *Mat_CompositeLink; 5 struct _Mat_CompositeLink { 6 Mat mat; 7 Vec work; 8 Mat_CompositeLink next,prev; 9 }; 10 11 typedef struct { 12 MatCompositeType type; 13 Mat_CompositeLink head,tail; 14 Vec work; 15 PetscScalar scale; /* scale factor supplied with MatScale() */ 16 Vec left,right; /* left and right diagonal scaling provided with MatDiagonalScale() */ 17 Vec leftwork,rightwork; 18 PetscInt nmat; 19 } Mat_Composite; 20 21 PetscErrorCode MatDestroy_Composite(Mat mat) 22 { 23 PetscErrorCode ierr; 24 Mat_Composite *shell = (Mat_Composite*)mat->data; 25 Mat_CompositeLink next = shell->head,oldnext; 26 27 PetscFunctionBegin; 28 while (next) { 29 ierr = MatDestroy(&next->mat);CHKERRQ(ierr); 30 if (next->work && (!next->next || next->work != next->next->work)) { 31 ierr = VecDestroy(&next->work);CHKERRQ(ierr); 32 } 33 oldnext = next; 34 next = next->next; 35 ierr = PetscFree(oldnext);CHKERRQ(ierr); 36 } 37 ierr = VecDestroy(&shell->work);CHKERRQ(ierr); 38 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 39 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 40 ierr = VecDestroy(&shell->leftwork);CHKERRQ(ierr); 41 ierr = VecDestroy(&shell->rightwork);CHKERRQ(ierr); 42 ierr = PetscFree(mat->data);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y) 47 { 48 Mat_Composite *shell = (Mat_Composite*)A->data; 49 Mat_CompositeLink next = shell->head; 50 PetscErrorCode ierr; 51 Vec in,out; 52 53 PetscFunctionBegin; 54 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 55 in = x; 56 if (shell->right) { 57 if (!shell->rightwork) { 58 ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr); 59 } 60 ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr); 61 in = shell->rightwork; 62 } 63 while (next->next) { 64 if (!next->work) { /* should reuse previous work if the same size */ 65 ierr = MatCreateVecs(next->mat,NULL,&next->work);CHKERRQ(ierr); 66 } 67 out = next->work; 68 ierr = MatMult(next->mat,in,out);CHKERRQ(ierr); 69 in = out; 70 next = next->next; 71 } 72 ierr = MatMult(next->mat,in,y);CHKERRQ(ierr); 73 if (shell->left) { 74 ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr); 75 } 76 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y) 81 { 82 Mat_Composite *shell = (Mat_Composite*)A->data; 83 Mat_CompositeLink tail = shell->tail; 84 PetscErrorCode ierr; 85 Vec in,out; 86 87 PetscFunctionBegin; 88 if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 89 in = x; 90 if (shell->left) { 91 if (!shell->leftwork) { 92 ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr); 93 } 94 ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr); 95 in = shell->leftwork; 96 } 97 while (tail->prev) { 98 if (!tail->prev->work) { /* should reuse previous work if the same size */ 99 ierr = MatCreateVecs(tail->mat,NULL,&tail->prev->work);CHKERRQ(ierr); 100 } 101 out = tail->prev->work; 102 ierr = MatMultTranspose(tail->mat,in,out);CHKERRQ(ierr); 103 in = out; 104 tail = tail->prev; 105 } 106 ierr = MatMultTranspose(tail->mat,in,y);CHKERRQ(ierr); 107 if (shell->right) { 108 ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr); 109 } 110 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 114 PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y) 115 { 116 Mat_Composite *shell = (Mat_Composite*)A->data; 117 Mat_CompositeLink next = shell->head; 118 PetscErrorCode ierr; 119 Vec in; 120 121 PetscFunctionBegin; 122 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 123 in = x; 124 if (shell->right) { 125 if (!shell->rightwork) { 126 ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr); 127 } 128 ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr); 129 in = shell->rightwork; 130 } 131 ierr = MatMult(next->mat,in,y);CHKERRQ(ierr); 132 while ((next = next->next)) { 133 ierr = MatMultAdd(next->mat,in,y,y);CHKERRQ(ierr); 134 } 135 if (shell->left) { 136 ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr); 137 } 138 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 142 PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y) 143 { 144 Mat_Composite *shell = (Mat_Composite*)A->data; 145 Mat_CompositeLink next = shell->head; 146 PetscErrorCode ierr; 147 Vec in; 148 149 PetscFunctionBegin; 150 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 151 in = x; 152 if (shell->left) { 153 if (!shell->leftwork) { 154 ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr); 155 } 156 ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr); 157 in = shell->leftwork; 158 } 159 ierr = MatMultTranspose(next->mat,in,y);CHKERRQ(ierr); 160 while ((next = next->next)) { 161 ierr = MatMultTransposeAdd(next->mat,in,y,y);CHKERRQ(ierr); 162 } 163 if (shell->right) { 164 ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr); 165 } 166 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 167 PetscFunctionReturn(0); 168 } 169 170 PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v) 171 { 172 Mat_Composite *shell = (Mat_Composite*)A->data; 173 Mat_CompositeLink next = shell->head; 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 178 if (shell->right || shell->left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get diagonal if left or right scaling"); 179 180 ierr = MatGetDiagonal(next->mat,v);CHKERRQ(ierr); 181 if (next->next && !shell->work) { 182 ierr = VecDuplicate(v,&shell->work);CHKERRQ(ierr); 183 } 184 while ((next = next->next)) { 185 ierr = MatGetDiagonal(next->mat,shell->work);CHKERRQ(ierr); 186 ierr = VecAXPY(v,1.0,shell->work);CHKERRQ(ierr); 187 } 188 ierr = VecScale(v,shell->scale);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191 192 PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t) 193 { 194 PetscErrorCode ierr; 195 PetscBool flg = PETSC_FALSE; 196 197 PetscFunctionBegin; 198 ierr = PetscOptionsGetBool(((PetscObject)Y)->options,((PetscObject)Y)->prefix,"-mat_composite_merge",&flg,NULL);CHKERRQ(ierr); 199 if (flg) { 200 ierr = MatCompositeMerge(Y);CHKERRQ(ierr); 201 } 202 PetscFunctionReturn(0); 203 } 204 205 PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha) 206 { 207 Mat_Composite *a = (Mat_Composite*)inA->data; 208 209 PetscFunctionBegin; 210 a->scale *= alpha; 211 PetscFunctionReturn(0); 212 } 213 214 PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right) 215 { 216 Mat_Composite *a = (Mat_Composite*)inA->data; 217 PetscErrorCode ierr; 218 219 PetscFunctionBegin; 220 if (left) { 221 if (!a->left) { 222 ierr = VecDuplicate(left,&a->left);CHKERRQ(ierr); 223 ierr = VecCopy(left,a->left);CHKERRQ(ierr); 224 } else { 225 ierr = VecPointwiseMult(a->left,left,a->left);CHKERRQ(ierr); 226 } 227 } 228 if (right) { 229 if (!a->right) { 230 ierr = VecDuplicate(right,&a->right);CHKERRQ(ierr); 231 ierr = VecCopy(right,a->right);CHKERRQ(ierr); 232 } else { 233 ierr = VecPointwiseMult(a->right,right,a->right);CHKERRQ(ierr); 234 } 235 } 236 PetscFunctionReturn(0); 237 } 238 239 /*@ 240 MatCreateComposite - Creates a matrix as the sum of zero or more matrices 241 242 Collective on MPI_Comm 243 244 Input Parameters: 245 + comm - MPI communicator 246 . nmat - number of matrices to put in 247 - mats - the matrices 248 249 Output Parameter: 250 . mat - the matrix 251 252 Level: advanced 253 254 Notes: 255 Alternative construction 256 $ MatCreate(comm,&mat); 257 $ MatSetSizes(mat,m,n,M,N); 258 $ MatSetType(mat,MATCOMPOSITE); 259 $ MatCompositeAddMat(mat,mats[0]); 260 $ .... 261 $ MatCompositeAddMat(mat,mats[nmat-1]); 262 $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 263 $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); 264 265 For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0] 266 267 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType 268 269 @*/ 270 PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat) 271 { 272 PetscErrorCode ierr; 273 PetscInt m,n,M,N,i; 274 275 PetscFunctionBegin; 276 if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix"); 277 PetscValidPointer(mat,3); 278 279 ierr = MatGetLocalSize(mats[0],PETSC_IGNORE,&n);CHKERRQ(ierr); 280 ierr = MatGetLocalSize(mats[nmat-1],&m,PETSC_IGNORE);CHKERRQ(ierr); 281 ierr = MatGetSize(mats[0],PETSC_IGNORE,&N);CHKERRQ(ierr); 282 ierr = MatGetSize(mats[nmat-1],&M,PETSC_IGNORE);CHKERRQ(ierr); 283 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 284 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 285 ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr); 286 for (i=0; i<nmat; i++) { 287 ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr); 288 } 289 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 290 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 295 static PetscErrorCode MatCompositeAddMat_Composite(Mat mat,Mat smat) 296 { 297 Mat_Composite *shell = (Mat_Composite*)mat->data; 298 Mat_CompositeLink ilink,next = shell->head; 299 PetscErrorCode ierr; 300 301 PetscFunctionBegin; 302 ierr = PetscNewLog(mat,&ilink);CHKERRQ(ierr); 303 ilink->next = 0; 304 ierr = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr); 305 ilink->mat = smat; 306 307 if (!next) shell->head = ilink; 308 else { 309 while (next->next) { 310 next = next->next; 311 } 312 next->next = ilink; 313 ilink->prev = next; 314 } 315 shell->tail = ilink; 316 shell->nmat += 1; 317 PetscFunctionReturn(0); 318 } 319 320 /*@ 321 MatCompositeAddMat - add another matrix to a composite matrix 322 323 Collective on Mat 324 325 Input Parameters: 326 + mat - the composite matrix 327 - smat - the partial matrix 328 329 Level: advanced 330 331 .seealso: MatCreateComposite() 332 @*/ 333 PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat) 334 { 335 PetscErrorCode ierr; 336 337 PetscFunctionBegin; 338 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 339 PetscValidHeaderSpecific(smat,MAT_CLASSID,2); 340 ierr = PetscUseMethod(mat,"MatCompositeAddMat_C",(Mat,Mat),(mat,smat));CHKERRQ(ierr); 341 PetscFunctionReturn(0); 342 } 343 344 static PetscErrorCode MatCompositeSetType_Composite(Mat mat,MatCompositeType type) 345 { 346 Mat_Composite *b = (Mat_Composite*)mat->data; 347 348 PetscFunctionBegin; 349 if (type == MAT_COMPOSITE_MULTIPLICATIVE) { 350 mat->ops->getdiagonal = 0; 351 mat->ops->mult = MatMult_Composite_Multiplicative; 352 mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative; 353 b->type = MAT_COMPOSITE_MULTIPLICATIVE; 354 } else { 355 mat->ops->getdiagonal = MatGetDiagonal_Composite; 356 mat->ops->mult = MatMult_Composite; 357 mat->ops->multtranspose = MatMultTranspose_Composite; 358 b->type = MAT_COMPOSITE_ADDITIVE; 359 } 360 PetscFunctionReturn(0); 361 } 362 363 /*@ 364 MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product 365 366 Collective on MPI_Comm 367 368 Input Parameters: 369 . mat - the composite matrix 370 371 372 Level: advanced 373 374 Notes: 375 The MatType of the resulting matrix will be the same as the MatType of the FIRST 376 matrix in the composite matrix. 377 378 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE 379 380 @*/ 381 PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type) 382 { 383 PetscErrorCode ierr; 384 385 PetscFunctionBegin; 386 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 387 ierr = PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));CHKERRQ(ierr); 388 PetscFunctionReturn(0); 389 } 390 391 static PetscErrorCode MatCompositeMerge_Composite(Mat mat) 392 { 393 Mat_Composite *shell = (Mat_Composite*)mat->data; 394 Mat_CompositeLink next = shell->head, prev = shell->tail; 395 PetscErrorCode ierr; 396 Mat tmat,newmat; 397 Vec left,right; 398 PetscScalar scale; 399 400 PetscFunctionBegin; 401 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 402 403 PetscFunctionBegin; 404 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 405 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 406 while ((next = next->next)) { 407 ierr = MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 408 } 409 } else { 410 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 411 while ((next = next->next)) { 412 ierr = MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 413 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 414 tmat = newmat; 415 } 416 } 417 418 scale = shell->scale; 419 if ((left = shell->left)) {ierr = PetscObjectReference((PetscObject)left);CHKERRQ(ierr);} 420 if ((right = shell->right)) {ierr = PetscObjectReference((PetscObject)right);CHKERRQ(ierr);} 421 422 ierr = MatHeaderReplace(mat,&tmat);CHKERRQ(ierr); 423 424 ierr = MatDiagonalScale(mat,left,right);CHKERRQ(ierr); 425 ierr = MatScale(mat,scale);CHKERRQ(ierr); 426 ierr = VecDestroy(&left);CHKERRQ(ierr); 427 ierr = VecDestroy(&right);CHKERRQ(ierr); 428 PetscFunctionReturn(0); 429 } 430 431 /*@ 432 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 433 by summing all the matrices inside the composite matrix. 434 435 Collective on MPI_Comm 436 437 Input Parameters: 438 . mat - the composite matrix 439 440 441 Options Database: 442 . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked) 443 444 Level: advanced 445 446 Notes: 447 The MatType of the resulting matrix will be the same as the MatType of the FIRST 448 matrix in the composite matrix. 449 450 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE 451 452 @*/ 453 PetscErrorCode MatCompositeMerge(Mat mat) 454 { 455 PetscErrorCode ierr; 456 457 PetscFunctionBegin; 458 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 459 ierr = PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));CHKERRQ(ierr); 460 PetscFunctionReturn(0); 461 } 462 463 static PetscErrorCode MatCompositeGetNMat_Composite(Mat mat,PetscInt *nmat) 464 { 465 Mat_Composite *shell = (Mat_Composite*)mat->data; 466 467 PetscFunctionBegin; 468 *nmat = shell->nmat; 469 PetscFunctionReturn(0); 470 } 471 472 /*@ 473 MatCompositeGetNMat - Returns the number of matrices in composite. 474 475 Not Collective 476 477 Input Parameter: 478 . mat - the composite matrix 479 480 Output Parameter: 481 . size - the local size 482 . nmat - number of matrices in composite 483 484 Level: advanced 485 486 .seealso: MatCreateComposite(), MatCompositeGetMat() 487 488 @*/ 489 PetscErrorCode MatCompositeGetNMat(Mat mat,PetscInt *nmat) 490 { 491 PetscErrorCode ierr; 492 493 PetscFunctionBegin; 494 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 495 PetscValidPointer(nmat,2); 496 ierr = PetscUseMethod(mat,"MatCompositeGetNMat_C",(Mat,PetscInt*),(mat,nmat));CHKERRQ(ierr); 497 PetscFunctionReturn(0); 498 } 499 500 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai) 501 { 502 Mat_Composite *shell = (Mat_Composite*)mat->data; 503 Mat_CompositeLink ilink; 504 PetscInt k; 505 506 PetscFunctionBegin; 507 if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat); 508 ilink = shell->head; 509 for (k=0; k<i; k++) { 510 ilink = ilink->next; 511 } 512 *Ai = ilink->mat; 513 PetscFunctionReturn(0); 514 } 515 516 /*@ 517 MatCompositeGetMat - Returns the ith matrix from composite. 518 519 Logically Collective on Mat 520 521 Input Parameter: 522 + mat - the composite matrix 523 - i - the number of requested matrix 524 525 Output Parameter: 526 . Ai - ith matrix in composite 527 528 Level: advanced 529 530 .seealso: MatCreateComposite(), MatCompositeGetNMat() 531 532 @*/ 533 PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai) 534 { 535 PetscErrorCode ierr; 536 537 PetscFunctionBegin; 538 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 539 PetscValidLogicalCollectiveInt(mat,i,2); 540 PetscValidPointer(Ai,3); 541 ierr = PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));CHKERRQ(ierr); 542 PetscFunctionReturn(0); 543 } 544 545 static struct _MatOps MatOps_Values = {0, 546 0, 547 0, 548 MatMult_Composite, 549 0, 550 /* 5*/ MatMultTranspose_Composite, 551 0, 552 0, 553 0, 554 0, 555 /* 10*/ 0, 556 0, 557 0, 558 0, 559 0, 560 /* 15*/ 0, 561 0, 562 MatGetDiagonal_Composite, 563 MatDiagonalScale_Composite, 564 0, 565 /* 20*/ 0, 566 MatAssemblyEnd_Composite, 567 0, 568 0, 569 /* 24*/ 0, 570 0, 571 0, 572 0, 573 0, 574 /* 29*/ 0, 575 0, 576 0, 577 0, 578 0, 579 /* 34*/ 0, 580 0, 581 0, 582 0, 583 0, 584 /* 39*/ 0, 585 0, 586 0, 587 0, 588 0, 589 /* 44*/ 0, 590 MatScale_Composite, 591 MatShift_Basic, 592 0, 593 0, 594 /* 49*/ 0, 595 0, 596 0, 597 0, 598 0, 599 /* 54*/ 0, 600 0, 601 0, 602 0, 603 0, 604 /* 59*/ 0, 605 MatDestroy_Composite, 606 0, 607 0, 608 0, 609 /* 64*/ 0, 610 0, 611 0, 612 0, 613 0, 614 /* 69*/ 0, 615 0, 616 0, 617 0, 618 0, 619 /* 74*/ 0, 620 0, 621 0, 622 0, 623 0, 624 /* 79*/ 0, 625 0, 626 0, 627 0, 628 0, 629 /* 84*/ 0, 630 0, 631 0, 632 0, 633 0, 634 /* 89*/ 0, 635 0, 636 0, 637 0, 638 0, 639 /* 94*/ 0, 640 0, 641 0, 642 0, 643 0, 644 /*99*/ 0, 645 0, 646 0, 647 0, 648 0, 649 /*104*/ 0, 650 0, 651 0, 652 0, 653 0, 654 /*109*/ 0, 655 0, 656 0, 657 0, 658 0, 659 /*114*/ 0, 660 0, 661 0, 662 0, 663 0, 664 /*119*/ 0, 665 0, 666 0, 667 0, 668 0, 669 /*124*/ 0, 670 0, 671 0, 672 0, 673 0, 674 /*129*/ 0, 675 0, 676 0, 677 0, 678 0, 679 /*134*/ 0, 680 0, 681 0, 682 0, 683 0, 684 /*139*/ 0, 685 0, 686 0 687 }; 688 689 /*MC 690 MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices. 691 The matrices need to have a correct size and parallel layout for the sum or product to be valid. 692 693 Notes: 694 to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE); 695 696 Level: advanced 697 698 .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType 699 M*/ 700 701 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A) 702 { 703 Mat_Composite *b; 704 PetscErrorCode ierr; 705 706 PetscFunctionBegin; 707 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 708 A->data = (void*)b; 709 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 710 711 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 712 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 713 714 A->assembled = PETSC_TRUE; 715 A->preallocated = PETSC_TRUE; 716 b->type = MAT_COMPOSITE_ADDITIVE; 717 b->scale = 1.0; 718 b->nmat = 0; 719 720 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);CHKERRQ(ierr); 721 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);CHKERRQ(ierr); 722 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);CHKERRQ(ierr); 723 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNMat_C",MatCompositeGetNMat_Composite);CHKERRQ(ierr); 724 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);CHKERRQ(ierr); 725 726 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 727 PetscFunctionReturn(0); 728 } 729 730