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