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