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