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