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 of zero 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() 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 if (type == MAT_COMPOSITE_MULTIPLICATIVE) { 391 mat->ops->getdiagonal = 0; 392 mat->ops->mult = MatMult_Composite_Multiplicative; 393 mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative; 394 b->type = MAT_COMPOSITE_MULTIPLICATIVE; 395 } else { 396 mat->ops->getdiagonal = MatGetDiagonal_Composite; 397 mat->ops->mult = MatMult_Composite; 398 mat->ops->multtranspose = MatMultTranspose_Composite; 399 b->type = MAT_COMPOSITE_ADDITIVE; 400 } 401 PetscFunctionReturn(0); 402 } 403 404 /*@ 405 MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product 406 407 Collective on MPI_Comm 408 409 Input Parameters: 410 . mat - the composite matrix 411 412 413 Level: advanced 414 415 Notes: 416 The MatType of the resulting matrix will be the same as the MatType of the FIRST 417 matrix in the composite matrix. 418 419 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeGetType(), MATCOMPOSITE 420 421 @*/ 422 PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type) 423 { 424 PetscErrorCode ierr; 425 426 PetscFunctionBegin; 427 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 428 ierr = PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));CHKERRQ(ierr); 429 PetscFunctionReturn(0); 430 } 431 432 static PetscErrorCode MatCompositeGetType_Composite(Mat mat,MatCompositeType *type) 433 { 434 Mat_Composite *b = (Mat_Composite*)mat->data; 435 436 PetscFunctionBegin; 437 *type = b->type; 438 PetscFunctionReturn(0); 439 } 440 441 /*@ 442 MatCompositeGetType - Returns type of composite. 443 444 Not Collective 445 446 Input Parameter: 447 . mat - the composite matrix 448 449 Output Parameter: 450 . type - type of composite 451 452 Level: advanced 453 454 .seealso: MatCreateComposite(), MatCompositeSetType(), MATCOMPOSITE 455 456 @*/ 457 PetscErrorCode MatCompositeGetType(Mat mat,MatCompositeType *type) 458 { 459 PetscErrorCode ierr; 460 461 PetscFunctionBegin; 462 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 463 PetscValidPointer(type,2); 464 ierr = PetscUseMethod(mat,"MatCompositeGetType_C",(Mat,MatCompositeType*),(mat,type));CHKERRQ(ierr); 465 PetscFunctionReturn(0); 466 } 467 468 static PetscErrorCode MatCompositeMerge_Composite(Mat mat) 469 { 470 Mat_Composite *shell = (Mat_Composite*)mat->data; 471 Mat_CompositeLink next = shell->head, prev = shell->tail; 472 PetscErrorCode ierr; 473 Mat tmat,newmat; 474 Vec left,right; 475 PetscScalar scale; 476 477 PetscFunctionBegin; 478 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 479 480 PetscFunctionBegin; 481 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 482 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 483 while ((next = next->next)) { 484 ierr = MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 485 } 486 } else { 487 if (shell->mergefromright) { 488 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 489 while ((next = next->next)) { 490 ierr = MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 491 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 492 tmat = newmat; 493 } 494 } else { 495 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 496 while ((prev = prev->prev)) { 497 ierr = MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 498 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 499 tmat = newmat; 500 } 501 } 502 } 503 504 scale = shell->scale; 505 if ((left = shell->left)) {ierr = PetscObjectReference((PetscObject)left);CHKERRQ(ierr);} 506 if ((right = shell->right)) {ierr = PetscObjectReference((PetscObject)right);CHKERRQ(ierr);} 507 508 ierr = MatHeaderReplace(mat,&tmat);CHKERRQ(ierr); 509 510 ierr = MatDiagonalScale(mat,left,right);CHKERRQ(ierr); 511 ierr = MatScale(mat,scale);CHKERRQ(ierr); 512 ierr = VecDestroy(&left);CHKERRQ(ierr); 513 ierr = VecDestroy(&right);CHKERRQ(ierr); 514 PetscFunctionReturn(0); 515 } 516 517 /*@ 518 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 519 by summing all the matrices inside the composite matrix. 520 521 Collective on MPI_Comm 522 523 Input Parameters: 524 . mat - the composite matrix 525 526 527 Options Database: 528 . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked) 529 530 Level: advanced 531 532 Notes: 533 The MatType of the resulting matrix will be the same as the MatType of the FIRST 534 matrix in the composite matrix. 535 536 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeSetMergeFromRight(), MATCOMPOSITE 537 538 @*/ 539 PetscErrorCode MatCompositeMerge(Mat mat) 540 { 541 PetscErrorCode ierr; 542 543 PetscFunctionBegin; 544 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 545 ierr = PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 static PetscErrorCode MatCompositeGetNMat_Composite(Mat mat,PetscInt *nmat) 550 { 551 Mat_Composite *shell = (Mat_Composite*)mat->data; 552 553 PetscFunctionBegin; 554 *nmat = shell->nmat; 555 PetscFunctionReturn(0); 556 } 557 558 /*@ 559 MatCompositeGetNMat - Returns the number of matrices in composite. 560 561 Not Collective 562 563 Input Parameter: 564 . mat - the composite matrix 565 566 Output Parameter: 567 + size - the local size 568 - nmat - number of matrices in composite 569 570 Level: advanced 571 572 .seealso: MatCreateComposite(), MatCompositeGetMat() 573 574 @*/ 575 PetscErrorCode MatCompositeGetNMat(Mat mat,PetscInt *nmat) 576 { 577 PetscErrorCode ierr; 578 579 PetscFunctionBegin; 580 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 581 PetscValidPointer(nmat,2); 582 ierr = PetscUseMethod(mat,"MatCompositeGetNMat_C",(Mat,PetscInt*),(mat,nmat));CHKERRQ(ierr); 583 PetscFunctionReturn(0); 584 } 585 586 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai) 587 { 588 Mat_Composite *shell = (Mat_Composite*)mat->data; 589 Mat_CompositeLink ilink; 590 PetscInt k; 591 592 PetscFunctionBegin; 593 if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat); 594 ilink = shell->head; 595 for (k=0; k<i; k++) { 596 ilink = ilink->next; 597 } 598 *Ai = ilink->mat; 599 PetscFunctionReturn(0); 600 } 601 602 /*@ 603 MatCompositeGetMat - Returns the ith matrix from composite. 604 605 Logically Collective on Mat 606 607 Input Parameter: 608 + mat - the composite matrix 609 - i - the number of requested matrix 610 611 Output Parameter: 612 . Ai - ith matrix in composite 613 614 Level: advanced 615 616 .seealso: MatCreateComposite(), MatCompositeGetNMat() 617 618 @*/ 619 PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai) 620 { 621 PetscErrorCode ierr; 622 623 PetscFunctionBegin; 624 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 625 PetscValidLogicalCollectiveInt(mat,i,2); 626 PetscValidPointer(Ai,3); 627 ierr = PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));CHKERRQ(ierr); 628 PetscFunctionReturn(0); 629 } 630 631 static PetscErrorCode MatCompositeSetMergeFromRight_Composite(Mat mat,PetscBool flg) 632 { 633 Mat_Composite *shell = (Mat_Composite*)mat->data; 634 635 PetscFunctionBegin; 636 shell->mergefromright = flg; 637 PetscFunctionReturn(0); 638 } 639 640 /*@ 641 MatCompositeSetMergeFromRight - Sets MatCompositeMerge() to start from right. 642 643 Logically Collective on Mat 644 645 Input Parameter: 646 + mat - the composite matrix 647 - flg - if true (default) the merge starts with the first added matrix (mat[0]) 648 649 Level: advanced 650 651 .seealso: MatCreateComposite(), MatCompositeMerge() 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(), MatCompositeGetNmat(), 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,"MatCompositeGetNMat_C",MatCompositeGetNMat_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