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