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