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