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