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