1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: matrix.c,v 1.333 1999/04/30 13:30:29 balay Exp balay $"; 3 #endif 4 5 /* 6 This is where the abstract matrix operations are defined 7 */ 8 9 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 10 #include "src/vec/vecimpl.h" 11 12 #undef __FUNC__ 13 #define __FUNC__ "MatGetRow" 14 /*@C 15 MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow() 16 for each row that you get to ensure that your application does 17 not bleed memory. 18 19 Not Collective 20 21 Input Parameters: 22 + mat - the matrix 23 - row - the row to get 24 25 Output Parameters: 26 + ncols - the number of nonzeros in the row 27 . cols - if nonzero, the column numbers 28 - vals - if nonzero, the values 29 30 Notes: 31 This routine is provided for people who need to have direct access 32 to the structure of a matrix. We hope that we provide enough 33 high-level matrix routines that few users will need it. 34 35 MatGetRow() always returns 0-based column indices, regardless of 36 whether the internal representation is 0-based (default) or 1-based. 37 38 For better efficiency, set cols and/or vals to PETSC_NULL if you do 39 not wish to extract these quantities. 40 41 The user can only examine the values extracted with MatGetRow(); 42 the values cannot be altered. To change the matrix entries, one 43 must use MatSetValues(). 44 45 You can only have one call to MatGetRow() outstanding for a particular 46 matrix at a time. 47 48 Fortran Notes: 49 The calling sequence from Fortran is 50 .vb 51 MatGetRow(matrix,row,ncols,cols,values,ierr) 52 Mat matrix (input) 53 integer row (input) 54 integer ncols (output) 55 integer cols(maxcols) (output) 56 double precision (or double complex) values(maxcols) output 57 .ve 58 where maxcols >= maximum nonzeros in any row of the matrix. 59 60 Caution: 61 Do not try to change the contents of the output arrays (cols and vals). 62 In some cases, this may corrupt the matrix. 63 64 Level: advanced 65 66 .keywords: matrix, row, get, extract 67 68 .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubmatrices(), MatGetDiagonal() 69 @*/ 70 int MatGetRow(Mat mat,int row,int *ncols,int **cols,Scalar **vals) 71 { 72 int ierr; 73 74 PetscFunctionBegin; 75 PetscValidHeaderSpecific(mat,MAT_COOKIE); 76 PetscValidIntPointer(ncols); 77 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 78 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 79 if (!mat->ops->getrow) SETERRQ(PETSC_ERR_SUP,0,""); 80 PLogEventBegin(MAT_GetRow,mat,0,0,0); 81 ierr = (*mat->ops->getrow)(mat,row,ncols,cols,vals); CHKERRQ(ierr); 82 PLogEventEnd(MAT_GetRow,mat,0,0,0); 83 PetscFunctionReturn(0); 84 } 85 86 #undef __FUNC__ 87 #define __FUNC__ "MatRestoreRow" 88 /*@C 89 MatRestoreRow - Frees any temporary space allocated by MatGetRow(). 90 91 Not Collective 92 93 Input Parameters: 94 + mat - the matrix 95 . row - the row to get 96 . ncols, cols - the number of nonzeros and their columns 97 - vals - if nonzero the column values 98 99 Notes: 100 This routine should be called after you have finished examining the entries. 101 102 Fortran Notes: 103 The calling sequence from Fortran is 104 .vb 105 MatRestoreRow(matrix,row,ncols,cols,values,ierr) 106 Mat matrix (input) 107 integer row (input) 108 integer ncols (output) 109 integer cols(maxcols) (output) 110 double precision (or double complex) values(maxcols) output 111 .ve 112 Where maxcols >= maximum nonzeros in any row of the matrix. 113 114 In Fortran MatRestoreRow() MUST be called after MatGetRow() 115 before another call to MatGetRow() can be made. 116 117 Level: advanced 118 119 .keywords: matrix, row, restore 120 121 .seealso: MatGetRow() 122 @*/ 123 int MatRestoreRow(Mat mat,int row,int *ncols,int **cols,Scalar **vals) 124 { 125 int ierr; 126 127 PetscFunctionBegin; 128 PetscValidHeaderSpecific(mat,MAT_COOKIE); 129 PetscValidIntPointer(ncols); 130 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 131 if (!mat->ops->restorerow) PetscFunctionReturn(0); 132 ierr = (*mat->ops->restorerow)(mat,row,ncols,cols,vals);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNC__ 137 #define __FUNC__ "MatView" 138 /*@C 139 MatView - Visualizes a matrix object. 140 141 Collective on Mat unless Viewer is VIEWER_STDOUT_SELF 142 143 Input Parameters: 144 + mat - the matrix 145 - ptr - visualization context 146 147 Notes: 148 The available visualization contexts include 149 + VIEWER_STDOUT_SELF - standard output (default) 150 . VIEWER_STDOUT_WORLD - synchronized standard 151 output where only the first processor opens 152 the file. All other processors send their 153 data to the first processor to print. 154 - VIEWER_DRAW_WORLD - graphical display of nonzero structure 155 156 The user can open alternative visualization contexts with 157 + ViewerASCIIOpen() - Outputs matrix to a specified file 158 . ViewerBinaryOpen() - Outputs matrix in binary to a 159 specified file; corresponding input uses MatLoad() 160 . ViewerDrawOpen() - Outputs nonzero matrix structure to 161 an X window display 162 - ViewerSocketOpen() - Outputs matrix to Socket viewer. 163 Currently only the sequential dense and AIJ 164 matrix types support the Socket viewer. 165 166 The user can call ViewerSetFormat() to specify the output 167 format of ASCII printed objects (when using VIEWER_STDOUT_SELF, 168 VIEWER_STDOUT_WORLD and ViewerASCIIOpen). Available formats include 169 + VIEWER_FORMAT_ASCII_DEFAULT - default, prints matrix contents 170 . VIEWER_FORMAT_ASCII_MATLAB - prints matrix contents in Matlab format 171 . VIEWER_FORMAT_ASCII_DENSE - prints entire matrix including zeros 172 . VIEWER_FORMAT_ASCII_COMMON - prints matrix contents, using a sparse 173 format common among all matrix types 174 . VIEWER_FORMAT_ASCII_IMPL - prints matrix contents, using an implementation-specific 175 format (which is in many cases the same as the default) 176 . VIEWER_FORMAT_ASCII_INFO - prints basic information about the matrix 177 size and structure (not the matrix entries) 178 - VIEWER_FORMAT_ASCII_INFO_LONG - prints more detailed information about 179 the matrix structure 180 181 Level: beginner 182 183 .keywords: matrix, view, visualize, output, print, write, draw 184 185 .seealso: ViewerSetFormat(), ViewerASCIIOpen(), ViewerDrawOpen(), 186 ViewerSocketOpen(), ViewerBinaryOpen(), MatLoad() 187 @*/ 188 int MatView(Mat mat,Viewer viewer) 189 { 190 int format, ierr, rows, cols; 191 char *cstr; 192 ViewerType vtype; 193 194 PetscFunctionBegin; 195 PetscValidHeaderSpecific(mat,MAT_COOKIE); 196 if (viewer) PetscValidHeaderSpecific(viewer,VIEWER_COOKIE); 197 198 if (!viewer) { 199 viewer = VIEWER_STDOUT_SELF; 200 } 201 202 ierr = ViewerGetType(viewer,&vtype); 203 if (PetscTypeCompare(vtype,ASCII_VIEWER)) { 204 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 205 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 206 ierr = ViewerASCIIPrintf(viewer,"Matrix Object:\n");CHKERRQ(ierr); 207 ierr = MatGetType(mat,PETSC_NULL,&cstr); CHKERRQ(ierr); 208 ierr = MatGetSize(mat,&rows,&cols); CHKERRQ(ierr); 209 ierr = ViewerASCIIPrintf(viewer," type=%s, rows=%d, cols=%d\n",cstr,rows,cols);CHKERRQ(ierr); 210 if (mat->ops->getinfo) { 211 MatInfo info; 212 ierr = MatGetInfo(mat,MAT_GLOBAL_SUM,&info); CHKERRQ(ierr); 213 ierr = ViewerASCIIPrintf(viewer," total: nonzeros=%d, allocated nonzeros=%d\n", 214 (int)info.nz_used,(int)info.nz_allocated);CHKERRQ(ierr); 215 } 216 } 217 } 218 if (mat->ops->view) { 219 ierr = ViewerASCIIPushTab(viewer);CHKERRQ(ierr); 220 ierr = (*mat->ops->view)(mat,viewer); CHKERRQ(ierr); 221 ierr = ViewerASCIIPopTab(viewer);CHKERRQ(ierr); 222 } 223 PetscFunctionReturn(0); 224 } 225 226 #undef __FUNC__ 227 #define __FUNC__ "MatScaleSystem" 228 /*@C 229 MatScaleSystem - Scale a vector solution and right hand side to 230 match the scaling of a scaled matrix. 231 232 Collective on Mat 233 234 Input Parameter: 235 + mat - the matrix 236 . x - solution vector (or PETSC_NULL) 237 + b - right hand side vector (or PETSC_NULL) 238 239 240 Notes: 241 For AIJ, BAIJ, and BDiag matrix formats, the matrices are not 242 internally scaled, so this does nothing. For MPIROWBS it 243 permutes and diagonally scales. 244 245 The SLES methods automatically call this routine when required 246 (via PCPreSolve()) so it is rarely used directly. 247 248 Level: Developer 249 250 .keywords: matrix, scale 251 252 .seealso: MatUseScaledForm(), MatUnScaleSystem() 253 @*/ 254 int MatScaleSystem(Mat mat,Vec x,Vec b) 255 { 256 int ierr; 257 258 PetscFunctionBegin; 259 PetscValidHeaderSpecific(mat,MAT_COOKIE); 260 if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);} 261 if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);} 262 if (mat->ops->scalesystem) { 263 ierr = (*mat->ops->scalesystem)(mat,x,b); CHKERRQ(ierr); 264 } 265 PetscFunctionReturn(0); 266 } 267 268 #undef __FUNC__ 269 #define __FUNC__ "MatUnScaleSystem" 270 /*@C 271 MatUnScaleSystem - Unscales a vector solution and right hand side to 272 match the original scaling of a scaled matrix. 273 274 Collective on Mat 275 276 Input Parameter: 277 + mat - the matrix 278 . x - solution vector (or PETSC_NULL) 279 + b - right hand side vector (or PETSC_NULL) 280 281 282 Notes: 283 For AIJ, BAIJ, and BDiag matrix formats, the matrices are not 284 internally scaled, so this does nothing. For MPIROWBS it 285 permutes and diagonally scales. 286 287 The SLES methods automatically call this routine when required 288 (via PCPreSolve()) so it is rarely used directly. 289 290 Level: Developer 291 292 .keywords: matrix, scale 293 294 .seealso: MatUseScaledForm(), MatScaleSystem() 295 @*/ 296 int MatUnScaleSystem(Mat mat,Vec x,Vec b) 297 { 298 int ierr; 299 300 PetscFunctionBegin; 301 PetscValidHeaderSpecific(mat,MAT_COOKIE); 302 if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);} 303 if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);} 304 if (mat->ops->unscalesystem) { 305 ierr = (*mat->ops->unscalesystem)(mat,x,b); CHKERRQ(ierr); 306 } 307 PetscFunctionReturn(0); 308 } 309 310 #undef __FUNC__ 311 #define __FUNC__ "MatUseScaledForm" 312 /*@C 313 MatUseScaledForm - For matrix storage formats that scale the 314 matrix (for example MPIRowBS matrices are diagonally scaled on 315 assembly) indicates matrix operations (MatMult() etc) are 316 applied using the scaled matrix. 317 318 Collective on Mat 319 320 Input Parameter: 321 + mat - the matrix 322 - scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for 323 applying the original matrix 324 325 Notes: 326 For scaled matrix formats, applying the original, unscaled matrix 327 will be slightly more expensive 328 329 Level: Developer 330 331 .keywords: matrix, scale 332 333 .seealso: MatScaleSystem(), MatUnScaleSystem() 334 @*/ 335 int MatUseScaledForm(Mat mat,PetscTruth scaled) 336 { 337 int ierr; 338 339 PetscFunctionBegin; 340 PetscValidHeaderSpecific(mat,MAT_COOKIE); 341 if (mat->ops->usescaledform) { 342 ierr = (*mat->ops->usescaledform)(mat,scaled); CHKERRQ(ierr); 343 } 344 PetscFunctionReturn(0); 345 } 346 347 #undef __FUNC__ 348 #define __FUNC__ "MatDestroy" 349 /*@C 350 MatDestroy - Frees space taken by a matrix. 351 352 Collective on Mat 353 354 Input Parameter: 355 . mat - the matrix 356 357 Level: beginner 358 359 .keywords: matrix, destroy 360 @*/ 361 int MatDestroy(Mat mat) 362 { 363 int ierr; 364 365 PetscFunctionBegin; 366 PetscValidHeaderSpecific(mat,MAT_COOKIE); 367 ierr = (*mat->ops->destroy)(mat); CHKERRQ(ierr); 368 PetscFunctionReturn(0); 369 } 370 371 #undef __FUNC__ 372 #define __FUNC__ "MatValid" 373 /*@ 374 MatValid - Checks whether a matrix object is valid. 375 376 Collective on Mat 377 378 Input Parameter: 379 . m - the matrix to check 380 381 Output Parameter: 382 flg - flag indicating matrix status, either 383 PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise. 384 385 Level: developer 386 387 .keywords: matrix, valid 388 @*/ 389 int MatValid(Mat m,PetscTruth *flg) 390 { 391 PetscFunctionBegin; 392 PetscValidIntPointer(flg); 393 if (!m) *flg = PETSC_FALSE; 394 else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE; 395 else *flg = PETSC_TRUE; 396 PetscFunctionReturn(0); 397 } 398 399 #undef __FUNC__ 400 #define __FUNC__ "MatSetValues" 401 /*@ 402 MatSetValues - Inserts or adds a block of values into a matrix. 403 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 404 MUST be called after all calls to MatSetValues() have been completed. 405 406 Not Collective 407 408 Input Parameters: 409 + mat - the matrix 410 . v - a logically two-dimensional array of values 411 . m, idxm - the number of rows and their global indices 412 . n, idxn - the number of columns and their global indices 413 - addv - either ADD_VALUES or INSERT_VALUES, where 414 ADD_VALUES adds values to any existing entries, and 415 INSERT_VALUES replaces existing entries with new values 416 417 Notes: 418 By default the values, v, are row-oriented and unsorted. 419 See MatSetOption() for other options. 420 421 Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES 422 options cannot be mixed without intervening calls to the assembly 423 routines. 424 425 MatSetValues() uses 0-based row and column numbers in Fortran 426 as well as in C. 427 428 Negative indices may be passed in idxm and idxn, these rows and columns are 429 simply ignored. This allows easily inserting element stiffness matrices 430 with homogeneous Dirchlet boundary conditions that you don't want represented 431 in the matrix. 432 433 Efficiency Alert: 434 The routine MatSetValuesBlocked() may offer much better efficiency 435 for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ). 436 437 Level: beginner 438 439 .keywords: matrix, insert, add, set, values 440 441 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal() 442 @*/ 443 int MatSetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv) 444 { 445 int ierr; 446 447 PetscFunctionBegin; 448 if (!m || !n) PetscFunctionReturn(0); /* no values to insert */ 449 PetscValidHeaderSpecific(mat,MAT_COOKIE); 450 PetscValidIntPointer(idxm); 451 PetscValidIntPointer(idxn); 452 PetscValidScalarPointer(v); 453 if (mat->insertmode == NOT_SET_VALUES) { 454 mat->insertmode = addv; 455 } 456 #if defined(USE_PETSC_BOPT_g) 457 else if (mat->insertmode != addv) { 458 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values"); 459 } 460 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 461 #endif 462 463 if (mat->assembled) { 464 mat->was_assembled = PETSC_TRUE; 465 mat->assembled = PETSC_FALSE; 466 } 467 PLogEventBegin(MAT_SetValues,mat,0,0,0); 468 if (!mat->ops->setvalues) SETERRQ(PETSC_ERR_SUP,1,"Not supported for this matrix type"); 469 ierr = (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr); 470 PLogEventEnd(MAT_SetValues,mat,0,0,0); 471 PetscFunctionReturn(0); 472 } 473 474 #undef __FUNC__ 475 #define __FUNC__ "MatSetValuesBlocked" 476 /*@ 477 MatSetValuesBlocked - Inserts or adds a block of values into a matrix. 478 479 Not Collective 480 481 Input Parameters: 482 + mat - the matrix 483 . v - a logically two-dimensional array of values 484 . m, idxm - the number of block rows and their global block indices 485 . n, idxn - the number of block columns and their global block indices 486 - addv - either ADD_VALUES or INSERT_VALUES, where 487 ADD_VALUES adds values to any existing entries, and 488 INSERT_VALUES replaces existing entries with new values 489 490 Notes: 491 By default the values, v, are row-oriented and unsorted. So the layout of 492 v is the same as for MatSetValues(). See MatSetOption() for other options. 493 494 Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES 495 options cannot be mixed without intervening calls to the assembly 496 routines. 497 498 MatSetValuesBlocked() uses 0-based row and column numbers in Fortran 499 as well as in C. 500 501 Negative indices may be passed in idxm and idxn, these rows and columns are 502 simply ignored. This allows easily inserting element stiffness matrices 503 with homogeneous Dirchlet boundary conditions that you don't want represented 504 in the matrix. 505 506 Each time an entry is set within a sparse matrix via MatSetValues(), 507 internal searching must be done to determine where to place the the 508 data in the matrix storage space. By instead inserting blocks of 509 entries via MatSetValuesBlocked(), the overhead of matrix assembly is 510 reduced. 511 512 Restrictions: 513 MatSetValuesBlocked() is currently supported only for the block AIJ 514 matrix format (MATSEQBAIJ and MATMPIBAIJ, which are created via 515 MatCreateSeqBAIJ() and MatCreateMPIBAIJ()). 516 517 Level: intermediate 518 519 .keywords: matrix, insert, add, set, values 520 521 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal() 522 @*/ 523 int MatSetValuesBlocked(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv) 524 { 525 int ierr; 526 527 PetscFunctionBegin; 528 if (!m || !n) PetscFunctionReturn(0); /* no values to insert */ 529 PetscValidHeaderSpecific(mat,MAT_COOKIE); 530 PetscValidIntPointer(idxm); 531 PetscValidIntPointer(idxn); 532 PetscValidScalarPointer(v); 533 if (mat->insertmode == NOT_SET_VALUES) { 534 mat->insertmode = addv; 535 } 536 #if defined(USE_PETSC_BOPT_g) 537 else if (mat->insertmode != addv) { 538 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values"); 539 } 540 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 541 #endif 542 543 if (mat->assembled) { 544 mat->was_assembled = PETSC_TRUE; 545 mat->assembled = PETSC_FALSE; 546 } 547 PLogEventBegin(MAT_SetValues,mat,0,0,0); 548 if (!mat->ops->setvaluesblocked) SETERRQ(PETSC_ERR_SUP,1,"Not supported for this matrix type"); 549 ierr = (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr); 550 PLogEventEnd(MAT_SetValues,mat,0,0,0); 551 PetscFunctionReturn(0); 552 } 553 554 /*MC 555 MatSetValue - Set a single entry into a matrix. 556 557 Synopsis: 558 void MatSetValue(Mat m,int row,int col,Scalar value,InsertMode mode); 559 560 Not collective 561 562 Input Parameters: 563 + m - the matrix 564 . row - the row location of the entry 565 . col - the column location of the entry 566 . value - the value to insert 567 - mode - either INSERT_VALUES or ADD_VALUES 568 569 Notes: 570 For efficiency one should use MatSetValues() and set several or many 571 values simultaneously if possible. 572 573 Note that VecSetValue() does NOT return an error code (since this 574 is checked internally). 575 576 Level: beginner 577 578 .seealso: MatSetValues() 579 M*/ 580 581 #undef __FUNC__ 582 #define __FUNC__ "MatGetValues" 583 /*@ 584 MatGetValues - Gets a block of values from a matrix. 585 586 Not Collective; currently only returns a local block 587 588 Input Parameters: 589 + mat - the matrix 590 . v - a logically two-dimensional array for storing the values 591 . m, idxm - the number of rows and their global indices 592 - n, idxn - the number of columns and their global indices 593 594 Notes: 595 The user must allocate space (m*n Scalars) for the values, v. 596 The values, v, are then returned in a row-oriented format, 597 analogous to that used by default in MatSetValues(). 598 599 MatGetValues() uses 0-based row and column numbers in 600 Fortran as well as in C. 601 602 MatGetValues() requires that the matrix has been assembled 603 with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to 604 MatSetValues() and MatGetValues() CANNOT be made in succession 605 without intermediate matrix assembly. 606 607 Level: advanced 608 609 .keywords: matrix, get, values 610 611 .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues() 612 @*/ 613 int MatGetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v) 614 { 615 int ierr; 616 617 PetscFunctionBegin; 618 PetscValidHeaderSpecific(mat,MAT_COOKIE); 619 PetscValidIntPointer(idxm); 620 PetscValidIntPointer(idxn); 621 PetscValidScalarPointer(v); 622 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 623 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 624 if (!mat->ops->getvalues) SETERRQ(PETSC_ERR_SUP,0,""); 625 626 PLogEventBegin(MAT_GetValues,mat,0,0,0); 627 ierr = (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v); CHKERRQ(ierr); 628 PLogEventEnd(MAT_GetValues,mat,0,0,0); 629 PetscFunctionReturn(0); 630 } 631 632 #undef __FUNC__ 633 #define __FUNC__ "MatSetLocalToGlobalMapping" 634 /*@ 635 MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by 636 the routine MatSetValuesLocal() to allow users to insert matrix entries 637 using a local (per-processor) numbering. 638 639 Not Collective 640 641 Input Parameters: 642 + x - the matrix 643 - mapping - mapping created with ISLocalToGlobalMappingCreate() 644 or ISLocalToGlobalMappingCreateIS() 645 646 Level: intermediate 647 648 .keywords: matrix, set, values, local, global, mapping 649 650 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal() 651 @*/ 652 int MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping) 653 { 654 PetscFunctionBegin; 655 PetscValidHeaderSpecific(x,MAT_COOKIE); 656 PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE); 657 658 if (x->mapping) { 659 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Mapping already set for matrix"); 660 } 661 662 x->mapping = mapping; 663 PetscObjectReference((PetscObject)mapping); 664 PetscFunctionReturn(0); 665 } 666 667 #undef __FUNC__ 668 #define __FUNC__ "MatSetLocalToGlobalMappingBlocked" 669 /*@ 670 MatSetLocalToGlobalMappingBlocked - Sets a local-to-global numbering for use 671 by the routine MatSetValuesBlockedLocal() to allow users to insert matrix 672 entries using a local (per-processor) numbering. 673 674 Not Collective 675 676 Input Parameters: 677 + x - the matrix 678 - mapping - mapping created with ISLocalToGlobalMappingCreate() or 679 ISLocalToGlobalMappingCreateIS() 680 681 Level: intermediate 682 683 .keywords: matrix, set, values, local ordering 684 685 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(), 686 MatSetValuesBlocked(), MatSetValuesLocal() 687 @*/ 688 int MatSetLocalToGlobalMappingBlocked(Mat x,ISLocalToGlobalMapping mapping) 689 { 690 PetscFunctionBegin; 691 PetscValidHeaderSpecific(x,MAT_COOKIE); 692 PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE); 693 694 if (x->bmapping) { 695 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Mapping already set for matrix"); 696 } 697 698 x->bmapping = mapping; 699 PetscObjectReference((PetscObject)mapping); 700 PetscFunctionReturn(0); 701 } 702 703 #undef __FUNC__ 704 #define __FUNC__ "MatSetValuesLocal" 705 /*@ 706 MatSetValuesLocal - Inserts or adds values into certain locations of a matrix, 707 using a local ordering of the nodes. 708 709 Not Collective 710 711 Input Parameters: 712 + x - the matrix 713 . nrow, irow - number of rows and their local indices 714 . ncol, icol - number of columns and their local indices 715 . y - a logically two-dimensional array of values 716 - addv - either INSERT_VALUES or ADD_VALUES, where 717 ADD_VALUES adds values to any existing entries, and 718 INSERT_VALUES replaces existing entries with new values 719 720 Notes: 721 Before calling MatSetValuesLocal(), the user must first set the 722 local-to-global mapping by calling MatSetLocalToGlobalMapping(). 723 724 Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES 725 options cannot be mixed without intervening calls to the assembly 726 routines. 727 728 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 729 MUST be called after all calls to MatSetValuesLocal() have been completed. 730 731 Level: intermediate 732 733 .keywords: matrix, set, values, local ordering 734 735 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping() 736 @*/ 737 int MatSetValuesLocal(Mat mat,int nrow,int *irow,int ncol, int *icol,Scalar *y,InsertMode addv) 738 { 739 int ierr,irowm[2048],icolm[2048]; 740 741 PetscFunctionBegin; 742 PetscValidHeaderSpecific(mat,MAT_COOKIE); 743 PetscValidIntPointer(irow); 744 PetscValidIntPointer(icol); 745 PetscValidScalarPointer(y); 746 747 if (mat->insertmode == NOT_SET_VALUES) { 748 mat->insertmode = addv; 749 } 750 #if defined(USE_PETSC_BOPT_g) 751 else if (mat->insertmode != addv) { 752 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values"); 753 } 754 if (!mat->mapping) { 755 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMapping()"); 756 } 757 if (nrow > 2048 || ncol > 2048) { 758 SETERRQ2(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol); 759 } 760 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 761 #endif 762 763 if (mat->assembled) { 764 mat->was_assembled = PETSC_TRUE; 765 mat->assembled = PETSC_FALSE; 766 } 767 PLogEventBegin(MAT_SetValues,mat,0,0,0); 768 ierr = ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm); CHKERRQ(ierr); 769 ierr = ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);CHKERRQ(ierr); 770 ierr = (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr); 771 PLogEventEnd(MAT_SetValues,mat,0,0,0); 772 PetscFunctionReturn(0); 773 } 774 775 #undef __FUNC__ 776 #define __FUNC__ "MatSetValuesBlockedLocal" 777 /*@ 778 MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix, 779 using a local ordering of the nodes a block at a time. 780 781 Not Collective 782 783 Input Parameters: 784 + x - the matrix 785 . nrow, irow - number of rows and their local indices 786 . ncol, icol - number of columns and their local indices 787 . y - a logically two-dimensional array of values 788 - addv - either INSERT_VALUES or ADD_VALUES, where 789 ADD_VALUES adds values to any existing entries, and 790 INSERT_VALUES replaces existing entries with new values 791 792 Notes: 793 Before calling MatSetValuesBlockedLocal(), the user must first set the 794 local-to-global mapping by calling MatSetLocalToGlobalMappingBlocked(), 795 where the mapping MUST be set for matrix blocks, not for matrix elements. 796 797 Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES 798 options cannot be mixed without intervening calls to the assembly 799 routines. 800 801 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 802 MUST be called after all calls to MatSetValuesBlockedLocal() have been completed. 803 804 Level: intermediate 805 806 .keywords: matrix, set, values, blocked, local 807 808 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlocked(), MatSetValuesBlocked() 809 @*/ 810 int MatSetValuesBlockedLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv) 811 { 812 int ierr,irowm[2048],icolm[2048]; 813 814 PetscFunctionBegin; 815 PetscValidHeaderSpecific(mat,MAT_COOKIE); 816 PetscValidIntPointer(irow); 817 PetscValidIntPointer(icol); 818 PetscValidScalarPointer(y); 819 if (mat->insertmode == NOT_SET_VALUES) { 820 mat->insertmode = addv; 821 } 822 #if defined(USE_PETSC_BOPT_g) 823 else if (mat->insertmode != addv) { 824 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Cannot mix add values and insert values"); 825 } 826 if (!mat->bmapping) { 827 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMappingBlocked()"); 828 } 829 if (nrow > 2048 || ncol > 2048) { 830 SETERRQ2(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol); 831 } 832 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 833 #endif 834 835 if (mat->assembled) { 836 mat->was_assembled = PETSC_TRUE; 837 mat->assembled = PETSC_FALSE; 838 } 839 PLogEventBegin(MAT_SetValues,mat,0,0,0); 840 ierr = ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm); CHKERRQ(ierr); 841 ierr = ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm); CHKERRQ(ierr); 842 ierr = (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr); 843 PLogEventEnd(MAT_SetValues,mat,0,0,0); 844 PetscFunctionReturn(0); 845 } 846 847 /* --------------------------------------------------------*/ 848 #undef __FUNC__ 849 #define __FUNC__ "MatMult" 850 /*@ 851 MatMult - Computes the matrix-vector product, y = Ax. 852 853 Collective on Mat and Vec 854 855 Input Parameters: 856 + mat - the matrix 857 - x - the vector to be multilplied 858 859 Output Parameters: 860 . y - the result 861 862 Notes: 863 The vectors x and y cannot be the same. I.e., one cannot 864 call MatMult(A,y,y). 865 866 Level: beginner 867 868 .keywords: matrix, multiply, matrix-vector product 869 870 .seealso: MatMultTrans(), MatMultAdd(), MatMultTransAdd() 871 @*/ 872 int MatMult(Mat mat,Vec x,Vec y) 873 { 874 int ierr; 875 876 PetscFunctionBegin; 877 PetscValidHeaderSpecific(mat,MAT_COOKIE); 878 PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 879 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 880 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 881 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors"); 882 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 883 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 884 if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: local dim %d %d",mat->m,y->n); 885 886 PLogEventBegin(MAT_Mult,mat,x,y,0); 887 ierr = (*mat->ops->mult)(mat,x,y); CHKERRQ(ierr); 888 PLogEventEnd(MAT_Mult,mat,x,y,0); 889 890 PetscFunctionReturn(0); 891 } 892 893 #undef __FUNC__ 894 #define __FUNC__ "MatMultTrans" 895 /*@ 896 MatMultTrans - Computes matrix transpose times a vector. 897 898 Collective on Mat and Vec 899 900 Input Parameters: 901 + mat - the matrix 902 - x - the vector to be multilplied 903 904 Output Parameters: 905 . y - the result 906 907 Notes: 908 The vectors x and y cannot be the same. I.e., one cannot 909 call MatMultTrans(A,y,y). 910 911 Level: beginner 912 913 .keywords: matrix, multiply, matrix-vector product, transpose 914 915 .seealso: MatMult(), MatMultAdd(), MatMultTransAdd() 916 @*/ 917 int MatMultTrans(Mat mat,Vec x,Vec y) 918 { 919 int ierr; 920 921 PetscFunctionBegin; 922 PetscValidHeaderSpecific(mat,MAT_COOKIE); 923 PetscValidHeaderSpecific(x,VEC_COOKIE); PetscValidHeaderSpecific(y,VEC_COOKIE); 924 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 925 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 926 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors"); 927 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 928 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->N,y->N); 929 PLogEventBegin(MAT_MultTrans,mat,x,y,0); 930 ierr = (*mat->ops->multtrans)(mat,x,y); CHKERRQ(ierr); 931 PLogEventEnd(MAT_MultTrans,mat,x,y,0); 932 PetscFunctionReturn(0); 933 } 934 935 #undef __FUNC__ 936 #define __FUNC__ "MatMultAdd" 937 /*@ 938 MatMultAdd - Computes v3 = v2 + A * v1. 939 940 Collective on Mat and Vec 941 942 Input Parameters: 943 + mat - the matrix 944 - v1, v2 - the vectors 945 946 Output Parameters: 947 . v3 - the result 948 949 Notes: 950 The vectors v1 and v3 cannot be the same. I.e., one cannot 951 call MatMultAdd(A,v1,v2,v1). 952 953 Level: beginner 954 955 .keywords: matrix, multiply, matrix-vector product, add 956 957 .seealso: MatMultTrans(), MatMult(), MatMultTransAdd() 958 @*/ 959 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3) 960 { 961 int ierr; 962 963 PetscFunctionBegin; 964 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v1,VEC_COOKIE); 965 PetscValidHeaderSpecific(v2,VEC_COOKIE); PetscValidHeaderSpecific(v3,VEC_COOKIE); 966 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 967 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 968 if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N); 969 if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N); 970 if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N); 971 if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n); 972 if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n); 973 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors"); 974 975 PLogEventBegin(MAT_MultAdd,mat,v1,v2,v3); 976 ierr = (*mat->ops->multadd)(mat,v1,v2,v3); CHKERRQ(ierr); 977 PLogEventEnd(MAT_MultAdd,mat,v1,v2,v3); 978 PetscFunctionReturn(0); 979 } 980 981 #undef __FUNC__ 982 #define __FUNC__ "MatMultTransAdd" 983 /*@ 984 MatMultTransAdd - Computes v3 = v2 + A' * v1. 985 986 Collective on Mat and Vec 987 988 Input Parameters: 989 + mat - the matrix 990 - v1, v2 - the vectors 991 992 Output Parameters: 993 . v3 - the result 994 995 Notes: 996 The vectors v1 and v3 cannot be the same. I.e., one cannot 997 call MatMultTransAdd(A,v1,v2,v1). 998 999 Level: beginner 1000 1001 .keywords: matrix, multiply, matrix-vector product, transpose, add 1002 1003 .seealso: MatMultTrans(), MatMultAdd(), MatMult() 1004 @*/ 1005 int MatMultTransAdd(Mat mat,Vec v1,Vec v2,Vec v3) 1006 { 1007 int ierr; 1008 1009 PetscFunctionBegin; 1010 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v1,VEC_COOKIE); 1011 PetscValidHeaderSpecific(v2,VEC_COOKIE);PetscValidHeaderSpecific(v3,VEC_COOKIE); 1012 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1013 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1014 if (!mat->ops->multtransadd) SETERRQ(PETSC_ERR_SUP,0,""); 1015 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors"); 1016 if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N); 1017 if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N); 1018 if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N); 1019 1020 PLogEventBegin(MAT_MultTransAdd,mat,v1,v2,v3); 1021 ierr = (*mat->ops->multtransadd)(mat,v1,v2,v3); CHKERRQ(ierr); 1022 PLogEventEnd(MAT_MultTransAdd,mat,v1,v2,v3); 1023 PetscFunctionReturn(0); 1024 } 1025 /* ------------------------------------------------------------*/ 1026 #undef __FUNC__ 1027 #define __FUNC__ "MatGetInfo" 1028 /*@C 1029 MatGetInfo - Returns information about matrix storage (number of 1030 nonzeros, memory, etc.). 1031 1032 Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used 1033 as the flag 1034 1035 Input Parameters: 1036 . mat - the matrix 1037 1038 Output Parameters: 1039 + flag - flag indicating the type of parameters to be returned 1040 (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors, 1041 MAT_GLOBAL_SUM - sum over all processors) 1042 - info - matrix information context 1043 1044 Notes: 1045 The MatInfo context contains a variety of matrix data, including 1046 number of nonzeros allocated and used, number of mallocs during 1047 matrix assembly, etc. Additional information for factored matrices 1048 is provided (such as the fill ratio, number of mallocs during 1049 factorization, etc.). Much of this info is printed to STDOUT 1050 when using the runtime options 1051 $ -log_info -mat_view_info 1052 1053 Example for C/C++ Users: 1054 See the file ${PETSC_DIR}/include/mat.h for a complete list of 1055 data within the MatInfo context. For example, 1056 .vb 1057 MatInfo info; 1058 Mat A; 1059 double mal, nz_a, nz_u; 1060 1061 MatGetInfo(A,MAT_LOCAL,&info); 1062 mal = info.mallocs; 1063 nz_a = info.nz_allocated; 1064 .ve 1065 1066 Example for Fortran Users: 1067 Fortran users should declare info as a double precision 1068 array of dimension MAT_INFO_SIZE, and then extract the parameters 1069 of interest. See the file ${PETSC_DIR}/include/finclude/mat.h 1070 a complete list of parameter names. 1071 .vb 1072 double precision info(MAT_INFO_SIZE) 1073 double precision mal, nz_a 1074 Mat A 1075 integer ierr 1076 1077 call MatGetInfo(A,MAT_LOCAL,info,ierr) 1078 mal = info(MAT_INFO_MALLOCS) 1079 nz_a = info(MAT_INFO_NZ_ALLOCATED) 1080 .ve 1081 1082 Level: intermediate 1083 1084 .keywords: matrix, get, info, storage, nonzeros, memory, fill 1085 @*/ 1086 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info) 1087 { 1088 int ierr; 1089 1090 PetscFunctionBegin; 1091 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1092 PetscValidPointer(info); 1093 if (!mat->ops->getinfo) SETERRQ(PETSC_ERR_SUP,0,""); 1094 ierr = (*mat->ops->getinfo)(mat,flag,info);CHKERRQ(ierr); 1095 PetscFunctionReturn(0); 1096 } 1097 1098 /* ----------------------------------------------------------*/ 1099 #undef __FUNC__ 1100 #define __FUNC__ "MatILUDTFactor" 1101 /*@ 1102 MatILUDTFactor - Performs a drop tolerance ILU factorization. 1103 1104 Collective on Mat 1105 1106 Input Parameters: 1107 + mat - the matrix 1108 . dt - the drop tolerance 1109 . maxnz - the maximum number of nonzeros per row allowed 1110 . row - row permutation 1111 - col - column permutation 1112 1113 Output Parameters: 1114 . fact - the factored matrix 1115 1116 Notes: 1117 Most users should employ the simplified SLES interface for linear solvers 1118 instead of working directly with matrix algebra routines such as this. 1119 See, e.g., SLESCreate(). 1120 1121 Level: developer 1122 1123 .keywords: matrix, factor, LU, in-place 1124 1125 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 1126 @*/ 1127 int MatILUDTFactor(Mat mat,double dt,int maxnz,IS row,IS col,Mat *fact) 1128 { 1129 int ierr; 1130 1131 PetscFunctionBegin; 1132 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1133 PetscValidPointer(fact); 1134 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1135 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1136 if (!mat->ops->iludtfactor) SETERRQ(PETSC_ERR_SUP,0,""); 1137 1138 PLogEventBegin(MAT_ILUFactor,mat,row,col,0); 1139 ierr = (*mat->ops->iludtfactor)(mat,dt,maxnz,row,col,fact); CHKERRQ(ierr); 1140 PLogEventEnd(MAT_ILUFactor,mat,row,col,0); 1141 1142 PetscFunctionReturn(0); 1143 } 1144 1145 #undef __FUNC__ 1146 #define __FUNC__ "MatLUFactor" 1147 /*@ 1148 MatLUFactor - Performs in-place LU factorization of matrix. 1149 1150 Collective on Mat 1151 1152 Input Parameters: 1153 + mat - the matrix 1154 . row - row permutation 1155 . col - column permutation 1156 - f - expected fill as ratio of original fill. 1157 1158 Notes: 1159 Most users should employ the simplified SLES interface for linear solvers 1160 instead of working directly with matrix algebra routines such as this. 1161 See, e.g., SLESCreate(). 1162 1163 Level: developer 1164 1165 .keywords: matrix, factor, LU, in-place 1166 1167 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 1168 @*/ 1169 int MatLUFactor(Mat mat,IS row,IS col,double f) 1170 { 1171 int ierr; 1172 1173 PetscFunctionBegin; 1174 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1175 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square"); 1176 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1177 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1178 if (!mat->ops->lufactor) SETERRQ(PETSC_ERR_SUP,0,""); 1179 1180 PLogEventBegin(MAT_LUFactor,mat,row,col,0); 1181 ierr = (*mat->ops->lufactor)(mat,row,col,f); CHKERRQ(ierr); 1182 PLogEventEnd(MAT_LUFactor,mat,row,col,0); 1183 PetscFunctionReturn(0); 1184 } 1185 1186 #undef __FUNC__ 1187 #define __FUNC__ "MatILUFactor" 1188 /*@ 1189 MatILUFactor - Performs in-place ILU factorization of matrix. 1190 1191 Collective on Mat 1192 1193 Input Parameters: 1194 + mat - the matrix 1195 . row - row permutation 1196 . col - column permutation 1197 - info - structure containing 1198 $ levels - number of levels of fill. 1199 $ expected fill - as ratio of original fill. 1200 $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices 1201 missing diagonal entries) 1202 1203 Notes: 1204 Probably really in-place only when level of fill is zero, otherwise allocates 1205 new space to store factored matrix and deletes previous memory. 1206 1207 Most users should employ the simplified SLES interface for linear solvers 1208 instead of working directly with matrix algebra routines such as this. 1209 See, e.g., SLESCreate(). 1210 1211 Level: developer 1212 1213 .keywords: matrix, factor, ILU, in-place 1214 1215 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 1216 @*/ 1217 int MatILUFactor(Mat mat,IS row,IS col,MatILUInfo *info) 1218 { 1219 int ierr; 1220 1221 PetscFunctionBegin; 1222 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1223 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square"); 1224 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1225 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1226 if (!mat->ops->ilufactor) SETERRQ(PETSC_ERR_SUP,0,""); 1227 1228 PLogEventBegin(MAT_ILUFactor,mat,row,col,0); 1229 ierr = (*mat->ops->ilufactor)(mat,row,col,info); CHKERRQ(ierr); 1230 PLogEventEnd(MAT_ILUFactor,mat,row,col,0); 1231 PetscFunctionReturn(0); 1232 } 1233 1234 #undef __FUNC__ 1235 #define __FUNC__ "MatLUFactorSymbolic" 1236 /*@ 1237 MatLUFactorSymbolic - Performs symbolic LU factorization of matrix. 1238 Call this routine before calling MatLUFactorNumeric(). 1239 1240 Collective on Mat 1241 1242 Input Parameters: 1243 + mat - the matrix 1244 . row, col - row and column permutations 1245 - f - expected fill as ratio of the original number of nonzeros, 1246 for example 3.0; choosing this parameter well can result in 1247 more efficient use of time and space. Run with the option -log_info 1248 to determine an optimal value to use 1249 1250 Output Parameter: 1251 . fact - new matrix that has been symbolically factored 1252 1253 Notes: 1254 See the users manual for additional information about 1255 choosing the fill factor for better efficiency. 1256 1257 Most users should employ the simplified SLES interface for linear solvers 1258 instead of working directly with matrix algebra routines such as this. 1259 See, e.g., SLESCreate(). 1260 1261 Level: developer 1262 1263 .keywords: matrix, factor, LU, symbolic, fill 1264 1265 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor() 1266 @*/ 1267 int MatLUFactorSymbolic(Mat mat,IS row,IS col,double f,Mat *fact) 1268 { 1269 int ierr; 1270 1271 PetscFunctionBegin; 1272 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1273 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square"); 1274 PetscValidPointer(fact); 1275 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1276 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1277 if (!mat->ops->lufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,""); 1278 1279 PLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0); 1280 ierr = (*mat->ops->lufactorsymbolic)(mat,row,col,f,fact); CHKERRQ(ierr); 1281 PLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0); 1282 PetscFunctionReturn(0); 1283 } 1284 1285 #undef __FUNC__ 1286 #define __FUNC__ "MatLUFactorNumeric" 1287 /*@ 1288 MatLUFactorNumeric - Performs numeric LU factorization of a matrix. 1289 Call this routine after first calling MatLUFactorSymbolic(). 1290 1291 Collective on Mat 1292 1293 Input Parameters: 1294 + mat - the matrix 1295 - row, col - row and column permutations 1296 1297 Output Parameters: 1298 . fact - symbolically factored matrix that must have been generated 1299 by MatLUFactorSymbolic() 1300 1301 Notes: 1302 See MatLUFactor() for in-place factorization. See 1303 MatCholeskyFactorNumeric() for the symmetric, positive definite case. 1304 1305 Most users should employ the simplified SLES interface for linear solvers 1306 instead of working directly with matrix algebra routines such as this. 1307 See, e.g., SLESCreate(). 1308 1309 Level: developer 1310 1311 .keywords: matrix, factor, LU, numeric 1312 1313 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor() 1314 @*/ 1315 int MatLUFactorNumeric(Mat mat,Mat *fact) 1316 { 1317 int ierr,flg; 1318 1319 PetscFunctionBegin; 1320 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1321 PetscValidPointer(fact); 1322 PetscValidHeaderSpecific(*fact,MAT_COOKIE); 1323 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1324 if (mat->M != (*fact)->M || mat->N != (*fact)->N) { 1325 SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d", 1326 mat->M,(*fact)->M,mat->N,(*fact)->N); 1327 } 1328 if (!(*fact)->ops->lufactornumeric) SETERRQ(PETSC_ERR_SUP,0,""); 1329 1330 PLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0); 1331 ierr = (*(*fact)->ops->lufactornumeric)(mat,fact); CHKERRQ(ierr); 1332 PLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0); 1333 ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg); CHKERRQ(ierr); 1334 if (flg) { 1335 ierr = OptionsHasName(0,"-mat_view_contour",&flg); CHKERRQ(ierr); 1336 if (flg) { 1337 ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr); 1338 } 1339 ierr = MatView(*fact,VIEWER_DRAW_(mat->comm)); CHKERRQ(ierr); 1340 ierr = ViewerFlush(VIEWER_DRAW_(mat->comm)); CHKERRQ(ierr); 1341 if (flg) { 1342 ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 1343 } 1344 } 1345 PetscFunctionReturn(0); 1346 } 1347 1348 #undef __FUNC__ 1349 #define __FUNC__ "MatCholeskyFactor" 1350 /*@ 1351 MatCholeskyFactor - Performs in-place Cholesky factorization of a 1352 symmetric matrix. 1353 1354 Collective on Mat 1355 1356 Input Parameters: 1357 + mat - the matrix 1358 . perm - row and column permutations 1359 - f - expected fill as ratio of original fill 1360 1361 Notes: 1362 See MatLUFactor() for the nonsymmetric case. See also 1363 MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric(). 1364 1365 Most users should employ the simplified SLES interface for linear solvers 1366 instead of working directly with matrix algebra routines such as this. 1367 See, e.g., SLESCreate(). 1368 1369 Level: developer 1370 1371 .keywords: matrix, factor, in-place, Cholesky 1372 1373 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric() 1374 @*/ 1375 int MatCholeskyFactor(Mat mat,IS perm,double f) 1376 { 1377 int ierr; 1378 1379 PetscFunctionBegin; 1380 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1381 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Matrix must be square"); 1382 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1383 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1384 if (!mat->ops->choleskyfactor) SETERRQ(PETSC_ERR_SUP,0,""); 1385 1386 PLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0); 1387 ierr = (*mat->ops->choleskyfactor)(mat,perm,f); CHKERRQ(ierr); 1388 PLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0); 1389 PetscFunctionReturn(0); 1390 } 1391 1392 #undef __FUNC__ 1393 #define __FUNC__ "MatCholeskyFactorSymbolic" 1394 /*@ 1395 MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization 1396 of a symmetric matrix. 1397 1398 Collective on Mat 1399 1400 Input Parameters: 1401 + mat - the matrix 1402 . perm - row and column permutations 1403 - f - expected fill as ratio of original 1404 1405 Output Parameter: 1406 . fact - the factored matrix 1407 1408 Notes: 1409 See MatLUFactorSymbolic() for the nonsymmetric case. See also 1410 MatCholeskyFactor() and MatCholeskyFactorNumeric(). 1411 1412 Most users should employ the simplified SLES interface for linear solvers 1413 instead of working directly with matrix algebra routines such as this. 1414 See, e.g., SLESCreate(). 1415 1416 Level: developer 1417 1418 .keywords: matrix, factor, factorization, symbolic, Cholesky 1419 1420 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric() 1421 @*/ 1422 int MatCholeskyFactorSymbolic(Mat mat,IS perm,double f,Mat *fact) 1423 { 1424 int ierr; 1425 1426 PetscFunctionBegin; 1427 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1428 PetscValidPointer(fact); 1429 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Matrix must be square"); 1430 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1431 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1432 if (!mat->ops->choleskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,""); 1433 1434 PLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0); 1435 ierr = (*mat->ops->choleskyfactorsymbolic)(mat,perm,f,fact); CHKERRQ(ierr); 1436 PLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0); 1437 PetscFunctionReturn(0); 1438 } 1439 1440 #undef __FUNC__ 1441 #define __FUNC__ "MatCholeskyFactorNumeric" 1442 /*@ 1443 MatCholeskyFactorNumeric - Performs numeric Cholesky factorization 1444 of a symmetric matrix. Call this routine after first calling 1445 MatCholeskyFactorSymbolic(). 1446 1447 Collective on Mat 1448 1449 Input Parameter: 1450 . mat - the initial matrix 1451 1452 Output Parameter: 1453 . fact - the factored matrix 1454 1455 Notes: 1456 Most users should employ the simplified SLES interface for linear solvers 1457 instead of working directly with matrix algebra routines such as this. 1458 See, e.g., SLESCreate(). 1459 1460 Level: developer 1461 1462 .keywords: matrix, factor, numeric, Cholesky 1463 1464 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric() 1465 @*/ 1466 int MatCholeskyFactorNumeric(Mat mat,Mat *fact) 1467 { 1468 int ierr; 1469 1470 PetscFunctionBegin; 1471 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1472 PetscValidPointer(fact); 1473 if (!mat->ops->choleskyfactornumeric) SETERRQ(PETSC_ERR_SUP,0,""); 1474 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1475 if (mat->M != (*fact)->M || mat->N != (*fact)->N) { 1476 SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d", 1477 mat->M,(*fact)->M,mat->N,(*fact)->N); 1478 } 1479 1480 PLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0); 1481 ierr = (*mat->ops->choleskyfactornumeric)(mat,fact); CHKERRQ(ierr); 1482 PLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0); 1483 PetscFunctionReturn(0); 1484 } 1485 1486 /* ----------------------------------------------------------------*/ 1487 #undef __FUNC__ 1488 #define __FUNC__ "MatSolve" 1489 /*@ 1490 MatSolve - Solves A x = b, given a factored matrix. 1491 1492 Collective on Mat and Vec 1493 1494 Input Parameters: 1495 + mat - the factored matrix 1496 - b - the right-hand-side vector 1497 1498 Output Parameter: 1499 . x - the result vector 1500 1501 Notes: 1502 The vectors b and x cannot be the same. I.e., one cannot 1503 call MatSolve(A,x,x). 1504 1505 Notes: 1506 Most users should employ the simplified SLES interface for linear solvers 1507 instead of working directly with matrix algebra routines such as this. 1508 See, e.g., SLESCreate(). 1509 1510 Level: developer 1511 1512 .keywords: matrix, linear system, solve, LU, Cholesky, triangular solve 1513 1514 .seealso: MatSolveAdd(), MatSolveTrans(), MatSolveTransAdd() 1515 @*/ 1516 int MatSolve(Mat mat,Vec b,Vec x) 1517 { 1518 int ierr; 1519 1520 PetscFunctionBegin; 1521 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1522 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1523 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1524 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1525 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1526 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1527 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1528 if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0); 1529 1530 if (!mat->ops->solve) SETERRQ(PETSC_ERR_SUP,0,""); 1531 PLogEventBegin(MAT_Solve,mat,b,x,0); 1532 ierr = (*mat->ops->solve)(mat,b,x); CHKERRQ(ierr); 1533 PLogEventEnd(MAT_Solve,mat,b,x,0); 1534 PetscFunctionReturn(0); 1535 } 1536 1537 #undef __FUNC__ 1538 #define __FUNC__ "MatForwardSolve" 1539 /* @ 1540 MatForwardSolve - Solves L x = b, given a factored matrix, A = LU. 1541 1542 Collective on Mat and Vec 1543 1544 Input Parameters: 1545 + mat - the factored matrix 1546 - b - the right-hand-side vector 1547 1548 Output Parameter: 1549 . x - the result vector 1550 1551 Notes: 1552 MatSolve() should be used for most applications, as it performs 1553 a forward solve followed by a backward solve. 1554 1555 The vectors b and x cannot be the same. I.e., one cannot 1556 call MatForwardSolve(A,x,x). 1557 1558 Most users should employ the simplified SLES interface for linear solvers 1559 instead of working directly with matrix algebra routines such as this. 1560 See, e.g., SLESCreate(). 1561 1562 Level: developer 1563 1564 .keywords: matrix, forward, LU, Cholesky, triangular solve 1565 1566 .seealso: MatSolve(), MatBackwardSolve() 1567 @ */ 1568 int MatForwardSolve(Mat mat,Vec b,Vec x) 1569 { 1570 int ierr; 1571 1572 PetscFunctionBegin; 1573 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1574 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1575 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1576 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1577 if (!mat->ops->forwardsolve) SETERRQ(PETSC_ERR_SUP,0,""); 1578 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1579 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1580 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1581 1582 PLogEventBegin(MAT_ForwardSolve,mat,b,x,0); 1583 ierr = (*mat->ops->forwardsolve)(mat,b,x); CHKERRQ(ierr); 1584 PLogEventEnd(MAT_ForwardSolve,mat,b,x,0); 1585 PetscFunctionReturn(0); 1586 } 1587 1588 #undef __FUNC__ 1589 #define __FUNC__ "MatBackwardSolve" 1590 /* @ 1591 MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU. 1592 1593 Collective on Mat and Vec 1594 1595 Input Parameters: 1596 + mat - the factored matrix 1597 - b - the right-hand-side vector 1598 1599 Output Parameter: 1600 . x - the result vector 1601 1602 Notes: 1603 MatSolve() should be used for most applications, as it performs 1604 a forward solve followed by a backward solve. 1605 1606 The vectors b and x cannot be the same. I.e., one cannot 1607 call MatBackwardSolve(A,x,x). 1608 1609 Most users should employ the simplified SLES interface for linear solvers 1610 instead of working directly with matrix algebra routines such as this. 1611 See, e.g., SLESCreate(). 1612 1613 Level: developer 1614 1615 .keywords: matrix, backward, LU, Cholesky, triangular solve 1616 1617 .seealso: MatSolve(), MatForwardSolve() 1618 @ */ 1619 int MatBackwardSolve(Mat mat,Vec b,Vec x) 1620 { 1621 int ierr; 1622 1623 PetscFunctionBegin; 1624 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1625 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1626 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1627 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1628 if (!mat->ops->backwardsolve) SETERRQ(PETSC_ERR_SUP,0,""); 1629 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1630 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1631 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1632 1633 PLogEventBegin(MAT_BackwardSolve,mat,b,x,0); 1634 ierr = (*mat->ops->backwardsolve)(mat,b,x); CHKERRQ(ierr); 1635 PLogEventEnd(MAT_BackwardSolve,mat,b,x,0); 1636 PetscFunctionReturn(0); 1637 } 1638 1639 #undef __FUNC__ 1640 #define __FUNC__ "MatSolveAdd" 1641 /*@ 1642 MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix. 1643 1644 Collective on Mat and Vec 1645 1646 Input Parameters: 1647 + mat - the factored matrix 1648 . b - the right-hand-side vector 1649 - y - the vector to be added to 1650 1651 Output Parameter: 1652 . x - the result vector 1653 1654 Notes: 1655 The vectors b and x cannot be the same. I.e., one cannot 1656 call MatSolveAdd(A,x,y,x). 1657 1658 Most users should employ the simplified SLES interface for linear solvers 1659 instead of working directly with matrix algebra routines such as this. 1660 See, e.g., SLESCreate(). 1661 1662 Level: developer 1663 1664 .keywords: matrix, linear system, solve, LU, Cholesky, add 1665 1666 .seealso: MatSolve(), MatSolveTrans(), MatSolveTransAdd() 1667 @*/ 1668 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x) 1669 { 1670 Scalar one = 1.0; 1671 Vec tmp; 1672 int ierr; 1673 1674 PetscFunctionBegin; 1675 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1676 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1677 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1678 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1679 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1680 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1681 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 1682 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1683 if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim %d %d",x->n,y->n); 1684 1685 PLogEventBegin(MAT_SolveAdd,mat,b,x,y); 1686 if (mat->ops->solveadd) { 1687 ierr = (*mat->ops->solveadd)(mat,b,y,x); CHKERRQ(ierr); 1688 } else { 1689 /* do the solve then the add manually */ 1690 if (x != y) { 1691 ierr = MatSolve(mat,b,x); CHKERRQ(ierr); 1692 ierr = VecAXPY(&one,y,x); CHKERRQ(ierr); 1693 } else { 1694 ierr = VecDuplicate(x,&tmp); CHKERRQ(ierr); 1695 PLogObjectParent(mat,tmp); 1696 ierr = VecCopy(x,tmp); CHKERRQ(ierr); 1697 ierr = MatSolve(mat,b,x); CHKERRQ(ierr); 1698 ierr = VecAXPY(&one,tmp,x); CHKERRQ(ierr); 1699 ierr = VecDestroy(tmp); CHKERRQ(ierr); 1700 } 1701 } 1702 PLogEventEnd(MAT_SolveAdd,mat,b,x,y); 1703 PetscFunctionReturn(0); 1704 } 1705 1706 #undef __FUNC__ 1707 #define __FUNC__ "MatSolveTrans" 1708 /*@ 1709 MatSolveTrans - Solves A' x = b, given a factored matrix. 1710 1711 Collective on Mat and Vec 1712 1713 Input Parameters: 1714 + mat - the factored matrix 1715 - b - the right-hand-side vector 1716 1717 Output Parameter: 1718 . x - the result vector 1719 1720 Notes: 1721 The vectors b and x cannot be the same. I.e., one cannot 1722 call MatSolveTrans(A,x,x). 1723 1724 Most users should employ the simplified SLES interface for linear solvers 1725 instead of working directly with matrix algebra routines such as this. 1726 See, e.g., SLESCreate(). 1727 1728 Level: developer 1729 1730 .keywords: matrix, linear system, solve, LU, Cholesky, transpose 1731 1732 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransAdd() 1733 @*/ 1734 int MatSolveTrans(Mat mat,Vec b,Vec x) 1735 { 1736 int ierr; 1737 1738 PetscFunctionBegin; 1739 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1740 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1741 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1742 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1743 if (!mat->ops->solvetrans) SETERRQ(PETSC_ERR_SUP,0,""); 1744 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 1745 if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->N,b->N); 1746 1747 PLogEventBegin(MAT_SolveTrans,mat,b,x,0); 1748 ierr = (*mat->ops->solvetrans)(mat,b,x); CHKERRQ(ierr); 1749 PLogEventEnd(MAT_SolveTrans,mat,b,x,0); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNC__ 1754 #define __FUNC__ "MatSolveTransAdd" 1755 /*@ 1756 MatSolveTransAdd - Computes x = y + inv(trans(A)) b, given a 1757 factored matrix. 1758 1759 Collective on Mat and Vec 1760 1761 Input Parameters: 1762 + mat - the factored matrix 1763 . b - the right-hand-side vector 1764 - y - the vector to be added to 1765 1766 Output Parameter: 1767 . x - the result vector 1768 1769 Notes: 1770 The vectors b and x cannot be the same. I.e., one cannot 1771 call MatSolveTransAdd(A,x,y,x). 1772 1773 Most users should employ the simplified SLES interface for linear solvers 1774 instead of working directly with matrix algebra routines such as this. 1775 See, e.g., SLESCreate(). 1776 1777 Level: developer 1778 1779 .keywords: matrix, linear system, solve, LU, Cholesky, transpose, add 1780 1781 .seealso: MatSolve(), MatSolveAdd(), MatSolveTrans() 1782 @*/ 1783 int MatSolveTransAdd(Mat mat,Vec b,Vec y,Vec x) 1784 { 1785 Scalar one = 1.0; 1786 int ierr; 1787 Vec tmp; 1788 1789 PetscFunctionBegin; 1790 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1791 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1792 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1793 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1794 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 1795 if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->N,b->N); 1796 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim %d %d",mat->N,y->N); 1797 if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim %d %d",x->n,y->n); 1798 1799 PLogEventBegin(MAT_SolveTransAdd,mat,b,x,y); 1800 if (mat->ops->solvetransadd) { 1801 ierr = (*mat->ops->solvetransadd)(mat,b,y,x); CHKERRQ(ierr); 1802 } else { 1803 /* do the solve then the add manually */ 1804 if (x != y) { 1805 ierr = MatSolveTrans(mat,b,x); CHKERRQ(ierr); 1806 ierr = VecAXPY(&one,y,x); CHKERRQ(ierr); 1807 } else { 1808 ierr = VecDuplicate(x,&tmp); CHKERRQ(ierr); 1809 PLogObjectParent(mat,tmp); 1810 ierr = VecCopy(x,tmp); CHKERRQ(ierr); 1811 ierr = MatSolveTrans(mat,b,x); CHKERRQ(ierr); 1812 ierr = VecAXPY(&one,tmp,x); CHKERRQ(ierr); 1813 ierr = VecDestroy(tmp); CHKERRQ(ierr); 1814 } 1815 } 1816 PLogEventEnd(MAT_SolveTransAdd,mat,b,x,y); 1817 PetscFunctionReturn(0); 1818 } 1819 /* ----------------------------------------------------------------*/ 1820 1821 #undef __FUNC__ 1822 #define __FUNC__ "MatRelax" 1823 /*@ 1824 MatRelax - Computes one relaxation sweep. 1825 1826 Collective on Mat and Vec 1827 1828 Input Parameters: 1829 + mat - the matrix 1830 . b - the right hand side 1831 . omega - the relaxation factor 1832 . flag - flag indicating the type of SOR (see below) 1833 . shift - diagonal shift 1834 - its - the number of iterations 1835 1836 Output Parameters: 1837 . x - the solution (can contain an initial guess) 1838 1839 SOR Flags: 1840 . SOR_FORWARD_SWEEP - forward SOR 1841 . SOR_BACKWARD_SWEEP - backward SOR 1842 . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR) 1843 . SOR_LOCAL_FORWARD_SWEEP - local forward SOR 1844 . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR 1845 . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR 1846 . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies 1847 upper/lower triangular part of matrix to 1848 vector (with omega) 1849 . SOR_ZERO_INITIAL_GUESS - zero initial guess 1850 1851 Notes: 1852 SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and 1853 SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings 1854 on each processor. 1855 1856 Application programmers will not generally use MatRelax() directly, 1857 but instead will employ the SLES/PC interface. 1858 1859 Notes for Advanced Users: 1860 The flags are implemented as bitwise inclusive or operations. 1861 For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP) 1862 to specify a zero initial guess for SSOR. 1863 1864 Most users should employ the simplified SLES interface for linear solvers 1865 instead of working directly with matrix algebra routines such as this. 1866 See, e.g., SLESCreate(). 1867 1868 Level: developer 1869 1870 .keywords: matrix, relax, relaxation, sweep 1871 @*/ 1872 int MatRelax(Mat mat,Vec b,double omega,MatSORType flag,double shift,int its,Vec x) 1873 { 1874 int ierr; 1875 1876 PetscFunctionBegin; 1877 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1878 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1879 if (!mat->ops->relax) SETERRQ(PETSC_ERR_SUP,0,""); 1880 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1881 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1882 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1883 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1884 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1885 1886 PLogEventBegin(MAT_Relax,mat,b,x,0); 1887 ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,x); CHKERRQ(ierr); 1888 PLogEventEnd(MAT_Relax,mat,b,x,0); 1889 PetscFunctionReturn(0); 1890 } 1891 1892 #undef __FUNC__ 1893 #define __FUNC__ "MatCopy_Basic" 1894 /* 1895 Default matrix copy routine. 1896 */ 1897 int MatCopy_Basic(Mat A,Mat B,MatStructure str) 1898 { 1899 int ierr,i,rstart,rend,nz,*cwork; 1900 Scalar *vwork; 1901 1902 PetscFunctionBegin; 1903 ierr = MatZeroEntries(B); CHKERRQ(ierr); 1904 ierr = MatGetOwnershipRange(A,&rstart,&rend); CHKERRQ(ierr); 1905 for (i=rstart; i<rend; i++) { 1906 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1907 ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 1908 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1909 } 1910 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1911 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1912 PetscFunctionReturn(0); 1913 } 1914 1915 #undef __FUNC__ 1916 #define __FUNC__ "MatCopy" 1917 /*@C 1918 MatCopy - Copys a matrix to another matrix. 1919 1920 Collective on Mat 1921 1922 Input Parameters: 1923 + A - the matrix 1924 - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN 1925 1926 Output Parameter: 1927 . B - where the copy is put 1928 1929 Notes: 1930 If you use SAME_NONZERO_PATTERN then the zero matrices had better have the 1931 same nonzero pattern or the routine will crash. 1932 1933 MatCopy() copies the matrix entries of a matrix to another existing 1934 matrix (after first zeroing the second matrix). A related routine is 1935 MatConvert(), which first creates a new matrix and then copies the data. 1936 1937 Level: intermediate 1938 1939 .keywords: matrix, copy, convert 1940 1941 .seealso: MatConvert() 1942 @*/ 1943 int MatCopy(Mat A,Mat B,MatStructure str) 1944 { 1945 int ierr; 1946 1947 PetscFunctionBegin; 1948 PetscValidHeaderSpecific(A,MAT_COOKIE); PetscValidHeaderSpecific(B,MAT_COOKIE); 1949 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1950 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1951 if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim %d %d",A->M,B->M, 1952 A->N,B->N); 1953 1954 PLogEventBegin(MAT_Copy,A,B,0,0); 1955 if (A->ops->copy) { 1956 ierr = (*A->ops->copy)(A,B,str); CHKERRQ(ierr); 1957 } else { /* generic conversion */ 1958 ierr = MatCopy_Basic(A,B,str); CHKERRQ(ierr); 1959 } 1960 PLogEventEnd(MAT_Copy,A,B,0,0); 1961 PetscFunctionReturn(0); 1962 } 1963 1964 static int MatConvertersSet = 0; 1965 static int (*MatConverters[MAX_MATRIX_TYPES][MAX_MATRIX_TYPES])(Mat,MatType,Mat*) = 1966 {{0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1967 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1968 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1969 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1970 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1971 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}}; 1972 1973 #undef __FUNC__ 1974 #define __FUNC__ "MatConvertRegister" 1975 /*@C 1976 MatConvertRegister - Allows one to register a routine that converts between 1977 two matrix types. 1978 1979 Not Collective 1980 1981 Input Parameters: 1982 + intype - the type of matrix (defined in include/mat.h), for example, MATSEQAIJ. 1983 - outtype - new matrix type, or MATSAME 1984 1985 Level: advanced 1986 1987 .seealso: MatConvertRegisterAll() 1988 @*/ 1989 int MatConvertRegister(MatType intype,MatType outtype,int (*converter)(Mat,MatType,Mat*)) 1990 { 1991 PetscFunctionBegin; 1992 MatConverters[intype][outtype] = converter; 1993 MatConvertersSet = 1; 1994 PetscFunctionReturn(0); 1995 } 1996 1997 #undef __FUNC__ 1998 #define __FUNC__ "MatConvert" 1999 /*@C 2000 MatConvert - Converts a matrix to another matrix, either of the same 2001 or different type. 2002 2003 Collective on Mat 2004 2005 Input Parameters: 2006 + mat - the matrix 2007 - newtype - new matrix type. Use MATSAME to create a new matrix of the 2008 same type as the original matrix. 2009 2010 Output Parameter: 2011 . M - pointer to place new matrix 2012 2013 Notes: 2014 MatConvert() first creates a new matrix and then copies the data from 2015 the first matrix. A related routine is MatCopy(), which copies the matrix 2016 entries of one matrix to another already existing matrix context. 2017 2018 Level: intermediate 2019 2020 .keywords: matrix, copy, convert 2021 2022 .seealso: MatCopy(), MatDuplicate() 2023 @*/ 2024 int MatConvert(Mat mat,MatType newtype,Mat *M) 2025 { 2026 int ierr; 2027 2028 PetscFunctionBegin; 2029 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2030 PetscValidPointer(M); 2031 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2032 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2033 2034 if (newtype > MAX_MATRIX_TYPES || newtype < -1) { 2035 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Not a valid matrix type"); 2036 } 2037 *M = 0; 2038 2039 if (!MatConvertersSet) { 2040 ierr = MatLoadRegisterAll(); CHKERRQ(ierr); 2041 } 2042 2043 PLogEventBegin(MAT_Convert,mat,0,0,0); 2044 if ((newtype == mat->type || newtype == MATSAME) && mat->ops->duplicate) { 2045 ierr = (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M); CHKERRQ(ierr); 2046 } else { 2047 if (!MatConvertersSet) { 2048 ierr = MatConvertRegisterAll(); CHKERRQ(ierr); 2049 } 2050 if (!MatConverters[mat->type][newtype]) { 2051 SETERRQ(PETSC_ERR_ARG_WRONG,1,"Invalid matrix type, or matrix converter not registered"); 2052 } 2053 ierr = (*MatConverters[mat->type][newtype])(mat,newtype,M); CHKERRQ(ierr); 2054 } 2055 PLogEventEnd(MAT_Convert,mat,0,0,0); 2056 PetscFunctionReturn(0); 2057 } 2058 2059 #undef __FUNC__ 2060 #define __FUNC__ "MatDuplicate" 2061 /*@C 2062 MatDuplicate - Duplicates a matrix including the non-zero structure. 2063 2064 Collective on Mat 2065 2066 Input Parameters: 2067 + mat - the matrix 2068 - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero 2069 values as well or not 2070 2071 Output Parameter: 2072 . M - pointer to place new matrix 2073 2074 Level: intermediate 2075 2076 .keywords: matrix, copy, convert, duplicate 2077 2078 .seealso: MatCopy(), MatConvert() 2079 @*/ 2080 int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M) 2081 { 2082 int ierr; 2083 2084 PetscFunctionBegin; 2085 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2086 PetscValidPointer(M); 2087 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2088 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2089 2090 *M = 0; 2091 PLogEventBegin(MAT_Convert,mat,0,0,0); 2092 if (!mat->ops->duplicate) { 2093 SETERRQ(PETSC_ERR_SUP,1,"Not written for this matrix type"); 2094 } 2095 ierr = (*mat->ops->duplicate)(mat,op,M); CHKERRQ(ierr); 2096 PLogEventEnd(MAT_Convert,mat,0,0,0); 2097 PetscFunctionReturn(0); 2098 } 2099 2100 #undef __FUNC__ 2101 #define __FUNC__ "MatGetDiagonal" 2102 /*@ 2103 MatGetDiagonal - Gets the diagonal of a matrix. 2104 2105 Collective on Mat and Vec 2106 2107 Input Parameters: 2108 + mat - the matrix 2109 - v - the vector for storing the diagonal 2110 2111 Output Parameter: 2112 . v - the diagonal of the matrix 2113 2114 Notes: 2115 For the SeqAIJ matrix format, this routine may also be called 2116 on a LU factored matrix; in that case it routines the reciprocal of 2117 the diagonal entries in U. It returns the entries permuted by the 2118 row and column permutation used during the symbolic factorization. 2119 2120 Level: intermediate 2121 2122 .keywords: matrix, get, diagonal 2123 2124 .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix() 2125 @*/ 2126 int MatGetDiagonal(Mat mat,Vec v) 2127 { 2128 int ierr; 2129 2130 PetscFunctionBegin; 2131 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v,VEC_COOKIE); 2132 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2133 if (!mat->ops->getdiagonal) SETERRQ(PETSC_ERR_SUP,0,""); 2134 ierr = (*mat->ops->getdiagonal)(mat,v);CHKERRQ(ierr); 2135 PetscFunctionReturn(0); 2136 } 2137 2138 #undef __FUNC__ 2139 #define __FUNC__ "MatTranspose" 2140 /*@C 2141 MatTranspose - Computes an in-place or out-of-place transpose of a matrix. 2142 2143 Collective on Mat 2144 2145 Input Parameter: 2146 . mat - the matrix to transpose 2147 2148 Output Parameters: 2149 . B - the transpose (or pass in PETSC_NULL for an in-place transpose) 2150 2151 Level: intermediate 2152 2153 .keywords: matrix, transpose 2154 2155 .seealso: MatMultTrans(), MatMultTransAdd() 2156 @*/ 2157 int MatTranspose(Mat mat,Mat *B) 2158 { 2159 int ierr; 2160 2161 PetscFunctionBegin; 2162 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2163 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2164 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2165 if (!mat->ops->transpose) SETERRQ(PETSC_ERR_SUP,0,""); 2166 ierr = (*mat->ops->transpose)(mat,B);CHKERRQ(ierr); 2167 PetscFunctionReturn(0); 2168 } 2169 2170 #undef __FUNC__ 2171 #define __FUNC__ "MatPermute" 2172 /*@C 2173 MatPermute - Creates a new matrix with rows and columns permuted from the 2174 original. 2175 2176 Collective on Mat 2177 2178 Input Parameters: 2179 + mat - the matrix to permute 2180 . row - row permutation 2181 - col - column permutation 2182 2183 Output Parameters: 2184 . B - the permuted matrix 2185 2186 Level: advanced 2187 2188 .keywords: matrix, transpose 2189 2190 .seealso: MatGetOrdering() 2191 @*/ 2192 int MatPermute(Mat mat,IS row,IS col,Mat *B) 2193 { 2194 int ierr; 2195 2196 PetscFunctionBegin; 2197 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2198 PetscValidHeaderSpecific(row,IS_COOKIE); 2199 PetscValidHeaderSpecific(col,IS_COOKIE); 2200 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2201 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2202 if (!mat->ops->permute) SETERRQ(PETSC_ERR_SUP,0,""); 2203 ierr = (*mat->ops->permute)(mat,row,col,B);CHKERRQ(ierr); 2204 PetscFunctionReturn(0); 2205 } 2206 2207 #undef __FUNC__ 2208 #define __FUNC__ "MatEqual" 2209 /*@ 2210 MatEqual - Compares two matrices. 2211 2212 Collective on Mat 2213 2214 Input Parameters: 2215 + A - the first matrix 2216 - B - the second matrix 2217 2218 Output Parameter: 2219 . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise. 2220 2221 Level: intermediate 2222 2223 .keywords: matrix, equal, equivalent 2224 @*/ 2225 int MatEqual(Mat A,Mat B,PetscTruth *flg) 2226 { 2227 int ierr; 2228 2229 PetscFunctionBegin; 2230 PetscValidHeaderSpecific(A,MAT_COOKIE); PetscValidHeaderSpecific(B,MAT_COOKIE); 2231 PetscValidIntPointer(flg); 2232 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2233 if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2234 if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim %d %d %d %d", 2235 A->M,B->M,A->N,B->N); 2236 if (!A->ops->equal) SETERRQ(PETSC_ERR_SUP,0,""); 2237 ierr = (*A->ops->equal)(A,B,flg);CHKERRQ(ierr); 2238 PetscFunctionReturn(0); 2239 } 2240 2241 #undef __FUNC__ 2242 #define __FUNC__ "MatDiagonalScale" 2243 /*@ 2244 MatDiagonalScale - Scales a matrix on the left and right by diagonal 2245 matrices that are stored as vectors. Either of the two scaling 2246 matrices can be PETSC_NULL. 2247 2248 Collective on Mat 2249 2250 Input Parameters: 2251 + mat - the matrix to be scaled 2252 . l - the left scaling vector (or PETSC_NULL) 2253 - r - the right scaling vector (or PETSC_NULL) 2254 2255 Notes: 2256 MatDiagonalScale() computes A = LAR, where 2257 L = a diagonal matrix, R = a diagonal matrix 2258 2259 Level: intermediate 2260 2261 .keywords: matrix, diagonal, scale 2262 2263 .seealso: MatScale() 2264 @*/ 2265 int MatDiagonalScale(Mat mat,Vec l,Vec r) 2266 { 2267 int ierr; 2268 2269 PetscFunctionBegin; 2270 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2271 if (!mat->ops->diagonalscale) SETERRQ(PETSC_ERR_SUP,0,""); 2272 if (l) PetscValidHeaderSpecific(l,VEC_COOKIE); 2273 if (r) PetscValidHeaderSpecific(r,VEC_COOKIE); 2274 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2275 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2276 2277 PLogEventBegin(MAT_Scale,mat,0,0,0); 2278 ierr = (*mat->ops->diagonalscale)(mat,l,r); CHKERRQ(ierr); 2279 PLogEventEnd(MAT_Scale,mat,0,0,0); 2280 PetscFunctionReturn(0); 2281 } 2282 2283 #undef __FUNC__ 2284 #define __FUNC__ "MatScale" 2285 /*@ 2286 MatScale - Scales all elements of a matrix by a given number. 2287 2288 Collective on Mat 2289 2290 Input Parameters: 2291 + mat - the matrix to be scaled 2292 - a - the scaling value 2293 2294 Output Parameter: 2295 . mat - the scaled matrix 2296 2297 Level: intermediate 2298 2299 .keywords: matrix, scale 2300 2301 .seealso: MatDiagonalScale() 2302 @*/ 2303 int MatScale(Scalar *a,Mat mat) 2304 { 2305 int ierr; 2306 2307 PetscFunctionBegin; 2308 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2309 PetscValidScalarPointer(a); 2310 if (!mat->ops->scale) SETERRQ(PETSC_ERR_SUP,0,""); 2311 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2312 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2313 2314 PLogEventBegin(MAT_Scale,mat,0,0,0); 2315 ierr = (*mat->ops->scale)(a,mat); CHKERRQ(ierr); 2316 PLogEventEnd(MAT_Scale,mat,0,0,0); 2317 PetscFunctionReturn(0); 2318 } 2319 2320 #undef __FUNC__ 2321 #define __FUNC__ "MatNorm" 2322 /*@ 2323 MatNorm - Calculates various norms of a matrix. 2324 2325 Collective on Mat 2326 2327 Input Parameters: 2328 + mat - the matrix 2329 - type - the type of norm, NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY 2330 2331 Output Parameters: 2332 . norm - the resulting norm 2333 2334 Level: intermediate 2335 2336 .keywords: matrix, norm, Frobenius 2337 @*/ 2338 int MatNorm(Mat mat,NormType type,double *norm) 2339 { 2340 int ierr; 2341 2342 PetscFunctionBegin; 2343 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2344 PetscValidScalarPointer(norm); 2345 2346 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2347 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2348 if (!mat->ops->norm) SETERRQ(PETSC_ERR_SUP,0,"Not for this matrix type"); 2349 ierr = (*mat->ops->norm)(mat,type,norm);CHKERRQ(ierr); 2350 PetscFunctionReturn(0); 2351 } 2352 2353 /* 2354 This variable is used to prevent counting of MatAssemblyBegin() that 2355 are called from within a MatAssemblyEnd(). 2356 */ 2357 static int MatAssemblyEnd_InUse = 0; 2358 #undef __FUNC__ 2359 #define __FUNC__ "MatAssemblyBegin" 2360 /*@ 2361 MatAssemblyBegin - Begins assembling the matrix. This routine should 2362 be called after completing all calls to MatSetValues(). 2363 2364 Collective on Mat 2365 2366 Input Parameters: 2367 + mat - the matrix 2368 - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 2369 2370 Notes: 2371 MatSetValues() generally caches the values. The matrix is ready to 2372 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 2373 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 2374 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 2375 using the matrix. 2376 2377 Level: beginner 2378 2379 .keywords: matrix, assembly, assemble, begin 2380 2381 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled() 2382 @*/ 2383 int MatAssemblyBegin(Mat mat,MatAssemblyType type) 2384 { 2385 int ierr; 2386 2387 PetscFunctionBegin; 2388 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2389 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix.\n did you forget to call MatSetUnfactored()?"); 2390 if (mat->assembled) { 2391 mat->was_assembled = PETSC_TRUE; 2392 mat->assembled = PETSC_FALSE; 2393 } 2394 if (!MatAssemblyEnd_InUse) { 2395 PLogEventBegin(MAT_AssemblyBegin,mat,0,0,0); 2396 if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);} 2397 PLogEventEnd(MAT_AssemblyBegin,mat,0,0,0); 2398 } else { 2399 if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);} 2400 } 2401 PetscFunctionReturn(0); 2402 } 2403 2404 #undef __FUNC__ 2405 #define __FUNC__ "MatAssembed" 2406 /*@ 2407 MatAssembled - Indicates if a matrix has been assembled and is ready for 2408 use; for example, in matrix-vector product. 2409 2410 Collective on Mat 2411 2412 Input Parameter: 2413 . mat - the matrix 2414 2415 Output Parameter: 2416 . assembled - PETSC_TRUE or PETSC_FALSE 2417 2418 Level: advanced 2419 2420 .keywords: matrix, assembly, assemble, begin 2421 2422 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin() 2423 @*/ 2424 int MatAssembled(Mat mat,PetscTruth *assembled) 2425 { 2426 int ierr; 2427 2428 PetscFunctionBegin; 2429 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2430 *assembled = mat->assembled; 2431 PetscFunctionReturn(0); 2432 } 2433 2434 #undef __FUNC__ 2435 #define __FUNC__ "MatView_Private" 2436 /* 2437 Processes command line options to determine if/how a matrix 2438 is to be viewed. Called by MatAssemblyEnd() and MatLoad(). 2439 */ 2440 int MatView_Private(Mat mat) 2441 { 2442 int ierr,flg; 2443 2444 PetscFunctionBegin; 2445 ierr = OptionsHasName(PETSC_NULL,"-mat_view_info",&flg); CHKERRQ(ierr); 2446 if (flg) { 2447 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO,0);CHKERRQ(ierr); 2448 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2449 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2450 } 2451 ierr = OptionsHasName(PETSC_NULL,"-mat_view_info_detailed",&flg);CHKERRQ(ierr); 2452 if (flg) { 2453 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO_LONG,0);CHKERRQ(ierr); 2454 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2455 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2456 } 2457 ierr = OptionsHasName(PETSC_NULL,"-mat_view",&flg); CHKERRQ(ierr); 2458 if (flg) { 2459 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2460 } 2461 ierr = OptionsHasName(PETSC_NULL,"-mat_view_matlab",&flg); CHKERRQ(ierr); 2462 if (flg) { 2463 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_MATLAB,"M");CHKERRQ(ierr); 2464 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2465 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2466 } 2467 ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg); CHKERRQ(ierr); 2468 if (flg) { 2469 ierr = OptionsHasName(0,"-mat_view_contour",&flg); CHKERRQ(ierr); 2470 if (flg) { 2471 ViewerPushFormat(VIEWER_DRAW_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr); 2472 } 2473 ierr = MatView(mat,VIEWER_DRAW_(mat->comm)); CHKERRQ(ierr); 2474 ierr = ViewerFlush(VIEWER_DRAW_(mat->comm)); CHKERRQ(ierr); 2475 if (flg) { 2476 ViewerPopFormat(VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 2477 } 2478 } 2479 ierr = OptionsHasName(PETSC_NULL,"-mat_view_socket",&flg); CHKERRQ(ierr); 2480 if (flg) { 2481 ierr = MatView(mat,VIEWER_SOCKET_(mat->comm)); CHKERRQ(ierr); 2482 ierr = ViewerFlush(VIEWER_SOCKET_(mat->comm)); CHKERRQ(ierr); 2483 } 2484 PetscFunctionReturn(0); 2485 } 2486 2487 #undef __FUNC__ 2488 #define __FUNC__ "MatAssemblyEnd" 2489 /*@ 2490 MatAssemblyEnd - Completes assembling the matrix. This routine should 2491 be called after MatAssemblyBegin(). 2492 2493 Collective on Mat 2494 2495 Input Parameters: 2496 + mat - the matrix 2497 - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 2498 2499 Options Database Keys: 2500 + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly() 2501 . -mat_view_info_detailed - Prints more detailed info 2502 . -mat_view - Prints matrix in ASCII format 2503 . -mat_view_matlab - Prints matrix in Matlab format 2504 . -mat_view_draw - Draws nonzero structure of matrix, using MatView() and DrawOpenX(). 2505 . -display <name> - Sets display name (default is host) 2506 - -draw_pause <sec> - Sets number of seconds to pause after display 2507 2508 Notes: 2509 MatSetValues() generally caches the values. The matrix is ready to 2510 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 2511 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 2512 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 2513 using the matrix. 2514 2515 Level: beginner 2516 2517 .keywords: matrix, assembly, assemble, end 2518 2519 .seealso: MatAssemblyBegin(), MatSetValues(), DrawOpenX(), MatView(), MatAssembled() 2520 @*/ 2521 int MatAssemblyEnd(Mat mat,MatAssemblyType type) 2522 { 2523 int ierr; 2524 static int inassm = 0; 2525 2526 PetscFunctionBegin; 2527 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2528 2529 inassm++; 2530 MatAssemblyEnd_InUse++; 2531 if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */ 2532 PLogEventBegin(MAT_AssemblyEnd,mat,0,0,0); 2533 if (mat->ops->assemblyend) { 2534 ierr = (*mat->ops->assemblyend)(mat,type); CHKERRQ(ierr); 2535 } 2536 PLogEventEnd(MAT_AssemblyEnd,mat,0,0,0); 2537 } else { 2538 if (mat->ops->assemblyend) { 2539 ierr = (*mat->ops->assemblyend)(mat,type); CHKERRQ(ierr); 2540 } 2541 } 2542 2543 /* Flush assembly is not a true assembly */ 2544 if (type != MAT_FLUSH_ASSEMBLY) { 2545 mat->assembled = PETSC_TRUE; mat->num_ass++; 2546 } 2547 mat->insertmode = NOT_SET_VALUES; 2548 MatAssemblyEnd_InUse--; 2549 2550 if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) { 2551 ierr = MatView_Private(mat); CHKERRQ(ierr); 2552 } 2553 inassm--; 2554 PetscFunctionReturn(0); 2555 } 2556 2557 2558 #undef __FUNC__ 2559 #define __FUNC__ "MatCompress" 2560 /*@ 2561 MatCompress - Tries to store the matrix in as little space as 2562 possible. May fail if memory is already fully used, since it 2563 tries to allocate new space. 2564 2565 Collective on Mat 2566 2567 Input Parameters: 2568 . mat - the matrix 2569 2570 Level: advanced 2571 2572 .keywords: matrix, compress 2573 @*/ 2574 int MatCompress(Mat mat) 2575 { 2576 int ierr; 2577 2578 PetscFunctionBegin; 2579 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2580 if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);} 2581 PetscFunctionReturn(0); 2582 } 2583 2584 #undef __FUNC__ 2585 #define __FUNC__ "MatSetOption" 2586 /*@ 2587 MatSetOption - Sets a parameter option for a matrix. Some options 2588 may be specific to certain storage formats. Some options 2589 determine how values will be inserted (or added). Sorted, 2590 row-oriented input will generally assemble the fastest. The default 2591 is row-oriented, nonsorted input. 2592 2593 Collective on Mat 2594 2595 Input Parameters: 2596 + mat - the matrix 2597 - option - the option, one of those listed below (and possibly others), 2598 e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR 2599 2600 Options Describing Matrix Structure: 2601 + MAT_SYMMETRIC - symmetric in terms of both structure and value 2602 - MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure 2603 2604 Options For Use with MatSetValues(): 2605 Insert a logically dense subblock, which can be 2606 + MAT_ROW_ORIENTED - row-oriented 2607 . MAT_COLUMN_ORIENTED - column-oriented 2608 . MAT_ROWS_SORTED - sorted by row 2609 . MAT_ROWS_UNSORTED - not sorted by row 2610 . MAT_COLUMNS_SORTED - sorted by column 2611 - MAT_COLUMNS_UNSORTED - not sorted by column 2612 2613 Not these options reflect the data you pass in with MatSetValues(); it has 2614 nothing to do with how the data is stored internally in the matrix 2615 data structure. 2616 2617 When (re)assembling a matrix, we can restrict the input for 2618 efficiency/debugging purposes. These options include 2619 + MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be 2620 allowed if they generate a new nonzero 2621 . MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed 2622 . MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if 2623 they generate a nonzero in a new diagonal (for block diagonal format only) 2624 . MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only) 2625 . MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries 2626 . MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry 2627 - MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly 2628 2629 Notes: 2630 Some options are relevant only for particular matrix types and 2631 are thus ignored by others. Other options are not supported by 2632 certain matrix types and will generate an error message if set. 2633 2634 If using a Fortran 77 module to compute a matrix, one may need to 2635 use the column-oriented option (or convert to the row-oriented 2636 format). 2637 2638 MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion 2639 that would generate a new entry in the nonzero structure is instead 2640 ignored. Thus, if memory has not alredy been allocated for this particular 2641 data, then the insertion is ignored. For dense matrices, in which 2642 the entire array is allocated, no entries are ever ignored. 2643 2644 MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion 2645 that would generate a new entry in the nonzero structure instead produces 2646 an error. (Currently supported for AIJ and BAIJ formats only.) 2647 This is a useful flag when using SAME_NONZERO_PATTERN in calling 2648 SLESSetOperators() to ensure that the nonzero pattern truely does 2649 remain unchanged. 2650 2651 MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion 2652 that would generate a new entry that has not been preallocated will 2653 instead produce an error. (Currently supported for AIJ and BAIJ formats 2654 only.) This is a useful flag when debugging matrix memory preallocation. 2655 2656 MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for 2657 other processors should be dropped, rather than stashed. 2658 This is useful if you know that the "owning" processor is also 2659 always generating the correct matrix entries, so that PETSc need 2660 not transfer duplicate entries generated on another processor. 2661 2662 MAT_USE_HASH_TABLE indicates that a hash table be used to improve the 2663 searches during matrix assembly. When this flag is set, the hash table 2664 is created during the first Matrix Assembly. This hash table is 2665 used the next time through, during MatSetVaules()/MatSetVaulesBlocked() 2666 to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag 2667 should be used with MAT_USE_HASH_TABLE flag. This option is currently 2668 supported by MATMPIBAIJ format only. 2669 2670 Level: intermediate 2671 2672 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero 2673 @*/ 2674 int MatSetOption(Mat mat,MatOption op) 2675 { 2676 int ierr; 2677 2678 PetscFunctionBegin; 2679 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2680 if (mat->ops->setoption) {ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr);} 2681 PetscFunctionReturn(0); 2682 } 2683 2684 #undef __FUNC__ 2685 #define __FUNC__ "MatZeroEntries" 2686 /*@ 2687 MatZeroEntries - Zeros all entries of a matrix. For sparse matrices 2688 this routine retains the old nonzero structure. 2689 2690 Collective on Mat 2691 2692 Input Parameters: 2693 . mat - the matrix 2694 2695 Level: intermediate 2696 2697 .keywords: matrix, zero, entries 2698 2699 .seealso: MatZeroRows() 2700 @*/ 2701 int MatZeroEntries(Mat mat) 2702 { 2703 int ierr; 2704 2705 PetscFunctionBegin; 2706 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2707 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2708 if (!mat->ops->zeroentries) SETERRQ(PETSC_ERR_SUP,0,""); 2709 2710 PLogEventBegin(MAT_ZeroEntries,mat,0,0,0); 2711 ierr = (*mat->ops->zeroentries)(mat); CHKERRQ(ierr); 2712 PLogEventEnd(MAT_ZeroEntries,mat,0,0,0); 2713 PetscFunctionReturn(0); 2714 } 2715 2716 #undef __FUNC__ 2717 #define __FUNC__ "MatZeroRows" 2718 /*@ 2719 MatZeroRows - Zeros all entries (except possibly the main diagonal) 2720 of a set of rows of a matrix. 2721 2722 Collective on Mat 2723 2724 Input Parameters: 2725 + mat - the matrix 2726 . is - index set of rows to remove 2727 - diag - pointer to value put in all diagonals of eliminated rows. 2728 Note that diag is not a pointer to an array, but merely a 2729 pointer to a single value. 2730 2731 Notes: 2732 For the AIJ matrix formats this removes the old nonzero structure, 2733 but does not release memory. For the dense and block diagonal 2734 formats this does not alter the nonzero structure. 2735 2736 The user can set a value in the diagonal entry (or for the AIJ and 2737 row formats can optionally remove the main diagonal entry from the 2738 nonzero structure as well, by passing a null pointer as the final 2739 argument). 2740 2741 For the parallel case, all processes that share the matrix (i.e., 2742 those in the communicator used for matrix creation) MUST call this 2743 routine, regardless of whether any rows being zeroed are owned by 2744 them. 2745 2746 Level: intermediate 2747 2748 .keywords: matrix, zero, rows, boundary conditions 2749 2750 .seealso: MatZeroEntries(), 2751 @*/ 2752 int MatZeroRows(Mat mat,IS is, Scalar *diag) 2753 { 2754 int ierr; 2755 2756 PetscFunctionBegin; 2757 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2758 PetscValidHeaderSpecific(is,IS_COOKIE); 2759 if (diag) PetscValidScalarPointer(diag); 2760 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2761 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2762 if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,""); 2763 2764 ierr = (*mat->ops->zerorows)(mat,is,diag); CHKERRQ(ierr); 2765 ierr = MatView_Private(mat); CHKERRQ(ierr); 2766 PetscFunctionReturn(0); 2767 } 2768 2769 #undef __FUNC__ 2770 #define __FUNC__ "MatZeroRowsLocal" 2771 /*@ 2772 MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal) 2773 of a set of rows of a matrix; using local numbering of rows. 2774 2775 Collective on Mat 2776 2777 Input Parameters: 2778 + mat - the matrix 2779 . is - index set of rows to remove 2780 - diag - pointer to value put in all diagonals of eliminated rows. 2781 Note that diag is not a pointer to an array, but merely a 2782 pointer to a single value. 2783 2784 Notes: 2785 For the AIJ matrix formats this removes the old nonzero structure, 2786 but does not release memory. For the dense and block diagonal 2787 formats this does not alter the nonzero structure. 2788 2789 The user can set a value in the diagonal entry (or for the AIJ and 2790 row formats can optionally remove the main diagonal entry from the 2791 nonzero structure as well, by passing a null pointer as the final 2792 argument). 2793 2794 Level: intermediate 2795 2796 .keywords: matrix, zero, rows, boundary conditions 2797 2798 .seealso: MatZeroEntries(), 2799 @*/ 2800 int MatZeroRowsLocal(Mat mat,IS is, Scalar *diag) 2801 { 2802 int ierr; 2803 IS newis; 2804 2805 PetscFunctionBegin; 2806 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2807 PetscValidHeaderSpecific(is,IS_COOKIE); 2808 if (diag) PetscValidScalarPointer(diag); 2809 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2810 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2811 if (!mat->ops->zerorows) SETERRQ(PETSC_ERR_SUP,0,""); 2812 2813 ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis); CHKERRQ(ierr); 2814 ierr = (*mat->ops->zerorows)(mat,newis,diag); CHKERRQ(ierr); 2815 ierr = ISDestroy(newis); 2816 PetscFunctionReturn(0); 2817 } 2818 2819 #undef __FUNC__ 2820 #define __FUNC__ "MatGetSize" 2821 /*@ 2822 MatGetSize - Returns the numbers of rows and columns in a matrix. 2823 2824 Not Collective 2825 2826 Input Parameter: 2827 . mat - the matrix 2828 2829 Output Parameters: 2830 + m - the number of global rows 2831 - n - the number of global columns 2832 2833 Level: beginner 2834 2835 .keywords: matrix, dimension, size, rows, columns, global, get 2836 2837 .seealso: MatGetLocalSize() 2838 @*/ 2839 int MatGetSize(Mat mat,int *m,int* n) 2840 { 2841 int ierr; 2842 2843 PetscFunctionBegin; 2844 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2845 ierr = (*mat->ops->getsize)(mat,m,n);CHKERRQ(ierr); 2846 PetscFunctionReturn(0); 2847 } 2848 2849 #undef __FUNC__ 2850 #define __FUNC__ "MatGetLocalSize" 2851 /*@ 2852 MatGetLocalSize - Returns the number of rows and columns in a matrix 2853 stored locally. This information may be implementation dependent, so 2854 use with care. 2855 2856 Not Collective 2857 2858 Input Parameters: 2859 . mat - the matrix 2860 2861 Output Parameters: 2862 + m - the number of local rows 2863 - n - the number of local columns 2864 2865 Level: beginner 2866 2867 .keywords: matrix, dimension, size, local, rows, columns, get 2868 2869 .seealso: MatGetSize() 2870 @*/ 2871 int MatGetLocalSize(Mat mat,int *m,int* n) 2872 { 2873 int ierr; 2874 2875 PetscFunctionBegin; 2876 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2877 ierr = (*mat->ops->getlocalsize)(mat,m,n);CHKERRQ(ierr); 2878 PetscFunctionReturn(0); 2879 } 2880 2881 #undef __FUNC__ 2882 #define __FUNC__ "MatGetOwnershipRange" 2883 /*@ 2884 MatGetOwnershipRange - Returns the range of matrix rows owned by 2885 this processor, assuming that the matrix is laid out with the first 2886 n1 rows on the first processor, the next n2 rows on the second, etc. 2887 For certain parallel layouts this range may not be well defined. 2888 2889 Not Collective 2890 2891 Input Parameters: 2892 . mat - the matrix 2893 2894 Output Parameters: 2895 + m - the global index of the first local row 2896 - n - one more than the global index of the last local row 2897 2898 Level: beginner 2899 2900 .keywords: matrix, get, range, ownership 2901 @*/ 2902 int MatGetOwnershipRange(Mat mat,int *m,int* n) 2903 { 2904 int ierr; 2905 2906 PetscFunctionBegin; 2907 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2908 PetscValidIntPointer(m); 2909 PetscValidIntPointer(n); 2910 if (!mat->ops->getownershiprange) SETERRQ(PETSC_ERR_SUP,0,""); 2911 ierr = (*mat->ops->getownershiprange)(mat,m,n);CHKERRQ(ierr); 2912 PetscFunctionReturn(0); 2913 } 2914 2915 #undef __FUNC__ 2916 #define __FUNC__ "MatILUFactorSymbolic" 2917 /*@ 2918 MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix. 2919 Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric() 2920 to complete the factorization. 2921 2922 Collective on Mat 2923 2924 Input Parameters: 2925 + mat - the matrix 2926 . row - row permutation 2927 . column - column permutation 2928 - info - structure containing 2929 $ levels - number of levels of fill. 2930 $ expected fill - as ratio of original fill. 2931 $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices 2932 missing diagonal entries) 2933 2934 Output Parameters: 2935 . fact - new matrix that has been symbolically factored 2936 2937 Notes: 2938 See the users manual for additional information about 2939 choosing the fill factor for better efficiency. 2940 2941 Most users should employ the simplified SLES interface for linear solvers 2942 instead of working directly with matrix algebra routines such as this. 2943 See, e.g., SLESCreate(). 2944 2945 Level: developer 2946 2947 .keywords: matrix, factor, incomplete, ILU, symbolic, fill 2948 2949 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric() 2950 @*/ 2951 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact) 2952 { 2953 int ierr; 2954 2955 PetscFunctionBegin; 2956 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2957 PetscValidPointer(fact); 2958 if (info && info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative %d",info->levels); 2959 if (info && info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",info->fill); 2960 if (!mat->ops->ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Only MatCreateMPIRowbs() matrices support parallel ILU"); 2961 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2962 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2963 2964 PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0); 2965 ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact); CHKERRQ(ierr); 2966 PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0); 2967 PetscFunctionReturn(0); 2968 } 2969 2970 #undef __FUNC__ 2971 #define __FUNC__ "MatIncompleteCholeskyFactorSymbolic" 2972 /*@ 2973 MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete 2974 Cholesky factorization for a symmetric matrix. Use 2975 MatCholeskyFactorNumeric() to complete the factorization. 2976 2977 Collective on Mat 2978 2979 Input Parameters: 2980 + mat - the matrix 2981 . perm - row and column permutation 2982 . fill - levels of fill 2983 - f - expected fill as ratio of original fill 2984 2985 Output Parameter: 2986 . fact - the factored matrix 2987 2988 Notes: 2989 Currently only no-fill factorization is supported. 2990 2991 Most users should employ the simplified SLES interface for linear solvers 2992 instead of working directly with matrix algebra routines such as this. 2993 See, e.g., SLESCreate(). 2994 2995 Level: developer 2996 2997 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill 2998 2999 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor() 3000 @*/ 3001 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,double f,int fill,Mat *fact) 3002 { 3003 int ierr; 3004 3005 PetscFunctionBegin; 3006 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3007 PetscValidPointer(fact); 3008 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 3009 if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative %d",fill); 3010 if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,0,"Expected fill less than 1.0 %g",f); 3011 if (!mat->ops->incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel"); 3012 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 3013 3014 PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0); 3015 ierr = (*mat->ops->incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr); 3016 PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0); 3017 PetscFunctionReturn(0); 3018 } 3019 3020 #undef __FUNC__ 3021 #define __FUNC__ "MatGetArray" 3022 /*@C 3023 MatGetArray - Returns a pointer to the element values in the matrix. 3024 The result of this routine is dependent on the underlying matrix data 3025 structure, and may not even work for certain matrix types. You MUST 3026 call MatRestoreArray() when you no longer need to access the array. 3027 3028 Not Collective 3029 3030 Input Parameter: 3031 . mat - the matrix 3032 3033 Output Parameter: 3034 . v - the location of the values 3035 3036 Currently returns an array only for the dense formats, giving access to 3037 the local portion of the matrix in the usual Fortran column-oriented format. 3038 3039 Fortran Note: 3040 This routine is used differently from Fortran, e.g., 3041 .vb 3042 Mat mat 3043 Scalar mat_array(1) 3044 PetscOffset i_mat 3045 int ierr 3046 call MatGetArray(mat,mat_array,i_mat,ierr) 3047 3048 C Access first local entry in matrix; note that array is 3049 C treated as one dimensional 3050 value = mat_array(i_mat + 1) 3051 3052 [... other code ...] 3053 call MatRestoreArray(mat,mat_array,i_mat,ierr) 3054 .ve 3055 3056 See the Fortran chapter of the users manual and 3057 petsc/src/mat/examples/tests for details. 3058 3059 Level: advanced 3060 3061 .keywords: matrix, array, elements, values 3062 3063 .seealso: MatRestoreArray(), MatGetArrayF90() 3064 @*/ 3065 int MatGetArray(Mat mat,Scalar **v) 3066 { 3067 int ierr; 3068 3069 PetscFunctionBegin; 3070 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3071 PetscValidPointer(v); 3072 if (!mat->ops->getarray) SETERRQ(PETSC_ERR_SUP,0,""); 3073 ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr); 3074 PetscFunctionReturn(0); 3075 } 3076 3077 #undef __FUNC__ 3078 #define __FUNC__ "MatRestoreArray" 3079 /*@C 3080 MatRestoreArray - Restores the matrix after MatGetArray() has been called. 3081 3082 Not Collective 3083 3084 Input Parameter: 3085 + mat - the matrix 3086 - v - the location of the values 3087 3088 Fortran Note: 3089 This routine is used differently from Fortran, e.g., 3090 .vb 3091 Mat mat 3092 Scalar mat_array(1) 3093 PetscOffset i_mat 3094 int ierr 3095 call MatGetArray(mat,mat_array,i_mat,ierr) 3096 3097 C Access first local entry in matrix; note that array is 3098 C treated as one dimensional 3099 value = mat_array(i_mat + 1) 3100 3101 [... other code ...] 3102 call MatRestoreArray(mat,mat_array,i_mat,ierr) 3103 .ve 3104 3105 See the Fortran chapter of the users manual and 3106 petsc/src/mat/examples/tests for details 3107 3108 Level: advanced 3109 3110 .keywords: matrix, array, elements, values, restore 3111 3112 .seealso: MatGetArray(), MatRestoreArrayF90() 3113 @*/ 3114 int MatRestoreArray(Mat mat,Scalar **v) 3115 { 3116 int ierr; 3117 3118 PetscFunctionBegin; 3119 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3120 PetscValidPointer(v); 3121 if (!mat->ops->restorearray) SETERRQ(PETSC_ERR_SUP,0,""); 3122 ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr); 3123 PetscFunctionReturn(0); 3124 } 3125 3126 #undef __FUNC__ 3127 #define __FUNC__ "MatGetSubMatrices" 3128 /*@C 3129 MatGetSubMatrices - Extracts several submatrices from a matrix. If submat 3130 points to an array of valid matrices, they may be reused to store the new 3131 submatrices. 3132 3133 Collective on Mat 3134 3135 Input Parameters: 3136 + mat - the matrix 3137 . n - the number of submatrixes to be extracted 3138 . irow, icol - index sets of rows and columns to extract 3139 - scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 3140 3141 Output Parameter: 3142 . submat - the array of submatrices 3143 3144 Notes: 3145 MatGetSubMatrices() can extract only sequential submatrices 3146 (from both sequential and parallel matrices). Use MatGetSubMatrix() 3147 to extract a parallel submatrix. 3148 3149 When extracting submatrices from a parallel matrix, each processor can 3150 form a different submatrix by setting the rows and columns of its 3151 individual index sets according to the local submatrix desired. 3152 3153 When finished using the submatrices, the user should destroy 3154 them with MatDestroySubMatrices(). 3155 3156 MAT_REUSE_MATRIX can only be used when the nonzero structure of the 3157 original matrix has not changed from that last call to MatGetSubMatrices() 3158 3159 Fortran Note: 3160 The Fortran interface is slightly different from that given below, it 3161 requires one to pass in as submat a Mat (integer) array of size at least m. 3162 3163 Level: advanced 3164 3165 .keywords: matrix, get, submatrix, submatrices 3166 3167 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal() 3168 @*/ 3169 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat) 3170 { 3171 int ierr; 3172 3173 PetscFunctionBegin; 3174 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3175 if (!mat->ops->getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,""); 3176 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 3177 3178 PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0); 3179 ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat); CHKERRQ(ierr); 3180 PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0); 3181 3182 PetscFunctionReturn(0); 3183 } 3184 3185 #undef __FUNC__ 3186 #define __FUNC__ "MatDestroyMatrices" 3187 /*@C 3188 MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices(). 3189 3190 Collective on Mat 3191 3192 Input Parameters: 3193 + n - the number of local matrices 3194 - mat - the matrices 3195 3196 Level: advanced 3197 3198 .keywords: matrix, destroy, submatrix, submatrices 3199 3200 .seealso: MatGetSubMatrices() 3201 @*/ 3202 int MatDestroyMatrices(int n,Mat **mat) 3203 { 3204 int ierr,i; 3205 3206 PetscFunctionBegin; 3207 if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices %d",n); 3208 PetscValidPointer(mat); 3209 for ( i=0; i<n; i++ ) { 3210 ierr = MatDestroy((*mat)[i]); CHKERRQ(ierr); 3211 } 3212 if (n) PetscFree(*mat); 3213 PetscFunctionReturn(0); 3214 } 3215 3216 #undef __FUNC__ 3217 #define __FUNC__ "MatIncreaseOverlap" 3218 /*@ 3219 MatIncreaseOverlap - Given a set of submatrices indicated by index sets, 3220 replaces the index sets by larger ones that represent submatrices with 3221 additional overlap. 3222 3223 Collective on Mat 3224 3225 Input Parameters: 3226 + mat - the matrix 3227 . n - the number of index sets 3228 . is - the array of pointers to index sets 3229 - ov - the additional overlap requested 3230 3231 Level: developer 3232 3233 .keywords: matrix, overlap, Schwarz 3234 3235 .seealso: MatGetSubMatrices() 3236 @*/ 3237 int MatIncreaseOverlap(Mat mat,int n, IS *is,int ov) 3238 { 3239 int ierr; 3240 3241 PetscFunctionBegin; 3242 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3243 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 3244 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 3245 3246 if (ov == 0) PetscFunctionReturn(0); 3247 if (!mat->ops->increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,""); 3248 PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0); 3249 ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov); CHKERRQ(ierr); 3250 PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0); 3251 PetscFunctionReturn(0); 3252 } 3253 3254 #undef __FUNC__ 3255 #define __FUNC__ "MatPrintHelp" 3256 /*@ 3257 MatPrintHelp - Prints all the options for the matrix. 3258 3259 Collective on Mat 3260 3261 Input Parameter: 3262 . mat - the matrix 3263 3264 Options Database Keys: 3265 + -help - Prints matrix options 3266 - -h - Prints matrix options 3267 3268 Level: developer 3269 3270 .keywords: mat, help 3271 3272 .seealso: MatCreate(), MatCreateXXX() 3273 @*/ 3274 int MatPrintHelp(Mat mat) 3275 { 3276 static int called = 0; 3277 int ierr; 3278 MPI_Comm comm; 3279 3280 PetscFunctionBegin; 3281 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3282 3283 comm = mat->comm; 3284 if (!called) { 3285 ierr = (*PetscHelpPrintf)(comm,"General matrix options:\n"); CHKERRQ(ierr); 3286 ierr = (*PetscHelpPrintf)(comm," -mat_view_info: view basic matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr); 3287 ierr = (*PetscHelpPrintf)(comm," -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr); 3288 ierr = (*PetscHelpPrintf)(comm," -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n");CHKERRQ(ierr); 3289 ierr = (*PetscHelpPrintf)(comm," -draw_pause <sec>: set seconds of display pause\n");CHKERRQ(ierr); 3290 ierr = (*PetscHelpPrintf)(comm," -display <name>: set alternate display\n");CHKERRQ(ierr); 3291 called = 1; 3292 } 3293 if (mat->ops->printhelp) { 3294 ierr = (*mat->ops->printhelp)(mat); CHKERRQ(ierr); 3295 } 3296 PetscFunctionReturn(0); 3297 } 3298 3299 #undef __FUNC__ 3300 #define __FUNC__ "MatGetBlockSize" 3301 /*@ 3302 MatGetBlockSize - Returns the matrix block size; useful especially for the 3303 block row and block diagonal formats. 3304 3305 Not Collective 3306 3307 Input Parameter: 3308 . mat - the matrix 3309 3310 Output Parameter: 3311 . bs - block size 3312 3313 Notes: 3314 Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG. 3315 Block row formats are MATSEQBAIJ, MATMPIBAIJ 3316 3317 Level: intermediate 3318 3319 .keywords: matrix, get, block, size 3320 3321 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag() 3322 @*/ 3323 int MatGetBlockSize(Mat mat,int *bs) 3324 { 3325 int ierr; 3326 3327 PetscFunctionBegin; 3328 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3329 PetscValidIntPointer(bs); 3330 if (!mat->ops->getblocksize) SETERRQ(PETSC_ERR_SUP,0,""); 3331 ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr); 3332 PetscFunctionReturn(0); 3333 } 3334 3335 #undef __FUNC__ 3336 #define __FUNC__ "MatGetRowIJ" 3337 /*@C 3338 MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices. 3339 EXPERTS ONLY. 3340 3341 Collective on Mat 3342 3343 Input Parameters: 3344 + mat - the matrix 3345 . shift - 0 or 1 indicating we want the indices starting at 0 or 1 3346 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 3347 symmetrized 3348 3349 Output Parameters: 3350 + n - number of rows in the (possibly compressed) matrix 3351 . ia - the row pointers 3352 . ja - the column indices 3353 - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned 3354 3355 Level: developer 3356 3357 .seealso: MatGetColumnIJ(), MatRestoreRowIJ() 3358 @*/ 3359 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 3360 { 3361 int ierr; 3362 3363 PetscFunctionBegin; 3364 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3365 if (ia) PetscValidIntPointer(ia); 3366 if (ja) PetscValidIntPointer(ja); 3367 PetscValidIntPointer(done); 3368 if (!mat->ops->getrowij) *done = PETSC_FALSE; 3369 else { 3370 *done = PETSC_TRUE; 3371 ierr = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 3372 } 3373 PetscFunctionReturn(0); 3374 } 3375 3376 #undef __FUNC__ 3377 #define __FUNC__ "MatGetColumnIJ" 3378 /*@C 3379 MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices. 3380 EXPERTS ONLY. 3381 3382 Collective on Mat 3383 3384 Input Parameters: 3385 + mat - the matrix 3386 . shift - 1 or zero indicating we want the indices starting at 0 or 1 3387 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 3388 symmetrized 3389 3390 Output Parameters: 3391 + n - number of columns in the (possibly compressed) matrix 3392 . ia - the column pointers 3393 . ja - the row indices 3394 - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned 3395 3396 Level: developer 3397 3398 .seealso: MatGetRowIJ(), MatRestoreColumnIJ() 3399 @*/ 3400 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 3401 { 3402 int ierr; 3403 3404 PetscFunctionBegin; 3405 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3406 if (ia) PetscValidIntPointer(ia); 3407 if (ja) PetscValidIntPointer(ja); 3408 PetscValidIntPointer(done); 3409 3410 if (!mat->ops->getcolumnij) *done = PETSC_FALSE; 3411 else { 3412 *done = PETSC_TRUE; 3413 ierr = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 3414 } 3415 PetscFunctionReturn(0); 3416 } 3417 3418 #undef __FUNC__ 3419 #define __FUNC__ "MatRestoreRowIJ" 3420 /*@C 3421 MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with 3422 MatGetRowIJ(). EXPERTS ONLY. 3423 3424 Collective on Mat 3425 3426 Input Parameters: 3427 + mat - the matrix 3428 . shift - 1 or zero indicating we want the indices starting at 0 or 1 3429 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 3430 symmetrized 3431 3432 Output Parameters: 3433 + n - size of (possibly compressed) matrix 3434 . ia - the row pointers 3435 . ja - the column indices 3436 - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 3437 3438 Level: developer 3439 3440 .seealso: MatGetRowIJ(), MatRestoreColumnIJ() 3441 @*/ 3442 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 3443 { 3444 int ierr; 3445 3446 PetscFunctionBegin; 3447 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3448 if (ia) PetscValidIntPointer(ia); 3449 if (ja) PetscValidIntPointer(ja); 3450 PetscValidIntPointer(done); 3451 3452 if (!mat->ops->restorerowij) *done = PETSC_FALSE; 3453 else { 3454 *done = PETSC_TRUE; 3455 ierr = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 3456 } 3457 PetscFunctionReturn(0); 3458 } 3459 3460 #undef __FUNC__ 3461 #define __FUNC__ "MatRestoreColumnIJ" 3462 /*@C 3463 MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with 3464 MatGetColumnIJ(). EXPERTS ONLY. 3465 3466 Collective on Mat 3467 3468 Input Parameters: 3469 + mat - the matrix 3470 . shift - 1 or zero indicating we want the indices starting at 0 or 1 3471 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 3472 symmetrized 3473 3474 Output Parameters: 3475 + n - size of (possibly compressed) matrix 3476 . ia - the column pointers 3477 . ja - the row indices 3478 - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 3479 3480 Level: developer 3481 3482 .seealso: MatGetColumnIJ(), MatRestoreRowIJ() 3483 @*/ 3484 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 3485 { 3486 int ierr; 3487 3488 PetscFunctionBegin; 3489 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3490 if (ia) PetscValidIntPointer(ia); 3491 if (ja) PetscValidIntPointer(ja); 3492 PetscValidIntPointer(done); 3493 3494 if (!mat->ops->restorecolumnij) *done = PETSC_FALSE; 3495 else { 3496 *done = PETSC_TRUE; 3497 ierr = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 3498 } 3499 PetscFunctionReturn(0); 3500 } 3501 3502 #undef __FUNC__ 3503 #define __FUNC__ "MatColoringPatch" 3504 /*@C 3505 MatColoringPatch - EXPERTS ONLY, used inside matrix coloring routines that 3506 use MatGetRowIJ() and/or MatGetColumnIJ(). 3507 3508 Collective on Mat 3509 3510 Input Parameters: 3511 + mat - the matrix 3512 . n - number of colors 3513 - colorarray - array indicating color for each column 3514 3515 Output Parameters: 3516 . iscoloring - coloring generated using colorarray information 3517 3518 Level: developer 3519 3520 .seealso: MatGetRowIJ(), MatGetColumnIJ() 3521 3522 .keywords: mat, coloring, patch 3523 @*/ 3524 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring) 3525 { 3526 int ierr; 3527 3528 PetscFunctionBegin; 3529 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3530 PetscValidIntPointer(colorarray); 3531 3532 if (!mat->ops->coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");} 3533 else { 3534 ierr = (*mat->ops->coloringpatch)(mat,n,colorarray,iscoloring); CHKERRQ(ierr); 3535 } 3536 PetscFunctionReturn(0); 3537 } 3538 3539 3540 #undef __FUNC__ 3541 #define __FUNC__ "MatSetUnfactored" 3542 /*@ 3543 MatSetUnfactored - Resets a factored matrix to be treated as unfactored. 3544 3545 Collective on Mat 3546 3547 Input Parameter: 3548 . mat - the factored matrix to be reset 3549 3550 Notes: 3551 This routine should be used only with factored matrices formed by in-place 3552 factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE 3553 format). This option can save memory, for example, when solving nonlinear 3554 systems with a matrix-free Newton-Krylov method and a matrix-based, in-place 3555 ILU(0) preconditioner. 3556 3557 Note that one can specify in-place ILU(0) factorization by calling 3558 .vb 3559 PCType(pc,PCILU); 3560 PCILUSeUseInPlace(pc); 3561 .ve 3562 or by using the options -pc_type ilu -pc_ilu_in_place 3563 3564 In-place factorization ILU(0) can also be used as a local 3565 solver for the blocks within the block Jacobi or additive Schwarz 3566 methods (runtime option: -sub_pc_ilu_in_place). See the discussion 3567 of these preconditioners in the users manual for details on setting 3568 local solver options. 3569 3570 Most users should employ the simplified SLES interface for linear solvers 3571 instead of working directly with matrix algebra routines such as this. 3572 See, e.g., SLESCreate(). 3573 3574 Level: developer 3575 3576 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace() 3577 3578 .keywords: matrix-free, in-place ILU, in-place LU 3579 @*/ 3580 int MatSetUnfactored(Mat mat) 3581 { 3582 int ierr; 3583 3584 PetscFunctionBegin; 3585 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3586 mat->factor = 0; 3587 if (!mat->ops->setunfactored) PetscFunctionReturn(0); 3588 ierr = (*mat->ops->setunfactored)(mat); CHKERRQ(ierr); 3589 PetscFunctionReturn(0); 3590 } 3591 3592 #undef __FUNC__ 3593 #define __FUNC__ "MatGetType" 3594 /*@C 3595 MatGetType - Gets the matrix type and name (as a string) from the matrix. 3596 3597 Not Collective 3598 3599 Input Parameter: 3600 . mat - the matrix 3601 3602 Output Parameter: 3603 + type - the matrix type (or use PETSC_NULL) 3604 - name - name of matrix type (or use PETSC_NULL) 3605 3606 Level: intermediate 3607 3608 .keywords: matrix, get, type, name 3609 @*/ 3610 int MatGetType(Mat mat,MatType *type,char **name) 3611 { 3612 int itype = (int)mat->type; 3613 char *matname[10]; 3614 3615 PetscFunctionBegin; 3616 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3617 3618 if (type) *type = (MatType) mat->type; 3619 if (name) { 3620 /* Note: Be sure that this list corresponds to the enum in mat.h */ 3621 matname[0] = "MATSEQDENSE"; 3622 matname[1] = "MATSEQAIJ"; 3623 matname[2] = "MATMPIAIJ"; 3624 matname[3] = "MATSHELL"; 3625 matname[4] = "MATMPIROWBS"; 3626 matname[5] = "MATSEQBDIAG"; 3627 matname[6] = "MATMPIBDIAG"; 3628 matname[7] = "MATMPIDENSE"; 3629 matname[8] = "MATSEQBAIJ"; 3630 matname[9] = "MATMPIBAIJ"; 3631 3632 if (itype < 0 || itype > 9) *name = "Unknown matrix type"; 3633 else *name = matname[itype]; 3634 } 3635 PetscFunctionReturn(0); 3636 } 3637 3638 /*MC 3639 MatGetArrayF90 - Accesses a matrix array from Fortran90. 3640 3641 Synopsis: 3642 MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 3643 3644 Not collective 3645 3646 Input Parameter: 3647 . x - matrix 3648 3649 Output Parameters: 3650 + xx_v - the Fortran90 pointer to the array 3651 - ierr - error code 3652 3653 Example of Usage: 3654 .vb 3655 Scalar, pointer xx_v(:) 3656 .... 3657 call MatGetArrayF90(x,xx_v,ierr) 3658 a = xx_v(3) 3659 call MatRestoreArrayF90(x,xx_v,ierr) 3660 .ve 3661 3662 Notes: 3663 Not yet supported for all F90 compilers 3664 3665 Level: advanced 3666 3667 .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray() 3668 3669 .keywords: matrix, array, f90 3670 M*/ 3671 3672 /*MC 3673 MatRestoreArrayF90 - Restores a matrix array that has been 3674 accessed with MatGetArrayF90(). 3675 3676 Synopsis: 3677 MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 3678 3679 Not collective 3680 3681 Input Parameters: 3682 + x - matrix 3683 - xx_v - the Fortran90 pointer to the array 3684 3685 Output Parameter: 3686 . ierr - error code 3687 3688 Example of Usage: 3689 .vb 3690 Scalar, pointer xx_v(:) 3691 .... 3692 call MatGetArrayF90(x,xx_v,ierr) 3693 a = xx_v(3) 3694 call MatRestoreArrayF90(x,xx_v,ierr) 3695 .ve 3696 3697 Notes: 3698 Not yet supported for all F90 compilers 3699 3700 Level: advanced 3701 3702 .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray() 3703 3704 .keywords: matrix, array, f90 3705 M*/ 3706 3707 3708 #undef __FUNC__ 3709 #define __FUNC__ "MatGetSubMatrix" 3710 /*@ 3711 MatGetSubMatrix - Gets a single submatrix on the same number of processors 3712 as the original matrix. 3713 3714 Collective on Mat 3715 3716 Input Parameters: 3717 + mat - the original matrix 3718 . isrow - rows this processor should obtain 3719 . iscol - columns for all processors you wish to keep 3720 . csize - number of columns "local" to this processor (does nothing for sequential 3721 matrices). This should match the result from VecGetLocalSize(x,...) if you 3722 plan to use the matrix in a A*x or use PETSC_DECIDE 3723 - cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 3724 3725 Output Parameter: 3726 . newmat - the new submatrix, of the same type as the old 3727 3728 Level: advanced 3729 3730 .keywords: matrix, get, submatrix, submatrices 3731 3732 .seealso: MatGetSubMatrices() 3733 @*/ 3734 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat) 3735 { 3736 int ierr, size; 3737 Mat *local; 3738 3739 PetscFunctionBegin; 3740 ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr); 3741 3742 /* if original matrix is on just one processor then use submatrix generated */ 3743 if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) { 3744 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr); 3745 PetscFunctionReturn(0); 3746 } else if (!mat->ops->getsubmatrix && size == 1) { 3747 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr); 3748 *newmat = *local; 3749 PetscFree(local); 3750 PetscFunctionReturn(0); 3751 } 3752 3753 if (!mat->ops->getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented"); 3754 ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr); 3755 PetscFunctionReturn(0); 3756 } 3757 3758 #undef __FUNC__ 3759 #define __FUNC__ "MatGetMaps" 3760 /*@C 3761 MatGetMaps - Returns the maps associated with the matrix. 3762 3763 Not Collective 3764 3765 Input Parameter: 3766 . mat - the matrix 3767 3768 Output Parameters: 3769 + rmap - the row (right) map 3770 - cmap - the column (left) map 3771 3772 Level: developer 3773 3774 .keywords: matrix, get, map 3775 @*/ 3776 int MatGetMaps(Mat mat,Map *rmap,Map *cmap) 3777 { 3778 int ierr; 3779 3780 PetscFunctionBegin; 3781 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3782 ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr); 3783 PetscFunctionReturn(0); 3784 } 3785 3786 /* 3787 Version that works for all PETSc matrices 3788 */ 3789 #undef __FUNC__ 3790 #define __FUNC__ "MatGetMaps_Petsc" 3791 int MatGetMaps_Petsc(Mat mat,Map *rmap,Map *cmap) 3792 { 3793 PetscFunctionBegin; 3794 if (rmap) *rmap = mat->rmap; 3795 if (cmap) *cmap = mat->cmap; 3796 PetscFunctionReturn(0); 3797 } 3798 3799 #undef __FUNC__ 3800 #define __FUNC__ "MatSetStashInitialSize" 3801 /*@ 3802 MatSetStashInitialSize - sets the sizes of the matrix stash, that is 3803 used during the assembly process to store values that belong to 3804 other processors. 3805 3806 Not Collective 3807 3808 Input Parameters: 3809 + mat - the matrix 3810 . size - the initial size of the stash. 3811 - bsize - the initial size of the block-stash(if used). 3812 3813 Options Database Keys: 3814 + -matstash_initial_size <size> or <size0,size1,...sizep-1> 3815 - -matstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1> 3816 3817 Level: intermediate 3818 3819 Notes: 3820 The block-stash is used for values set with VecSetValuesBlocked() while 3821 the stash is used for values set with VecSetValues() 3822 3823 Run with the option -log_info and look for output of the form 3824 MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs. 3825 to determine the appropriate value, MM, to use for size and 3826 MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs. 3827 to determine the value, BMM to use for bsize 3828 3829 .keywords: matrix, stash, assembly 3830 @*/ 3831 int MatSetStashInitialSize(Mat mat,int size, int bsize) 3832 { 3833 int ierr; 3834 3835 PetscFunctionBegin; 3836 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3837 ierr = MatStashSetInitialSize_Private(&mat->stash,size); CHKERRQ(ierr); 3838 ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize); CHKERRQ(ierr); 3839 PetscFunctionReturn(0); 3840 } 3841