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