1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 const char *const MatCompositeMergeTypes[] = {"left","right","MatCompositeMergeType","MAT_COMPOSITE_",0}; 5 6 typedef struct _Mat_CompositeLink *Mat_CompositeLink; 7 struct _Mat_CompositeLink { 8 Mat mat; 9 Vec work; 10 Mat_CompositeLink next,prev; 11 }; 12 13 typedef struct { 14 MatCompositeType type; 15 Mat_CompositeLink head,tail; 16 Vec work; 17 PetscScalar scale; /* scale factor supplied with MatScale() */ 18 Vec left,right; /* left and right diagonal scaling provided with MatDiagonalScale() */ 19 Vec leftwork,rightwork; 20 PetscInt nmat; 21 PetscBool merge; 22 MatCompositeMergeType mergetype; 23 MatStructure structure; 24 } Mat_Composite; 25 26 PetscErrorCode MatDestroy_Composite(Mat mat) 27 { 28 PetscErrorCode ierr; 29 Mat_Composite *shell = (Mat_Composite*)mat->data; 30 Mat_CompositeLink next = shell->head,oldnext; 31 32 PetscFunctionBegin; 33 while (next) { 34 ierr = MatDestroy(&next->mat);CHKERRQ(ierr); 35 if (next->work && (!next->next || next->work != next->next->work)) { 36 ierr = VecDestroy(&next->work);CHKERRQ(ierr); 37 } 38 oldnext = next; 39 next = next->next; 40 ierr = PetscFree(oldnext);CHKERRQ(ierr); 41 } 42 ierr = VecDestroy(&shell->work);CHKERRQ(ierr); 43 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 44 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 45 ierr = VecDestroy(&shell->leftwork);CHKERRQ(ierr); 46 ierr = VecDestroy(&shell->rightwork);CHKERRQ(ierr); 47 ierr = PetscFree(mat->data);CHKERRQ(ierr); 48 PetscFunctionReturn(0); 49 } 50 51 PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y) 52 { 53 Mat_Composite *shell = (Mat_Composite*)A->data; 54 Mat_CompositeLink next = shell->head; 55 PetscErrorCode ierr; 56 Vec in,out; 57 58 PetscFunctionBegin; 59 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 60 in = x; 61 if (shell->right) { 62 if (!shell->rightwork) { 63 ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr); 64 } 65 ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr); 66 in = shell->rightwork; 67 } 68 while (next->next) { 69 if (!next->work) { /* should reuse previous work if the same size */ 70 ierr = MatCreateVecs(next->mat,NULL,&next->work);CHKERRQ(ierr); 71 } 72 out = next->work; 73 ierr = MatMult(next->mat,in,out);CHKERRQ(ierr); 74 in = out; 75 next = next->next; 76 } 77 ierr = MatMult(next->mat,in,y);CHKERRQ(ierr); 78 if (shell->left) { 79 ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr); 80 } 81 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y) 86 { 87 Mat_Composite *shell = (Mat_Composite*)A->data; 88 Mat_CompositeLink tail = shell->tail; 89 PetscErrorCode ierr; 90 Vec in,out; 91 92 PetscFunctionBegin; 93 if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 94 in = x; 95 if (shell->left) { 96 if (!shell->leftwork) { 97 ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr); 98 } 99 ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr); 100 in = shell->leftwork; 101 } 102 while (tail->prev) { 103 if (!tail->prev->work) { /* should reuse previous work if the same size */ 104 ierr = MatCreateVecs(tail->mat,NULL,&tail->prev->work);CHKERRQ(ierr); 105 } 106 out = tail->prev->work; 107 ierr = MatMultTranspose(tail->mat,in,out);CHKERRQ(ierr); 108 in = out; 109 tail = tail->prev; 110 } 111 ierr = MatMultTranspose(tail->mat,in,y);CHKERRQ(ierr); 112 if (shell->right) { 113 ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr); 114 } 115 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y) 120 { 121 Mat_Composite *shell = (Mat_Composite*)A->data; 122 Mat_CompositeLink next = shell->head; 123 PetscErrorCode ierr; 124 Vec in; 125 126 PetscFunctionBegin; 127 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 128 in = x; 129 if (shell->right) { 130 if (!shell->rightwork) { 131 ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr); 132 } 133 ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr); 134 in = shell->rightwork; 135 } 136 ierr = MatMult(next->mat,in,y);CHKERRQ(ierr); 137 while ((next = next->next)) { 138 ierr = MatMultAdd(next->mat,in,y,y);CHKERRQ(ierr); 139 } 140 if (shell->left) { 141 ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr); 142 } 143 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 144 PetscFunctionReturn(0); 145 } 146 147 PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y) 148 { 149 Mat_Composite *shell = (Mat_Composite*)A->data; 150 Mat_CompositeLink next = shell->head; 151 PetscErrorCode ierr; 152 Vec in; 153 154 PetscFunctionBegin; 155 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 156 in = x; 157 if (shell->left) { 158 if (!shell->leftwork) { 159 ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr); 160 } 161 ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr); 162 in = shell->leftwork; 163 } 164 ierr = MatMultTranspose(next->mat,in,y);CHKERRQ(ierr); 165 while ((next = next->next)) { 166 ierr = MatMultTransposeAdd(next->mat,in,y,y);CHKERRQ(ierr); 167 } 168 if (shell->right) { 169 ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr); 170 } 171 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 PetscErrorCode MatMultAdd_Composite(Mat A,Vec x,Vec y,Vec z) 176 { 177 Mat_Composite *shell = (Mat_Composite*)A->data; 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 if (y != z) { 182 ierr = MatMult(A,x,z);CHKERRQ(ierr); 183 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 184 } else { 185 if (!shell->leftwork) { 186 ierr = VecDuplicate(z,&shell->leftwork);CHKERRQ(ierr); 187 } 188 ierr = MatMult(A,x,shell->leftwork);CHKERRQ(ierr); 189 ierr = VecCopy(y,z);CHKERRQ(ierr); 190 ierr = VecAXPY(z,1.0,shell->leftwork);CHKERRQ(ierr); 191 } 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode MatMultTransposeAdd_Composite(Mat A,Vec x,Vec y, Vec z) 196 { 197 Mat_Composite *shell = (Mat_Composite*)A->data; 198 PetscErrorCode ierr; 199 200 PetscFunctionBegin; 201 if (y != z) { 202 ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr); 203 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 204 } else { 205 if (!shell->rightwork) { 206 ierr = VecDuplicate(z,&shell->rightwork);CHKERRQ(ierr); 207 } 208 ierr = MatMultTranspose(A,x,shell->rightwork);CHKERRQ(ierr); 209 ierr = VecCopy(y,z);CHKERRQ(ierr); 210 ierr = VecAXPY(z,1.0,shell->rightwork);CHKERRQ(ierr); 211 } 212 PetscFunctionReturn(0); 213 } 214 215 PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v) 216 { 217 Mat_Composite *shell = (Mat_Composite*)A->data; 218 Mat_CompositeLink next = shell->head; 219 PetscErrorCode ierr; 220 221 PetscFunctionBegin; 222 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 223 if (shell->right || shell->left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get diagonal if left or right scaling"); 224 225 ierr = MatGetDiagonal(next->mat,v);CHKERRQ(ierr); 226 if (next->next && !shell->work) { 227 ierr = VecDuplicate(v,&shell->work);CHKERRQ(ierr); 228 } 229 while ((next = next->next)) { 230 ierr = MatGetDiagonal(next->mat,shell->work);CHKERRQ(ierr); 231 ierr = VecAXPY(v,1.0,shell->work);CHKERRQ(ierr); 232 } 233 ierr = VecScale(v,shell->scale);CHKERRQ(ierr); 234 PetscFunctionReturn(0); 235 } 236 237 PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t) 238 { 239 Mat_Composite *shell = (Mat_Composite*)Y->data; 240 PetscErrorCode ierr; 241 242 PetscFunctionBegin; 243 if (shell->merge) { 244 ierr = MatCompositeMerge(Y);CHKERRQ(ierr); 245 } 246 PetscFunctionReturn(0); 247 } 248 249 PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha) 250 { 251 Mat_Composite *a = (Mat_Composite*)inA->data; 252 253 PetscFunctionBegin; 254 a->scale *= alpha; 255 PetscFunctionReturn(0); 256 } 257 258 PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right) 259 { 260 Mat_Composite *a = (Mat_Composite*)inA->data; 261 PetscErrorCode ierr; 262 263 PetscFunctionBegin; 264 if (left) { 265 if (!a->left) { 266 ierr = VecDuplicate(left,&a->left);CHKERRQ(ierr); 267 ierr = VecCopy(left,a->left);CHKERRQ(ierr); 268 } else { 269 ierr = VecPointwiseMult(a->left,left,a->left);CHKERRQ(ierr); 270 } 271 } 272 if (right) { 273 if (!a->right) { 274 ierr = VecDuplicate(right,&a->right);CHKERRQ(ierr); 275 ierr = VecCopy(right,a->right);CHKERRQ(ierr); 276 } else { 277 ierr = VecPointwiseMult(a->right,right,a->right);CHKERRQ(ierr); 278 } 279 } 280 PetscFunctionReturn(0); 281 } 282 283 PetscErrorCode MatSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,Mat A) 284 { 285 Mat_Composite *a = (Mat_Composite*)A->data; 286 PetscErrorCode ierr; 287 288 PetscFunctionBegin; 289 ierr = PetscOptionsHead(PetscOptionsObject,"MATCOMPOSITE options");CHKERRQ(ierr); 290 ierr = PetscOptionsBool("-mat_composite_merge","Merge at MatAssemblyEnd","MatCompositeMerge",a->merge,&a->merge,NULL);CHKERRQ(ierr); 291 ierr = PetscOptionsEnum("-mat_composite_merge_type","Set composite merge direction","MatCompositeSetMergeType",MatCompositeMergeTypes,(PetscEnum)a->mergetype,(PetscEnum*)&a->mergetype,NULL);CHKERRQ(ierr); 292 ierr = PetscOptionsTail();CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 /*@ 297 MatCreateComposite - Creates a matrix as the sum or product of one or more matrices 298 299 Collective on MPI_Comm 300 301 Input Parameters: 302 + comm - MPI communicator 303 . nmat - number of matrices to put in 304 - mats - the matrices 305 306 Output Parameter: 307 . mat - the matrix 308 309 Options Database Keys: 310 + -mat_composite_merge - merge in MatAssemblyEnd() 311 - -mat_composite_merge_type - set merge direction 312 313 Level: advanced 314 315 Notes: 316 Alternative construction 317 $ MatCreate(comm,&mat); 318 $ MatSetSizes(mat,m,n,M,N); 319 $ MatSetType(mat,MATCOMPOSITE); 320 $ MatCompositeAddMat(mat,mats[0]); 321 $ .... 322 $ MatCompositeAddMat(mat,mats[nmat-1]); 323 $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 324 $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); 325 326 For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0] 327 328 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeGetMat(), MatCompositeMerge(), MatCompositeSetType(), MATCOMPOSITE 329 330 @*/ 331 PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat) 332 { 333 PetscErrorCode ierr; 334 PetscInt m,n,M,N,i; 335 336 PetscFunctionBegin; 337 if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix"); 338 PetscValidPointer(mat,3); 339 340 ierr = MatGetLocalSize(mats[0],PETSC_IGNORE,&n);CHKERRQ(ierr); 341 ierr = MatGetLocalSize(mats[nmat-1],&m,PETSC_IGNORE);CHKERRQ(ierr); 342 ierr = MatGetSize(mats[0],PETSC_IGNORE,&N);CHKERRQ(ierr); 343 ierr = MatGetSize(mats[nmat-1],&M,PETSC_IGNORE);CHKERRQ(ierr); 344 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 345 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 346 ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr); 347 for (i=0; i<nmat; i++) { 348 ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr); 349 } 350 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 351 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 352 PetscFunctionReturn(0); 353 } 354 355 356 static PetscErrorCode MatCompositeAddMat_Composite(Mat mat,Mat smat) 357 { 358 Mat_Composite *shell = (Mat_Composite*)mat->data; 359 Mat_CompositeLink ilink,next = shell->head; 360 PetscErrorCode ierr; 361 362 PetscFunctionBegin; 363 ierr = PetscNewLog(mat,&ilink);CHKERRQ(ierr); 364 ilink->next = 0; 365 ierr = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr); 366 ilink->mat = smat; 367 368 if (!next) shell->head = ilink; 369 else { 370 while (next->next) { 371 next = next->next; 372 } 373 next->next = ilink; 374 ilink->prev = next; 375 } 376 shell->tail = ilink; 377 shell->nmat += 1; 378 PetscFunctionReturn(0); 379 } 380 381 /*@ 382 MatCompositeAddMat - Add another matrix to a composite matrix. 383 384 Collective on Mat 385 386 Input Parameters: 387 + mat - the composite matrix 388 - smat - the partial matrix 389 390 Level: advanced 391 392 .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE 393 @*/ 394 PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat) 395 { 396 PetscErrorCode ierr; 397 398 PetscFunctionBegin; 399 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 400 PetscValidHeaderSpecific(smat,MAT_CLASSID,2); 401 ierr = PetscUseMethod(mat,"MatCompositeAddMat_C",(Mat,Mat),(mat,smat));CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 405 static PetscErrorCode MatCompositeSetType_Composite(Mat mat,MatCompositeType type) 406 { 407 Mat_Composite *b = (Mat_Composite*)mat->data; 408 409 PetscFunctionBegin; 410 b->type = type; 411 if (type == MAT_COMPOSITE_MULTIPLICATIVE) { 412 mat->ops->getdiagonal = 0; 413 mat->ops->mult = MatMult_Composite_Multiplicative; 414 mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative; 415 } else { 416 mat->ops->getdiagonal = MatGetDiagonal_Composite; 417 mat->ops->mult = MatMult_Composite; 418 mat->ops->multtranspose = MatMultTranspose_Composite; 419 } 420 PetscFunctionReturn(0); 421 } 422 423 /*@ 424 MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product. 425 426 Logically Collective on Mat 427 428 Input Parameters: 429 . mat - the composite matrix 430 431 Level: advanced 432 433 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeGetType(), MATCOMPOSITE 434 435 @*/ 436 PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type) 437 { 438 PetscErrorCode ierr; 439 440 PetscFunctionBegin; 441 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 442 PetscValidLogicalCollectiveEnum(mat,type,2); 443 ierr = PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));CHKERRQ(ierr); 444 PetscFunctionReturn(0); 445 } 446 447 static PetscErrorCode MatCompositeGetType_Composite(Mat mat,MatCompositeType *type) 448 { 449 Mat_Composite *b = (Mat_Composite*)mat->data; 450 451 PetscFunctionBegin; 452 *type = b->type; 453 PetscFunctionReturn(0); 454 } 455 456 /*@ 457 MatCompositeGetType - Returns type of composite. 458 459 Not Collective 460 461 Input Parameter: 462 . mat - the composite matrix 463 464 Output Parameter: 465 . type - type of composite 466 467 Level: advanced 468 469 .seealso: MatCreateComposite(), MatCompositeSetType(), MATCOMPOSITE 470 471 @*/ 472 PetscErrorCode MatCompositeGetType(Mat mat,MatCompositeType *type) 473 { 474 PetscErrorCode ierr; 475 476 PetscFunctionBegin; 477 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 478 PetscValidPointer(type,2); 479 ierr = PetscUseMethod(mat,"MatCompositeGetType_C",(Mat,MatCompositeType*),(mat,type));CHKERRQ(ierr); 480 PetscFunctionReturn(0); 481 } 482 483 static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat,MatStructure str) 484 { 485 Mat_Composite *b = (Mat_Composite*)mat->data; 486 487 PetscFunctionBegin; 488 b->structure = str; 489 PetscFunctionReturn(0); 490 } 491 492 /*@ 493 MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix. 494 495 Not Collective 496 497 Input Parameters: 498 + mat - the composite matrix 499 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN (default) or SUBSET_NONZERO_PATTERN 500 501 Level: advanced 502 503 Notes: 504 Information about the matrices structure is used in MatCompositeMerge() for additive composite matrix. 505 506 .seealso: MatAXPY(), MatCreateComposite(), MatCompositeMerge() MatCompositeGetMatStructure(), MATCOMPOSITE 507 508 @*/ 509 PetscErrorCode MatCompositeSetMatStructure(Mat mat,MatStructure str) 510 { 511 PetscErrorCode ierr; 512 513 PetscFunctionBegin; 514 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 515 ierr = PetscUseMethod(mat,"MatCompositeSetMatStructure_C",(Mat,MatStructure),(mat,str));CHKERRQ(ierr); 516 PetscFunctionReturn(0); 517 } 518 519 static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat,MatStructure *str) 520 { 521 Mat_Composite *b = (Mat_Composite*)mat->data; 522 523 PetscFunctionBegin; 524 *str = b->structure; 525 PetscFunctionReturn(0); 526 } 527 528 /*@ 529 MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix. 530 531 Not Collective 532 533 Input Parameter: 534 . mat - the composite matrix 535 536 Output Parameter: 537 . str - structure of the matrices 538 539 Level: advanced 540 541 .seealso: MatCreateComposite(), MatCompositeSetMatStructure(), MATCOMPOSITE 542 543 @*/ 544 PetscErrorCode MatCompositeGetMatStructure(Mat mat,MatStructure *str) 545 { 546 PetscErrorCode ierr; 547 548 PetscFunctionBegin; 549 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 550 PetscValidPointer(str,2); 551 ierr = PetscUseMethod(mat,"MatCompositeGetMatStructure_C",(Mat,MatStructure*),(mat,str));CHKERRQ(ierr); 552 PetscFunctionReturn(0); 553 } 554 555 static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat,MatCompositeMergeType type) 556 { 557 Mat_Composite *shell = (Mat_Composite*)mat->data; 558 559 PetscFunctionBegin; 560 shell->mergetype = type; 561 PetscFunctionReturn(0); 562 } 563 564 /*@ 565 MatCompositeSetMergeType - Sets order of MatCompositeMerge(). 566 567 Logically Collective on Mat 568 569 Input Parameter: 570 + mat - the composite matrix 571 - type - MAT_COMPOSITE_MERGE RIGHT (default) to start merge from right with the first added matrix (mat[0]), 572 MAT_COMPOSITE_MERGE_LEFT to start merge from left with the last added matrix (mat[nmat-1]) 573 574 Level: advanced 575 576 Notes: 577 The resulting matrix is the same regardles of the MergeType. Only the order of operation is changed. 578 If set to MAT_COMPOSITE_MERGE_RIGHT the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0]))) 579 otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0]. 580 581 .seealso: MatCreateComposite(), MatCompositeMerge(), MATCOMPOSITE 582 583 @*/ 584 PetscErrorCode MatCompositeSetMergeType(Mat mat,MatCompositeMergeType type) 585 { 586 PetscErrorCode ierr; 587 588 PetscFunctionBegin; 589 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 590 PetscValidLogicalCollectiveEnum(mat,type,2); 591 ierr = PetscUseMethod(mat,"MatCompositeSetMergeType_C",(Mat,MatCompositeMergeType),(mat,type));CHKERRQ(ierr); 592 PetscFunctionReturn(0); 593 } 594 595 static PetscErrorCode MatCompositeMerge_Composite(Mat mat) 596 { 597 Mat_Composite *shell = (Mat_Composite*)mat->data; 598 Mat_CompositeLink next = shell->head, prev = shell->tail; 599 PetscErrorCode ierr; 600 Mat tmat,newmat; 601 Vec left,right; 602 PetscScalar scale; 603 604 PetscFunctionBegin; 605 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 606 607 PetscFunctionBegin; 608 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 609 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 610 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 611 while ((next = next->next)) { 612 ierr = MatAXPY(tmat,1.0,next->mat,shell->structure);CHKERRQ(ierr); 613 } 614 } else { 615 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 616 while ((prev = prev->prev)) { 617 ierr = MatAXPY(tmat,1.0,prev->mat,shell->structure);CHKERRQ(ierr); 618 } 619 } 620 } else { 621 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 622 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 623 while ((next = next->next)) { 624 ierr = MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 625 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 626 tmat = newmat; 627 } 628 } else { 629 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 630 while ((prev = prev->prev)) { 631 ierr = MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 632 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 633 tmat = newmat; 634 } 635 } 636 } 637 638 scale = shell->scale; 639 if ((left = shell->left)) {ierr = PetscObjectReference((PetscObject)left);CHKERRQ(ierr);} 640 if ((right = shell->right)) {ierr = PetscObjectReference((PetscObject)right);CHKERRQ(ierr);} 641 642 ierr = MatHeaderReplace(mat,&tmat);CHKERRQ(ierr); 643 644 ierr = MatDiagonalScale(mat,left,right);CHKERRQ(ierr); 645 ierr = MatScale(mat,scale);CHKERRQ(ierr); 646 ierr = VecDestroy(&left);CHKERRQ(ierr); 647 ierr = VecDestroy(&right);CHKERRQ(ierr); 648 PetscFunctionReturn(0); 649 } 650 651 /*@ 652 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 653 by summing or computing the product of all the matrices inside the composite matrix. 654 655 Collective 656 657 Input Parameters: 658 . mat - the composite matrix 659 660 661 Options Database Keys: 662 + -mat_composite_merge - merge in MatAssemblyEnd() 663 - -mat_composite_merge_type - set merge direction 664 665 Level: advanced 666 667 Notes: 668 The MatType of the resulting matrix will be the same as the MatType of the FIRST 669 matrix in the composite matrix. 670 671 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeSetMatStructure(), MatCompositeSetMergeType(), MATCOMPOSITE 672 673 @*/ 674 PetscErrorCode MatCompositeMerge(Mat mat) 675 { 676 PetscErrorCode ierr; 677 678 PetscFunctionBegin; 679 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 680 ierr = PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));CHKERRQ(ierr); 681 PetscFunctionReturn(0); 682 } 683 684 static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat,PetscInt *nmat) 685 { 686 Mat_Composite *shell = (Mat_Composite*)mat->data; 687 688 PetscFunctionBegin; 689 *nmat = shell->nmat; 690 PetscFunctionReturn(0); 691 } 692 693 /*@ 694 MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix. 695 696 Not Collective 697 698 Input Parameter: 699 . mat - the composite matrix 700 701 Output Parameter: 702 . nmat - number of matrices in the composite matrix 703 704 Level: advanced 705 706 .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE 707 708 @*/ 709 PetscErrorCode MatCompositeGetNumberMat(Mat mat,PetscInt *nmat) 710 { 711 PetscErrorCode ierr; 712 713 PetscFunctionBegin; 714 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 715 PetscValidPointer(nmat,2); 716 ierr = PetscUseMethod(mat,"MatCompositeGetNumberMat_C",(Mat,PetscInt*),(mat,nmat));CHKERRQ(ierr); 717 PetscFunctionReturn(0); 718 } 719 720 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai) 721 { 722 Mat_Composite *shell = (Mat_Composite*)mat->data; 723 Mat_CompositeLink ilink; 724 PetscInt k; 725 726 PetscFunctionBegin; 727 if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat); 728 ilink = shell->head; 729 for (k=0; k<i; k++) { 730 ilink = ilink->next; 731 } 732 *Ai = ilink->mat; 733 PetscFunctionReturn(0); 734 } 735 736 /*@ 737 MatCompositeGetMat - Returns the ith matrix from the composite matrix. 738 739 Logically Collective on Mat 740 741 Input Parameter: 742 + mat - the composite matrix 743 - i - the number of requested matrix 744 745 Output Parameter: 746 . Ai - ith matrix in composite 747 748 Level: advanced 749 750 .seealso: MatCreateComposite(), MatCompositeGetNumberMat(), MatCompositeAddMat(), MATCOMPOSITE 751 752 @*/ 753 PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai) 754 { 755 PetscErrorCode ierr; 756 757 PetscFunctionBegin; 758 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 759 PetscValidLogicalCollectiveInt(mat,i,2); 760 PetscValidPointer(Ai,3); 761 ierr = PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));CHKERRQ(ierr); 762 PetscFunctionReturn(0); 763 } 764 765 static struct _MatOps MatOps_Values = {0, 766 0, 767 0, 768 MatMult_Composite, 769 MatMultAdd_Composite, 770 /* 5*/ MatMultTranspose_Composite, 771 MatMultTransposeAdd_Composite, 772 0, 773 0, 774 0, 775 /* 10*/ 0, 776 0, 777 0, 778 0, 779 0, 780 /* 15*/ 0, 781 0, 782 MatGetDiagonal_Composite, 783 MatDiagonalScale_Composite, 784 0, 785 /* 20*/ 0, 786 MatAssemblyEnd_Composite, 787 0, 788 0, 789 /* 24*/ 0, 790 0, 791 0, 792 0, 793 0, 794 /* 29*/ 0, 795 0, 796 0, 797 0, 798 0, 799 /* 34*/ 0, 800 0, 801 0, 802 0, 803 0, 804 /* 39*/ 0, 805 0, 806 0, 807 0, 808 0, 809 /* 44*/ 0, 810 MatScale_Composite, 811 MatShift_Basic, 812 0, 813 0, 814 /* 49*/ 0, 815 0, 816 0, 817 0, 818 0, 819 /* 54*/ 0, 820 0, 821 0, 822 0, 823 0, 824 /* 59*/ 0, 825 MatDestroy_Composite, 826 0, 827 0, 828 0, 829 /* 64*/ 0, 830 0, 831 0, 832 0, 833 0, 834 /* 69*/ 0, 835 0, 836 0, 837 0, 838 0, 839 /* 74*/ 0, 840 0, 841 MatSetFromOptions_Composite, 842 0, 843 0, 844 /* 79*/ 0, 845 0, 846 0, 847 0, 848 0, 849 /* 84*/ 0, 850 0, 851 0, 852 0, 853 0, 854 /* 89*/ 0, 855 0, 856 0, 857 0, 858 0, 859 /* 94*/ 0, 860 0, 861 0, 862 0, 863 0, 864 /*99*/ 0, 865 0, 866 0, 867 0, 868 0, 869 /*104*/ 0, 870 0, 871 0, 872 0, 873 0, 874 /*109*/ 0, 875 0, 876 0, 877 0, 878 0, 879 /*114*/ 0, 880 0, 881 0, 882 0, 883 0, 884 /*119*/ 0, 885 0, 886 0, 887 0, 888 0, 889 /*124*/ 0, 890 0, 891 0, 892 0, 893 0, 894 /*129*/ 0, 895 0, 896 0, 897 0, 898 0, 899 /*134*/ 0, 900 0, 901 0, 902 0, 903 0, 904 /*139*/ 0, 905 0, 906 0 907 }; 908 909 /*MC 910 MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices. 911 The matrices need to have a correct size and parallel layout for the sum or product to be valid. 912 913 Notes: 914 to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE); 915 916 Level: advanced 917 918 .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeSetType(), MatCompositeGetType(), MatCompositeSetMatStructure(), MatCompositeGetMatStructure(), MatCompositeMerge(), MatCompositeSetMergeType(), MatCompositeGetNumberMat(), MatCompositeGetMat() 919 M*/ 920 921 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A) 922 { 923 Mat_Composite *b; 924 PetscErrorCode ierr; 925 926 PetscFunctionBegin; 927 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 928 A->data = (void*)b; 929 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 930 931 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 932 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 933 934 A->assembled = PETSC_TRUE; 935 A->preallocated = PETSC_TRUE; 936 b->type = MAT_COMPOSITE_ADDITIVE; 937 b->scale = 1.0; 938 b->nmat = 0; 939 b->merge = PETSC_FALSE; 940 b->mergetype = MAT_COMPOSITE_MERGE_RIGHT; 941 b->structure = DIFFERENT_NONZERO_PATTERN; 942 943 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);CHKERRQ(ierr); 944 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);CHKERRQ(ierr); 945 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetType_C",MatCompositeGetType_Composite);CHKERRQ(ierr); 946 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMergeType_C",MatCompositeSetMergeType_Composite);CHKERRQ(ierr); 947 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMatStructure_C",MatCompositeSetMatStructure_Composite);CHKERRQ(ierr); 948 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMatStructure_C",MatCompositeGetMatStructure_Composite);CHKERRQ(ierr); 949 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);CHKERRQ(ierr); 950 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNumberMat_C",MatCompositeGetNumberMat_Composite);CHKERRQ(ierr); 951 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);CHKERRQ(ierr); 952 953 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 954 PetscFunctionReturn(0); 955 } 956 957