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