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 } Mat_Composite; 19 20 #undef __FUNCT__ 21 #define __FUNCT__ "MatDestroy_Composite" 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 #undef __FUNCT__ 48 #define __FUNCT__ "MatMult_Composite_Multiplicative" 49 PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y) 50 { 51 Mat_Composite *shell = (Mat_Composite*)A->data; 52 Mat_CompositeLink next = shell->head; 53 PetscErrorCode ierr; 54 Vec in,out; 55 56 PetscFunctionBegin; 57 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 58 in = x; 59 if (shell->right) { 60 if (!shell->rightwork) { 61 ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr); 62 } 63 ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr); 64 in = shell->rightwork; 65 } 66 while (next->next) { 67 if (!next->work) { /* should reuse previous work if the same size */ 68 ierr = MatGetVecs(next->mat,NULL,&next->work);CHKERRQ(ierr); 69 } 70 out = next->work; 71 ierr = MatMult(next->mat,in,out);CHKERRQ(ierr); 72 in = out; 73 next = next->next; 74 } 75 ierr = MatMult(next->mat,in,y);CHKERRQ(ierr); 76 if (shell->left) { 77 ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr); 78 } 79 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "MatMultTranspose_Composite_Multiplicative" 85 PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y) 86 { 87 Mat_Composite *shell = (Mat_Composite*)A->data; 88 Mat_CompositeLink tail = shell->tail; 89 PetscErrorCode ierr; 90 Vec in,out; 91 92 PetscFunctionBegin; 93 if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 94 in = x; 95 if (shell->left) { 96 if (!shell->leftwork) { 97 ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr); 98 } 99 ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr); 100 in = shell->leftwork; 101 } 102 while (tail->prev) { 103 if (!tail->prev->work) { /* should reuse previous work if the same size */ 104 ierr = MatGetVecs(tail->mat,NULL,&tail->prev->work);CHKERRQ(ierr); 105 } 106 out = tail->prev->work; 107 ierr = MatMultTranspose(tail->mat,in,out);CHKERRQ(ierr); 108 in = out; 109 tail = tail->prev; 110 } 111 ierr = MatMultTranspose(tail->mat,in,y);CHKERRQ(ierr); 112 if (shell->right) { 113 ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr); 114 } 115 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "MatMult_Composite" 121 PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y) 122 { 123 Mat_Composite *shell = (Mat_Composite*)A->data; 124 Mat_CompositeLink next = shell->head; 125 PetscErrorCode ierr; 126 Vec in; 127 128 PetscFunctionBegin; 129 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 130 in = x; 131 if (shell->right) { 132 if (!shell->rightwork) { 133 ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr); 134 } 135 ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr); 136 in = shell->rightwork; 137 } 138 ierr = MatMult(next->mat,in,y);CHKERRQ(ierr); 139 while ((next = next->next)) { 140 ierr = MatMultAdd(next->mat,in,y,y);CHKERRQ(ierr); 141 } 142 if (shell->left) { 143 ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr); 144 } 145 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNCT__ 150 #define __FUNCT__ "MatMultTranspose_Composite" 151 PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y) 152 { 153 Mat_Composite *shell = (Mat_Composite*)A->data; 154 Mat_CompositeLink next = shell->head; 155 PetscErrorCode ierr; 156 Vec in; 157 158 PetscFunctionBegin; 159 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 160 in = x; 161 if (shell->left) { 162 if (!shell->leftwork) { 163 ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr); 164 } 165 ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr); 166 in = shell->leftwork; 167 } 168 ierr = MatMultTranspose(next->mat,in,y);CHKERRQ(ierr); 169 while ((next = next->next)) { 170 ierr = MatMultTransposeAdd(next->mat,in,y,y);CHKERRQ(ierr); 171 } 172 if (shell->right) { 173 ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr); 174 } 175 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 179 #undef __FUNCT__ 180 #define __FUNCT__ "MatGetDiagonal_Composite" 181 PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v) 182 { 183 Mat_Composite *shell = (Mat_Composite*)A->data; 184 Mat_CompositeLink next = shell->head; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 189 if (shell->right || shell->left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get diagonal if left or right scaling"); 190 191 ierr = MatGetDiagonal(next->mat,v);CHKERRQ(ierr); 192 if (next->next && !shell->work) { 193 ierr = VecDuplicate(v,&shell->work);CHKERRQ(ierr); 194 } 195 while ((next = next->next)) { 196 ierr = MatGetDiagonal(next->mat,shell->work);CHKERRQ(ierr); 197 ierr = VecAXPY(v,1.0,shell->work);CHKERRQ(ierr); 198 } 199 ierr = VecScale(v,shell->scale);CHKERRQ(ierr); 200 PetscFunctionReturn(0); 201 } 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "MatAssemblyEnd_Composite" 205 PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t) 206 { 207 PetscErrorCode ierr; 208 PetscBool flg = PETSC_FALSE; 209 210 PetscFunctionBegin; 211 ierr = PetscOptionsGetBool(((PetscObject)Y)->prefix,"-mat_composite_merge",&flg,NULL);CHKERRQ(ierr); 212 if (flg) { 213 ierr = MatCompositeMerge(Y);CHKERRQ(ierr); 214 } 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "MatScale_Composite" 220 PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha) 221 { 222 Mat_Composite *a = (Mat_Composite*)inA->data; 223 224 PetscFunctionBegin; 225 a->scale *= alpha; 226 PetscFunctionReturn(0); 227 } 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "MatDiagonalScale_Composite" 231 PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right) 232 { 233 Mat_Composite *a = (Mat_Composite*)inA->data; 234 PetscErrorCode ierr; 235 236 PetscFunctionBegin; 237 if (left) { 238 if (!a->left) { 239 ierr = VecDuplicate(left,&a->left);CHKERRQ(ierr); 240 ierr = VecCopy(left,a->left);CHKERRQ(ierr); 241 } else { 242 ierr = VecPointwiseMult(a->left,left,a->left);CHKERRQ(ierr); 243 } 244 } 245 if (right) { 246 if (!a->right) { 247 ierr = VecDuplicate(right,&a->right);CHKERRQ(ierr); 248 ierr = VecCopy(right,a->right);CHKERRQ(ierr); 249 } else { 250 ierr = VecPointwiseMult(a->right,right,a->right);CHKERRQ(ierr); 251 } 252 } 253 PetscFunctionReturn(0); 254 } 255 256 static struct _MatOps MatOps_Values = {0, 257 0, 258 0, 259 MatMult_Composite, 260 0, 261 /* 5*/ MatMultTranspose_Composite, 262 0, 263 0, 264 0, 265 0, 266 /* 10*/ 0, 267 0, 268 0, 269 0, 270 0, 271 /* 15*/ 0, 272 0, 273 MatGetDiagonal_Composite, 274 MatDiagonalScale_Composite, 275 0, 276 /* 20*/ 0, 277 MatAssemblyEnd_Composite, 278 0, 279 0, 280 /* 24*/ 0, 281 0, 282 0, 283 0, 284 0, 285 /* 29*/ 0, 286 0, 287 0, 288 0, 289 0, 290 /* 34*/ 0, 291 0, 292 0, 293 0, 294 0, 295 /* 39*/ 0, 296 0, 297 0, 298 0, 299 0, 300 /* 44*/ 0, 301 MatScale_Composite, 302 0, 303 0, 304 0, 305 /* 49*/ 0, 306 0, 307 0, 308 0, 309 0, 310 /* 54*/ 0, 311 0, 312 0, 313 0, 314 0, 315 /* 59*/ 0, 316 MatDestroy_Composite, 317 0, 318 0, 319 0, 320 /* 64*/ 0, 321 0, 322 0, 323 0, 324 0, 325 /* 69*/ 0, 326 0, 327 0, 328 0, 329 0, 330 /* 74*/ 0, 331 0, 332 0, 333 0, 334 0, 335 /* 79*/ 0, 336 0, 337 0, 338 0, 339 0, 340 /* 84*/ 0, 341 0, 342 0, 343 0, 344 0, 345 /* 89*/ 0, 346 0, 347 0, 348 0, 349 0, 350 /* 94*/ 0, 351 0, 352 0, 353 0, 354 0, 355 /*99*/ 0, 356 0, 357 0, 358 0, 359 0, 360 /*104*/ 0, 361 0, 362 0, 363 0, 364 0, 365 /*109*/ 0, 366 0, 367 0, 368 0, 369 0, 370 /*114*/ 0, 371 0, 372 0, 373 0, 374 0, 375 /*119*/ 0, 376 0, 377 0, 378 0, 379 0, 380 /*124*/ 0, 381 0, 382 0, 383 0, 384 0, 385 /*129*/ 0, 386 0, 387 0, 388 0, 389 0, 390 /*134*/ 0, 391 0, 392 0, 393 0, 394 0, 395 /*139*/ 0, 396 0 397 }; 398 399 /*MC 400 MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices (all matrices are of same size and parallel layout). 401 402 Notes: to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE); 403 404 Level: advanced 405 406 .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType 407 M*/ 408 409 #undef __FUNCT__ 410 #define __FUNCT__ "MatCreate_Composite" 411 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A) 412 { 413 Mat_Composite *b; 414 PetscErrorCode ierr; 415 416 PetscFunctionBegin; 417 ierr = PetscNewLog(A,Mat_Composite,&b);CHKERRQ(ierr); 418 A->data = (void*)b; 419 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 420 421 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 422 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 423 424 A->assembled = PETSC_TRUE; 425 A->preallocated = PETSC_TRUE; 426 b->type = MAT_COMPOSITE_ADDITIVE; 427 b->scale = 1.0; 428 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 429 PetscFunctionReturn(0); 430 } 431 432 #undef __FUNCT__ 433 #define __FUNCT__ "MatCreateComposite" 434 /*@C 435 MatCreateComposite - Creates a matrix as the sum of zero or more matrices 436 437 Collective on MPI_Comm 438 439 Input Parameters: 440 + comm - MPI communicator 441 . nmat - number of matrices to put in 442 - mats - the matrices 443 444 Output Parameter: 445 . A - the matrix 446 447 Level: advanced 448 449 Notes: 450 Alternative construction 451 $ MatCreate(comm,&mat); 452 $ MatSetSizes(mat,m,n,M,N); 453 $ MatSetType(mat,MATCOMPOSITE); 454 $ MatCompositeAddMat(mat,mats[0]); 455 $ .... 456 $ MatCompositeAddMat(mat,mats[nmat-1]); 457 $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 458 $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); 459 460 For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0] 461 462 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType 463 464 @*/ 465 PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat) 466 { 467 PetscErrorCode ierr; 468 PetscInt m,n,M,N,i; 469 470 PetscFunctionBegin; 471 if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix"); 472 PetscValidPointer(mat,3); 473 474 ierr = MatGetLocalSize(mats[0],&m,&n);CHKERRQ(ierr); 475 ierr = MatGetSize(mats[0],&M,&N);CHKERRQ(ierr); 476 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 477 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 478 ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr); 479 for (i=0; i<nmat; i++) { 480 ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr); 481 } 482 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 483 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 484 PetscFunctionReturn(0); 485 } 486 487 #undef __FUNCT__ 488 #define __FUNCT__ "MatCompositeAddMat" 489 /*@ 490 MatCompositeAddMat - add another matrix to a composite matrix 491 492 Collective on Mat 493 494 Input Parameters: 495 + mat - the composite matrix 496 - smat - the partial matrix 497 498 Level: advanced 499 500 .seealso: MatCreateComposite() 501 @*/ 502 PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat) 503 { 504 Mat_Composite *shell; 505 PetscErrorCode ierr; 506 Mat_CompositeLink ilink,next; 507 508 PetscFunctionBegin; 509 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 510 PetscValidHeaderSpecific(smat,MAT_CLASSID,2); 511 ierr = PetscNewLog(mat,struct _Mat_CompositeLink,&ilink);CHKERRQ(ierr); 512 ilink->next = 0; 513 ierr = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr); 514 ilink->mat = smat; 515 516 shell = (Mat_Composite*)mat->data; 517 next = shell->head; 518 if (!next) shell->head = ilink; 519 else { 520 while (next->next) { 521 next = next->next; 522 } 523 next->next = ilink; 524 ilink->prev = next; 525 } 526 shell->tail = ilink; 527 PetscFunctionReturn(0); 528 } 529 530 #undef __FUNCT__ 531 #define __FUNCT__ "MatCompositeSetType" 532 /*@C 533 MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product 534 535 Collective on MPI_Comm 536 537 Input Parameters: 538 . mat - the composite matrix 539 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(), MATCOMPOSITE 548 549 @*/ 550 PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type) 551 { 552 Mat_Composite *b = (Mat_Composite*)mat->data; 553 PetscBool flg; 554 PetscErrorCode ierr; 555 556 PetscFunctionBegin; 557 ierr = PetscObjectTypeCompare((PetscObject)mat,MATCOMPOSITE,&flg);CHKERRQ(ierr); 558 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Can only use with composite matrix"); 559 if (type == MAT_COMPOSITE_MULTIPLICATIVE) { 560 mat->ops->getdiagonal = 0; 561 mat->ops->mult = MatMult_Composite_Multiplicative; 562 mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative; 563 b->type = MAT_COMPOSITE_MULTIPLICATIVE; 564 } else { 565 mat->ops->getdiagonal = MatGetDiagonal_Composite; 566 mat->ops->mult = MatMult_Composite; 567 mat->ops->multtranspose = MatMultTranspose_Composite; 568 b->type = MAT_COMPOSITE_ADDITIVE; 569 } 570 PetscFunctionReturn(0); 571 } 572 573 574 #undef __FUNCT__ 575 #define __FUNCT__ "MatCompositeMerge" 576 /*@C 577 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 578 by summing all the matrices inside the composite matrix. 579 580 Collective on MPI_Comm 581 582 Input Parameters: 583 . mat - the composite matrix 584 585 586 Options Database: 587 . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked) 588 589 Level: advanced 590 591 Notes: 592 The MatType of the resulting matrix will be the same as the MatType of the FIRST 593 matrix in the composite matrix. 594 595 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE 596 597 @*/ 598 PetscErrorCode MatCompositeMerge(Mat mat) 599 { 600 Mat_Composite *shell = (Mat_Composite*)mat->data; 601 Mat_CompositeLink next = shell->head, prev = shell->tail; 602 PetscErrorCode ierr; 603 Mat tmat,newmat; 604 Vec left,right; 605 PetscScalar scale; 606 607 PetscFunctionBegin; 608 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 609 610 PetscFunctionBegin; 611 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 612 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 613 while ((next = next->next)) { 614 ierr = MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 615 } 616 } else { 617 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 618 while ((prev = prev->prev)) { 619 ierr = MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 620 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 621 tmat = newmat; 622 } 623 } 624 625 scale = shell->scale; 626 if ((left = shell->left)) {ierr = PetscObjectReference((PetscObject)left);CHKERRQ(ierr);} 627 if ((right = shell->right)) {ierr = PetscObjectReference((PetscObject)right);CHKERRQ(ierr);} 628 629 ierr = MatHeaderReplace(mat,tmat);CHKERRQ(ierr); 630 631 ierr = MatDiagonalScale(mat,left,right);CHKERRQ(ierr); 632 ierr = MatScale(mat,scale);CHKERRQ(ierr); 633 ierr = VecDestroy(&left);CHKERRQ(ierr); 634 ierr = VecDestroy(&right);CHKERRQ(ierr); 635 PetscFunctionReturn(0); 636 } 637