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 Logically Collective on Mat 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 PetscValidLogicalCollectiveEnum(mat,type,2); 442 ierr = PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));CHKERRQ(ierr); 443 PetscFunctionReturn(0); 444 } 445 446 static PetscErrorCode MatCompositeGetType_Composite(Mat mat,MatCompositeType *type) 447 { 448 Mat_Composite *b = (Mat_Composite*)mat->data; 449 450 PetscFunctionBegin; 451 *type = b->type; 452 PetscFunctionReturn(0); 453 } 454 455 /*@ 456 MatCompositeGetType - Returns type of composite. 457 458 Not Collective 459 460 Input Parameter: 461 . mat - the composite matrix 462 463 Output Parameter: 464 . type - type of composite 465 466 Level: advanced 467 468 .seealso: MatCreateComposite(), MatCompositeSetType(), MATCOMPOSITE 469 470 @*/ 471 PetscErrorCode MatCompositeGetType(Mat mat,MatCompositeType *type) 472 { 473 PetscErrorCode ierr; 474 475 PetscFunctionBegin; 476 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 477 PetscValidPointer(type,2); 478 ierr = PetscUseMethod(mat,"MatCompositeGetType_C",(Mat,MatCompositeType*),(mat,type));CHKERRQ(ierr); 479 PetscFunctionReturn(0); 480 } 481 482 static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat,MatCompositeMergeType type) 483 { 484 Mat_Composite *shell = (Mat_Composite*)mat->data; 485 486 PetscFunctionBegin; 487 shell->mergetype = type; 488 PetscFunctionReturn(0); 489 } 490 491 /*@ 492 MatCompositeSetMergeType - Sets order of MatCompositeMerge(). 493 494 Logically Collective on Mat 495 496 Input Parameter: 497 + mat - the composite matrix 498 - type - MAT_COMPOSITE_MERGE RIGHT (default) to start merge from right with the first added matrix (mat[0]), 499 MAT_COMPOSITE_MERGE_LEFT to start merge from left with the last added matrix (mat[nmat-1]) 500 501 Level: advanced 502 503 Notes: 504 The resulting matrix is the same regardles of the MergeType. Only the order of operation is changed. 505 If set to MAT_COMPOSITE_MERGE_RIGHT the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0]))) 506 otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0]. 507 508 .seealso: MatCreateComposite(), MatCompositeMerge(), MATCOMPOSITE 509 510 @*/ 511 PetscErrorCode MatCompositeSetMergeType(Mat mat,MatCompositeMergeType type) 512 { 513 PetscErrorCode ierr; 514 515 PetscFunctionBegin; 516 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 517 PetscValidLogicalCollectiveEnum(mat,type,2); 518 ierr = PetscUseMethod(mat,"MatCompositeSetMergeType_C",(Mat,MatCompositeMergeType),(mat,type));CHKERRQ(ierr); 519 PetscFunctionReturn(0); 520 } 521 522 static PetscErrorCode MatCompositeMerge_Composite(Mat mat) 523 { 524 Mat_Composite *shell = (Mat_Composite*)mat->data; 525 Mat_CompositeLink next = shell->head, prev = shell->tail; 526 PetscErrorCode ierr; 527 Mat tmat,newmat; 528 Vec left,right; 529 PetscScalar scale; 530 531 PetscFunctionBegin; 532 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 533 534 PetscFunctionBegin; 535 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 536 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 537 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 538 while ((next = next->next)) { 539 ierr = MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 540 } 541 } else { 542 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 543 while ((prev = prev->prev)) { 544 ierr = MatAXPY(tmat,1.0,prev->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 545 } 546 } 547 } else { 548 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 549 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 550 while ((next = next->next)) { 551 ierr = MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 552 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 553 tmat = newmat; 554 } 555 } else { 556 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 557 while ((prev = prev->prev)) { 558 ierr = MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 559 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 560 tmat = newmat; 561 } 562 } 563 } 564 565 scale = shell->scale; 566 if ((left = shell->left)) {ierr = PetscObjectReference((PetscObject)left);CHKERRQ(ierr);} 567 if ((right = shell->right)) {ierr = PetscObjectReference((PetscObject)right);CHKERRQ(ierr);} 568 569 ierr = MatHeaderReplace(mat,&tmat);CHKERRQ(ierr); 570 571 ierr = MatDiagonalScale(mat,left,right);CHKERRQ(ierr); 572 ierr = MatScale(mat,scale);CHKERRQ(ierr); 573 ierr = VecDestroy(&left);CHKERRQ(ierr); 574 ierr = VecDestroy(&right);CHKERRQ(ierr); 575 PetscFunctionReturn(0); 576 } 577 578 /*@ 579 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 580 by summing or computing the product of all the matrices inside the composite matrix. 581 582 Collective 583 584 Input Parameters: 585 . mat - the composite matrix 586 587 588 Options Database Keys: 589 + -mat_composite_merge - merge in MatAssemblyEnd() 590 - -mat_composite_merge_type - set merge direction 591 592 Level: advanced 593 594 Notes: 595 The MatType of the resulting matrix will be the same as the MatType of the FIRST 596 matrix in the composite matrix. 597 598 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeSetMergeType(), MATCOMPOSITE 599 600 @*/ 601 PetscErrorCode MatCompositeMerge(Mat mat) 602 { 603 PetscErrorCode ierr; 604 605 PetscFunctionBegin; 606 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 607 ierr = PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));CHKERRQ(ierr); 608 PetscFunctionReturn(0); 609 } 610 611 static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat,PetscInt *nmat) 612 { 613 Mat_Composite *shell = (Mat_Composite*)mat->data; 614 615 PetscFunctionBegin; 616 *nmat = shell->nmat; 617 PetscFunctionReturn(0); 618 } 619 620 /*@ 621 MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix. 622 623 Not Collective 624 625 Input Parameter: 626 . mat - the composite matrix 627 628 Output Parameter: 629 . nmat - number of matrices in the composite matrix 630 631 Level: advanced 632 633 .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE 634 635 @*/ 636 PetscErrorCode MatCompositeGetNumberMat(Mat mat,PetscInt *nmat) 637 { 638 PetscErrorCode ierr; 639 640 PetscFunctionBegin; 641 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 642 PetscValidPointer(nmat,2); 643 ierr = PetscUseMethod(mat,"MatCompositeGetNumberMat_C",(Mat,PetscInt*),(mat,nmat));CHKERRQ(ierr); 644 PetscFunctionReturn(0); 645 } 646 647 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai) 648 { 649 Mat_Composite *shell = (Mat_Composite*)mat->data; 650 Mat_CompositeLink ilink; 651 PetscInt k; 652 653 PetscFunctionBegin; 654 if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat); 655 ilink = shell->head; 656 for (k=0; k<i; k++) { 657 ilink = ilink->next; 658 } 659 *Ai = ilink->mat; 660 PetscFunctionReturn(0); 661 } 662 663 /*@ 664 MatCompositeGetMat - Returns the ith matrix from the composite matrix. 665 666 Logically Collective on Mat 667 668 Input Parameter: 669 + mat - the composite matrix 670 - i - the number of requested matrix 671 672 Output Parameter: 673 . Ai - ith matrix in composite 674 675 Level: advanced 676 677 .seealso: MatCreateComposite(), MatCompositeGetNumberMat(), MatCompositeAddMat(), MATCOMPOSITE 678 679 @*/ 680 PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai) 681 { 682 PetscErrorCode ierr; 683 684 PetscFunctionBegin; 685 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 686 PetscValidLogicalCollectiveInt(mat,i,2); 687 PetscValidPointer(Ai,3); 688 ierr = PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));CHKERRQ(ierr); 689 PetscFunctionReturn(0); 690 } 691 692 static struct _MatOps MatOps_Values = {0, 693 0, 694 0, 695 MatMult_Composite, 696 MatMultAdd_Composite, 697 /* 5*/ MatMultTranspose_Composite, 698 MatMultTransposeAdd_Composite, 699 0, 700 0, 701 0, 702 /* 10*/ 0, 703 0, 704 0, 705 0, 706 0, 707 /* 15*/ 0, 708 0, 709 MatGetDiagonal_Composite, 710 MatDiagonalScale_Composite, 711 0, 712 /* 20*/ 0, 713 MatAssemblyEnd_Composite, 714 0, 715 0, 716 /* 24*/ 0, 717 0, 718 0, 719 0, 720 0, 721 /* 29*/ 0, 722 0, 723 0, 724 0, 725 0, 726 /* 34*/ 0, 727 0, 728 0, 729 0, 730 0, 731 /* 39*/ 0, 732 0, 733 0, 734 0, 735 0, 736 /* 44*/ 0, 737 MatScale_Composite, 738 MatShift_Basic, 739 0, 740 0, 741 /* 49*/ 0, 742 0, 743 0, 744 0, 745 0, 746 /* 54*/ 0, 747 0, 748 0, 749 0, 750 0, 751 /* 59*/ 0, 752 MatDestroy_Composite, 753 0, 754 0, 755 0, 756 /* 64*/ 0, 757 0, 758 0, 759 0, 760 0, 761 /* 69*/ 0, 762 0, 763 0, 764 0, 765 0, 766 /* 74*/ 0, 767 0, 768 MatSetFromOptions_Composite, 769 0, 770 0, 771 /* 79*/ 0, 772 0, 773 0, 774 0, 775 0, 776 /* 84*/ 0, 777 0, 778 0, 779 0, 780 0, 781 /* 89*/ 0, 782 0, 783 0, 784 0, 785 0, 786 /* 94*/ 0, 787 0, 788 0, 789 0, 790 0, 791 /*99*/ 0, 792 0, 793 0, 794 0, 795 0, 796 /*104*/ 0, 797 0, 798 0, 799 0, 800 0, 801 /*109*/ 0, 802 0, 803 0, 804 0, 805 0, 806 /*114*/ 0, 807 0, 808 0, 809 0, 810 0, 811 /*119*/ 0, 812 0, 813 0, 814 0, 815 0, 816 /*124*/ 0, 817 0, 818 0, 819 0, 820 0, 821 /*129*/ 0, 822 0, 823 0, 824 0, 825 0, 826 /*134*/ 0, 827 0, 828 0, 829 0, 830 0, 831 /*139*/ 0, 832 0, 833 0 834 }; 835 836 /*MC 837 MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices. 838 The matrices need to have a correct size and parallel layout for the sum or product to be valid. 839 840 Notes: 841 to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE); 842 843 Level: advanced 844 845 .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeSetType(), MatCompositeGetType(), MatCompositeMerge(), MatCompositeSetMergeType(), MatCompositeGetNumberMat(), MatCompositeGetMat() 846 M*/ 847 848 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A) 849 { 850 Mat_Composite *b; 851 PetscErrorCode ierr; 852 853 PetscFunctionBegin; 854 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 855 A->data = (void*)b; 856 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 857 858 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 859 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 860 861 A->assembled = PETSC_TRUE; 862 A->preallocated = PETSC_TRUE; 863 b->type = MAT_COMPOSITE_ADDITIVE; 864 b->scale = 1.0; 865 b->nmat = 0; 866 b->merge = PETSC_FALSE; 867 b->mergetype = MAT_COMPOSITE_MERGE_RIGHT; 868 869 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);CHKERRQ(ierr); 870 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);CHKERRQ(ierr); 871 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetType_C",MatCompositeGetType_Composite);CHKERRQ(ierr); 872 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMergeType_C",MatCompositeSetMergeType_Composite);CHKERRQ(ierr); 873 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);CHKERRQ(ierr); 874 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNumberMat_C",MatCompositeGetNumberMat_Composite);CHKERRQ(ierr); 875 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);CHKERRQ(ierr); 876 877 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 878 PetscFunctionReturn(0); 879 } 880 881