1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: matrix.c,v 1.277 1998/01/29 15:44:49 curfman Exp curfman $"; 3 #endif 4 5 /* 6 This is where the abstract matrix operations are defined 7 */ 8 9 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 10 #include "src/vec/vecimpl.h" 11 #include "pinclude/pviewer.h" 12 13 #undef __FUNC__ 14 #define __FUNC__ "MatGetRow" 15 /*@C 16 MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow() 17 for each row that you get to ensure that your application does 18 not bleed memory. 19 20 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->view) {ierr = (*mat->view)((PetscObject)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->destroy)((PetscObject)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 MatSetOptions() 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: MatSetOptions(), 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 MatSetOptions() 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: MatSetOptions(), 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 2200 Notes: 2201 Some options are relevant only for particular matrix types and 2202 are thus ignored by others. Other options are not supported by 2203 certain matrix types and will generate an error message if set. 2204 2205 If using a Fortran 77 module to compute a matrix, one may need to 2206 use the column-oriented option (or convert to the row-oriented 2207 format). 2208 2209 MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion 2210 that would generate a new entry in the nonzero structure is instead 2211 ignored. Thus, if memory has not alredy been allocated for this particular 2212 data, then the insertion is ignored. For dense matrices, in which 2213 the entire array is allocated, no entries are ever ignored. 2214 2215 MAT_NEW_NONZERO_LOCATION_ERROR indicates that any add or insertion 2216 that would generate a new entry in the nonzero structure instead produces 2217 an error. (Currently supported for AIJ and BAIJ formats only.) 2218 This is a useful flag when using SAME_NONZERO_PATTERN in calling 2219 SLESSetOperators() to ensure that the nonzero pattern truely does 2220 remain unchanged. 2221 2222 MAT_NEW_NONZERO_ALLOCATION_ERROR indicates that any add or insertion 2223 that would generate a new entry that has not been preallocated will 2224 instead produce an error. (Currently supported for AIJ and BAIJ formats 2225 only.) This is a useful flag when debugging matrix memory preallocation. 2226 2227 MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for 2228 other processors should be dropped, rather than stashed. 2229 This is useful if you know that the "owning" processor is also 2230 always generating the correct matrix entries, so that PETSc need 2231 not transfer duplicate entries generated on another processor. 2232 2233 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero 2234 @*/ 2235 int MatSetOption(Mat mat,MatOption op) 2236 { 2237 int ierr; 2238 2239 PetscFunctionBegin; 2240 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2241 if (mat->ops.setoption) {ierr = (*mat->ops.setoption)(mat,op);CHKERRQ(ierr);} 2242 PetscFunctionReturn(0); 2243 } 2244 2245 #undef __FUNC__ 2246 #define __FUNC__ "MatZeroEntries" 2247 /*@ 2248 MatZeroEntries - Zeros all entries of a matrix. For sparse matrices 2249 this routine retains the old nonzero structure. 2250 2251 Input Parameters: 2252 . mat - the matrix 2253 2254 .keywords: matrix, zero, entries 2255 2256 .seealso: MatZeroRows() 2257 @*/ 2258 int MatZeroEntries(Mat mat) 2259 { 2260 int ierr; 2261 2262 PetscFunctionBegin; 2263 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2264 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2265 if (!mat->ops.zeroentries) SETERRQ(PETSC_ERR_SUP,0,""); 2266 2267 PLogEventBegin(MAT_ZeroEntries,mat,0,0,0); 2268 ierr = (*mat->ops.zeroentries)(mat); CHKERRQ(ierr); 2269 PLogEventEnd(MAT_ZeroEntries,mat,0,0,0); 2270 PetscFunctionReturn(0); 2271 } 2272 2273 #undef __FUNC__ 2274 #define __FUNC__ "MatZeroRows" 2275 /*@ 2276 MatZeroRows - Zeros all entries (except possibly the main diagonal) 2277 of a set of rows of a matrix. 2278 2279 Input Parameters: 2280 . mat - the matrix 2281 . is - index set of rows to remove 2282 . diag - pointer to value put in all diagonals of eliminated rows. 2283 Note that diag is not a pointer to an array, but merely a 2284 pointer to a single value. 2285 2286 Notes: 2287 For the AIJ matrix formats this removes the old nonzero structure, 2288 but does not release memory. For the dense and block diagonal 2289 formats this does not alter the nonzero structure. 2290 2291 The user can set a value in the diagonal entry (or for the AIJ and 2292 row formats can optionally remove the main diagonal entry from the 2293 nonzero structure as well, by passing a null pointer as the final 2294 argument). 2295 2296 For the parallel case, all processes that share the matrix (i.e., 2297 those in the communicator used for matrix creation) MUST call this 2298 routine, regardless of whether any rows being zeroed are owned by 2299 them. 2300 2301 .keywords: matrix, zero, rows, boundary conditions 2302 2303 .seealso: MatZeroEntries(), 2304 @*/ 2305 int MatZeroRows(Mat mat,IS is, Scalar *diag) 2306 { 2307 int ierr; 2308 2309 PetscFunctionBegin; 2310 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2311 PetscValidHeaderSpecific(is,IS_COOKIE); 2312 if (diag) PetscValidScalarPointer(diag); 2313 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2314 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2315 if (!mat->ops.zerorows) SETERRQ(PETSC_ERR_SUP,0,""); 2316 2317 ierr = (*mat->ops.zerorows)(mat,is,diag); CHKERRQ(ierr); 2318 ierr = MatView_Private(mat); CHKERRQ(ierr); 2319 PetscFunctionReturn(0); 2320 } 2321 2322 #undef __FUNC__ 2323 #define __FUNC__ "MatZeroRowsLocal" 2324 /*@ 2325 MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal) 2326 of a set of rows of a matrix; using local numbering of rows. 2327 2328 Input Parameters: 2329 . mat - the matrix 2330 . is - index set of rows to remove 2331 . diag - pointer to value put in all diagonals of eliminated rows. 2332 Note that diag is not a pointer to an array, but merely a 2333 pointer to a single value. 2334 2335 Notes: 2336 For the AIJ matrix formats this removes the old nonzero structure, 2337 but does not release memory. For the dense and block diagonal 2338 formats this does not alter the nonzero structure. 2339 2340 The user can set a value in the diagonal entry (or for the AIJ and 2341 row formats can optionally remove the main diagonal entry from the 2342 nonzero structure as well, by passing a null pointer as the final 2343 argument). 2344 2345 .keywords: matrix, zero, rows, boundary conditions 2346 2347 .seealso: MatZeroEntries(), 2348 @*/ 2349 int MatZeroRowsLocal(Mat mat,IS is, Scalar *diag) 2350 { 2351 int ierr; 2352 IS newis; 2353 2354 PetscFunctionBegin; 2355 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2356 PetscValidHeaderSpecific(is,IS_COOKIE); 2357 if (diag) PetscValidScalarPointer(diag); 2358 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2359 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2360 if (!mat->ops.zerorows) SETERRQ(PETSC_ERR_SUP,0,""); 2361 2362 ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis); CHKERRQ(ierr); 2363 ierr = (*mat->ops.zerorows)(mat,newis,diag); CHKERRQ(ierr); 2364 ierr = ISDestroy(newis); 2365 PetscFunctionReturn(0); 2366 } 2367 2368 #undef __FUNC__ 2369 #define __FUNC__ "MatGetSize" 2370 /*@ 2371 MatGetSize - Returns the numbers of rows and columns in a matrix. 2372 2373 Input Parameter: 2374 . mat - the matrix 2375 2376 Output Parameters: 2377 . m - the number of global rows 2378 . n - the number of global columns 2379 2380 .keywords: matrix, dimension, size, rows, columns, global, get 2381 2382 .seealso: MatGetLocalSize() 2383 @*/ 2384 int MatGetSize(Mat mat,int *m,int* n) 2385 { 2386 int ierr; 2387 2388 PetscFunctionBegin; 2389 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2390 ierr = (*mat->ops.getsize)(mat,m,n);CHKERRQ(ierr); 2391 PetscFunctionReturn(0); 2392 } 2393 2394 #undef __FUNC__ 2395 #define __FUNC__ "MatGetLocalSize" 2396 /*@ 2397 MatGetLocalSize - Returns the number of rows and columns in a matrix 2398 stored locally. This information may be implementation dependent, so 2399 use with care. 2400 2401 Input Parameters: 2402 . mat - the matrix 2403 2404 Output Parameters: 2405 . m - the number of local rows 2406 . n - the number of local columns 2407 2408 .keywords: matrix, dimension, size, local, rows, columns, get 2409 2410 .seealso: MatGetSize() 2411 @*/ 2412 int MatGetLocalSize(Mat mat,int *m,int* n) 2413 { 2414 int ierr; 2415 2416 PetscFunctionBegin; 2417 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2418 ierr = (*mat->ops.getlocalsize)(mat,m,n);CHKERRQ(ierr); 2419 PetscFunctionReturn(0); 2420 } 2421 2422 #undef __FUNC__ 2423 #define __FUNC__ "MatGetOwnershipRange" 2424 /*@ 2425 MatGetOwnershipRange - Returns the range of matrix rows owned by 2426 this processor, assuming that the matrix is laid out with the first 2427 n1 rows on the first processor, the next n2 rows on the second, etc. 2428 For certain parallel layouts this range may not be well defined. 2429 2430 Input Parameters: 2431 . mat - the matrix 2432 2433 Output Parameters: 2434 . m - the global index of the first local row 2435 . n - one more than the global index of the last local row 2436 2437 .keywords: matrix, get, range, ownership 2438 @*/ 2439 int MatGetOwnershipRange(Mat mat,int *m,int* n) 2440 { 2441 int ierr; 2442 2443 PetscFunctionBegin; 2444 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2445 PetscValidIntPointer(m); 2446 PetscValidIntPointer(n); 2447 if (!mat->ops.getownershiprange) SETERRQ(PETSC_ERR_SUP,0,""); 2448 ierr = (*mat->ops.getownershiprange)(mat,m,n);CHKERRQ(ierr); 2449 PetscFunctionReturn(0); 2450 } 2451 2452 #undef __FUNC__ 2453 #define __FUNC__ "MatILUFactorSymbolic" 2454 /*@ 2455 MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix. 2456 Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric() 2457 to complete the factorization. 2458 2459 Input Parameters: 2460 . mat - the matrix 2461 . row - row permutation 2462 . column - column permutation 2463 . fill - number of levels of fill 2464 . f - expected fill as ratio of the original number of nonzeros, 2465 for example 3.0; choosing this parameter well can result in 2466 more efficient use of time and space. Run your code with -log_info 2467 to determine an optimal value to use. 2468 2469 Output Parameters: 2470 . fact - new matrix that has been symbolically factored 2471 2472 Notes: 2473 See the file $(PETSC_DIR)/Performace for additional information about 2474 choosing the fill factor for better efficiency. 2475 2476 .keywords: matrix, factor, incomplete, ILU, symbolic, fill 2477 2478 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric() 2479 @*/ 2480 int MatILUFactorSymbolic(Mat mat,IS row,IS col,double f,int fill,Mat *fact) 2481 { 2482 int ierr; 2483 2484 PetscFunctionBegin; 2485 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2486 PetscValidPointer(fact); 2487 if (fill < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative"); 2488 if (!mat->ops.ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support parallel ILU"); 2489 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2490 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2491 2492 PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0); 2493 ierr = (*mat->ops.ilufactorsymbolic)(mat,row,col,f,fill,fact); CHKERRQ(ierr); 2494 PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0); 2495 PetscFunctionReturn(0); 2496 } 2497 2498 #undef __FUNC__ 2499 #define __FUNC__ "MatIncompleteCholeskyFactorSymbolic" 2500 /*@ 2501 MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete 2502 Cholesky factorization for a symmetric matrix. Use 2503 MatCholeskyFactorNumeric() to complete the factorization. 2504 2505 Input Parameters: 2506 . mat - the matrix 2507 . perm - row and column permutation 2508 . fill - levels of fill 2509 . f - expected fill as ratio of original fill 2510 2511 Output Parameter: 2512 . fact - the factored matrix 2513 2514 Note: Currently only no-fill factorization is supported. 2515 2516 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill 2517 2518 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor() 2519 @*/ 2520 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,double f,int fill,Mat *fact) 2521 { 2522 int ierr; 2523 2524 PetscFunctionBegin; 2525 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2526 PetscValidPointer(fact); 2527 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2528 if (fill < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative"); 2529 if (!mat->ops.incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel"); 2530 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2531 2532 PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0); 2533 ierr = (*mat->ops.incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr); 2534 PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0); 2535 PetscFunctionReturn(0); 2536 } 2537 2538 #undef __FUNC__ 2539 #define __FUNC__ "MatGetArray" 2540 /*@C 2541 MatGetArray - Returns a pointer to the element values in the matrix. 2542 This routine is implementation dependent, and may not even work for 2543 certain matrix types. You MUST call MatRestoreArray() when you no 2544 longer need to access the array. 2545 2546 Input Parameter: 2547 . mat - the matrix 2548 2549 Output Parameter: 2550 . v - the location of the values 2551 2552 Fortran Note: 2553 The Fortran interface is slightly different from that given below. 2554 See the Fortran chapter of the users manual and 2555 petsc/src/mat/examples for details. 2556 2557 .keywords: matrix, array, elements, values 2558 2559 .seealso: MatRestoreArray() 2560 @*/ 2561 int MatGetArray(Mat mat,Scalar **v) 2562 { 2563 int ierr; 2564 2565 PetscFunctionBegin; 2566 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2567 PetscValidPointer(v); 2568 if (!mat->ops.getarray) SETERRQ(PETSC_ERR_SUP,0,""); 2569 ierr = (*mat->ops.getarray)(mat,v);CHKERRQ(ierr); 2570 PetscFunctionReturn(0); 2571 } 2572 2573 #undef __FUNC__ 2574 #define __FUNC__ "MatRestoreArray" 2575 /*@C 2576 MatRestoreArray - Restores the matrix after MatGetArray has been called. 2577 2578 Input Parameter: 2579 . mat - the matrix 2580 . v - the location of the values 2581 2582 Fortran Note: 2583 The Fortran interface is slightly different from that given below. 2584 See the users manual and petsc/src/mat/examples for details. 2585 2586 .keywords: matrix, array, elements, values, restore 2587 2588 .seealso: MatGetArray() 2589 @*/ 2590 int MatRestoreArray(Mat mat,Scalar **v) 2591 { 2592 int ierr; 2593 2594 PetscFunctionBegin; 2595 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2596 PetscValidPointer(v); 2597 if (!mat->ops.restorearray) SETERRQ(PETSC_ERR_SUP,0,""); 2598 ierr = (*mat->ops.restorearray)(mat,v);CHKERRQ(ierr); 2599 PetscFunctionReturn(0); 2600 } 2601 2602 #undef __FUNC__ 2603 #define __FUNC__ "MatGetSubMatrices" 2604 /*@C 2605 MatGetSubMatrices - Extracts several submatrices from a matrix. If submat 2606 points to an array of valid matrices, they may be reused to store the new 2607 submatrices. 2608 2609 Input Parameters: 2610 . mat - the matrix 2611 . n - the number of submatrixes to be extracted 2612 . irow, icol - index sets of rows and columns to extract 2613 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 2614 2615 Output Parameter: 2616 . submat - the array of submatrices 2617 2618 Notes: 2619 MatGetSubMatrices() can extract only sequential submatrices 2620 (from both sequential and parallel matrices). Use MatGetSubMatrix() 2621 to extract a parallel submatrix. 2622 2623 When extracting submatrices from a parallel matrix, each processor can 2624 form a different submatrix by setting the rows and columns of its 2625 individual index sets according to the local submatrix desired. 2626 2627 When finished using the submatrices, the user should destroy 2628 them with MatDestroySubMatrices(). 2629 2630 .keywords: matrix, get, submatrix, submatrices 2631 2632 .seealso: MatDestroyMatrices(), MatGetSubMatrix() 2633 @*/ 2634 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **submat) 2635 { 2636 int ierr; 2637 2638 PetscFunctionBegin; 2639 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2640 if (!mat->ops.getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,""); 2641 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2642 2643 PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0); 2644 ierr = (*mat->ops.getsubmatrices)(mat,n,irow,icol,scall,submat); CHKERRQ(ierr); 2645 PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0); 2646 2647 PetscFunctionReturn(0); 2648 } 2649 2650 #undef __FUNC__ 2651 #define __FUNC__ "MatDestroyMatrices" 2652 /*@C 2653 MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices(). 2654 2655 Input Parameters: 2656 . n - the number of local matrices 2657 . mat - the matrices 2658 2659 .keywords: matrix, destroy, submatrix, submatrices 2660 2661 .seealso: MatGetSubMatrices() 2662 @*/ 2663 int MatDestroyMatrices(int n,Mat **mat) 2664 { 2665 int ierr,i; 2666 2667 PetscFunctionBegin; 2668 if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices"); 2669 PetscValidPointer(mat); 2670 for ( i=0; i<n; i++ ) { 2671 ierr = MatDestroy((*mat)[i]); CHKERRQ(ierr); 2672 } 2673 if (n) PetscFree(*mat); 2674 PetscFunctionReturn(0); 2675 } 2676 2677 #undef __FUNC__ 2678 #define __FUNC__ "MatIncreaseOverlap" 2679 /*@ 2680 MatIncreaseOverlap - Given a set of submatrices indicated by index sets, 2681 replaces the index sets by larger ones that represent submatrices with 2682 additional overlap. 2683 2684 Input Parameters: 2685 . mat - the matrix 2686 . n - the number of index sets 2687 . is - the array of pointers to index sets 2688 . ov - the additional overlap requested 2689 2690 .keywords: matrix, overlap, Schwarz 2691 2692 .seealso: MatGetSubMatrices() 2693 @*/ 2694 int MatIncreaseOverlap(Mat mat,int n, IS *is,int ov) 2695 { 2696 int ierr; 2697 2698 PetscFunctionBegin; 2699 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2700 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2701 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2702 2703 if (ov == 0) PetscFunctionReturn(0); 2704 if (!mat->ops.increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,""); 2705 PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0); 2706 ierr = (*mat->ops.increaseoverlap)(mat,n,is,ov); CHKERRQ(ierr); 2707 PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0); 2708 PetscFunctionReturn(0); 2709 } 2710 2711 #undef __FUNC__ 2712 #define __FUNC__ "MatPrintHelp" 2713 /*@ 2714 MatPrintHelp - Prints all the options for the matrix. 2715 2716 Input Parameter: 2717 . mat - the matrix 2718 2719 Options Database Keys: 2720 $ -help, -h 2721 2722 .keywords: mat, help 2723 2724 .seealso: MatCreate(), MatCreateXXX() 2725 @*/ 2726 int MatPrintHelp(Mat mat) 2727 { 2728 static int called = 0; 2729 int ierr; 2730 MPI_Comm comm; 2731 2732 PetscFunctionBegin; 2733 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2734 2735 comm = mat->comm; 2736 if (!called) { 2737 (*PetscHelpPrintf)(comm,"General matrix options:\n"); 2738 (*PetscHelpPrintf)(comm," -mat_view_info: view basic matrix info during MatAssemblyEnd()\n"); 2739 (*PetscHelpPrintf)(comm," -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n"); 2740 (*PetscHelpPrintf)(comm," -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n"); 2741 (*PetscHelpPrintf)(comm," -draw_pause <sec>: set seconds of display pause\n"); 2742 (*PetscHelpPrintf)(comm," -display <name>: set alternate display\n"); 2743 called = 1; 2744 } 2745 if (mat->ops.printhelp) { 2746 ierr = (*mat->ops.printhelp)(mat); CHKERRQ(ierr); 2747 } 2748 PetscFunctionReturn(0); 2749 } 2750 2751 #undef __FUNC__ 2752 #define __FUNC__ "MatGetBlockSize" 2753 /*@ 2754 MatGetBlockSize - Returns the matrix block size; useful especially for the 2755 block row and block diagonal formats. 2756 2757 Input Parameter: 2758 . mat - the matrix 2759 2760 Output Parameter: 2761 . bs - block size 2762 2763 Notes: 2764 $ block diagonal formats: MATSEQBDIAG, MATMPIBDIAG 2765 $ block row formats: MATSEQBAIJ, MATMPIBAIJ 2766 2767 .keywords: matrix, get, block, size 2768 2769 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag() 2770 @*/ 2771 int MatGetBlockSize(Mat mat,int *bs) 2772 { 2773 int ierr; 2774 2775 PetscFunctionBegin; 2776 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2777 PetscValidIntPointer(bs); 2778 if (!mat->ops.getblocksize) SETERRQ(PETSC_ERR_SUP,0,""); 2779 ierr = (*mat->ops.getblocksize)(mat,bs);CHKERRQ(ierr); 2780 PetscFunctionReturn(0); 2781 } 2782 2783 #undef __FUNC__ 2784 #define __FUNC__ "MatGetRowIJ" 2785 /*@C 2786 MatGetRowIJ - Returns the compress row storage i and j indices for sequential matrices. 2787 EXPERTS ONLY. 2788 2789 Input Parameters: 2790 . mat - the matrix 2791 . shift - 0 or 1 indicating we want the indices starting at 0 or 1 2792 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2793 symmetrized 2794 2795 Output Parameters: 2796 . n - number of rows and columns in the (possibly compressed) matrix 2797 . ia - the row indices 2798 . ja - the column indices 2799 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2800 @*/ 2801 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2802 { 2803 int ierr; 2804 2805 PetscFunctionBegin; 2806 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2807 if (ia) PetscValidIntPointer(ia); 2808 if (ja) PetscValidIntPointer(ja); 2809 PetscValidIntPointer(done); 2810 if (!mat->ops.getrowij) *done = PETSC_FALSE; 2811 else { 2812 *done = PETSC_TRUE; 2813 ierr = (*mat->ops.getrowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2814 } 2815 PetscFunctionReturn(0); 2816 } 2817 2818 #undef __FUNC__ 2819 #define __FUNC__ "MatGetColumnIJ" 2820 /*@C 2821 MatGetColumnIJ - Returns the compress Column storage i and j indices for sequential matrices. 2822 EXPERTS ONLY. 2823 2824 Input Parameters: 2825 . mat - the matrix 2826 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2827 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2828 symmetrized 2829 2830 Output Parameters: 2831 . n - number of Columns and columns in the (possibly compressed) matrix 2832 . ia - the Column indices 2833 . ja - the column indices 2834 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2835 @*/ 2836 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2837 { 2838 int ierr; 2839 2840 PetscFunctionBegin; 2841 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2842 if (ia) PetscValidIntPointer(ia); 2843 if (ja) PetscValidIntPointer(ja); 2844 PetscValidIntPointer(done); 2845 2846 if (!mat->ops.getcolumnij) *done = PETSC_FALSE; 2847 else { 2848 *done = PETSC_TRUE; 2849 ierr = (*mat->ops.getcolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2850 } 2851 PetscFunctionReturn(0); 2852 } 2853 2854 #undef __FUNC__ 2855 #define __FUNC__ "MatRestoreRowIJ" 2856 /*@C 2857 MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with 2858 MatGetRowIJ(). EXPERTS ONLY. 2859 2860 Input Parameters: 2861 . mat - the matrix 2862 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2863 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2864 symmetrized 2865 2866 Output Parameters: 2867 . n - size of (possibly compressed) matrix 2868 . ia - the row indices 2869 . ja - the column indices 2870 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2871 @*/ 2872 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2873 { 2874 int ierr; 2875 2876 PetscFunctionBegin; 2877 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2878 if (ia) PetscValidIntPointer(ia); 2879 if (ja) PetscValidIntPointer(ja); 2880 PetscValidIntPointer(done); 2881 2882 if (!mat->ops.restorerowij) *done = PETSC_FALSE; 2883 else { 2884 *done = PETSC_TRUE; 2885 ierr = (*mat->ops.restorerowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2886 } 2887 PetscFunctionReturn(0); 2888 } 2889 2890 #undef __FUNC__ 2891 #define __FUNC__ "MatRestoreColumnIJ" 2892 /*@C 2893 MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with 2894 MatGetColumnIJ(). EXPERTS ONLY. 2895 2896 Input Parameters: 2897 . mat - the matrix 2898 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2899 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2900 symmetrized 2901 2902 Output Parameters: 2903 . n - size of (possibly compressed) matrix 2904 . ia - the Column indices 2905 . ja - the column indices 2906 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2907 @*/ 2908 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2909 { 2910 int ierr; 2911 2912 PetscFunctionBegin; 2913 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2914 if (ia) PetscValidIntPointer(ia); 2915 if (ja) PetscValidIntPointer(ja); 2916 PetscValidIntPointer(done); 2917 2918 if (!mat->ops.restorecolumnij) *done = PETSC_FALSE; 2919 else { 2920 *done = PETSC_TRUE; 2921 ierr = (*mat->ops.restorecolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2922 } 2923 PetscFunctionReturn(0); 2924 } 2925 2926 #undef __FUNC__ 2927 #define __FUNC__ "MatColoringPatch" 2928 /*@C 2929 MatColoringPatch - EXPERTS ONLY, used inside matrix coloring routines that 2930 use matGetRowIJ() and/or MatGetColumnIJ(). 2931 2932 Input Parameters: 2933 . mat - the matrix 2934 . n - number of colors 2935 . colorarray - array indicating color for each column 2936 2937 Output Parameters: 2938 . iscoloring - coloring generated using colorarray information 2939 2940 @*/ 2941 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring) 2942 { 2943 int ierr; 2944 2945 PetscFunctionBegin; 2946 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2947 PetscValidIntPointer(colorarray); 2948 2949 if (!mat->ops.coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");} 2950 else { 2951 ierr = (*mat->ops.coloringpatch)(mat,n,colorarray,iscoloring); CHKERRQ(ierr); 2952 } 2953 PetscFunctionReturn(0); 2954 } 2955 2956 2957 #undef __FUNC__ 2958 #define __FUNC__ "MatSetUnfactored" 2959 /*@ 2960 MatSetUnfactored - Resets a factored matrix to be treated as unfactored. 2961 2962 Input Paramter: 2963 . mat - the factored matrix to be reset 2964 2965 Notes: 2966 This routine should be used only with factored matrices formed by in-place 2967 factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE 2968 format). This option can save memory, for example, when solving nonlinear 2969 systems with a matrix-free Newton-Krylov method and a matrix-based, in-place 2970 ILU(0) preconditioner. 2971 2972 Note that one can specify in-place ILU(0) factorization by calling 2973 $ PCType(pc,PCILU); 2974 $ PCILUSeUseInPlace(pc); 2975 or by using the options -pc_type ilu -pc_ilu_in_place 2976 2977 In-place factorization ILU(0) can also be used as a local 2978 solver for the blocks within the block Jacobi or additive Schwarz 2979 methods (runtime option: -sub_pc_ilu_in_place). See the discussion 2980 of these preconditioners in the users manual for details on setting 2981 local solver options. 2982 2983 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace() 2984 2985 .keywords: matrix-free, in-place ILU, in-place LU 2986 @*/ 2987 int MatSetUnfactored(Mat mat) 2988 { 2989 int ierr; 2990 2991 PetscFunctionBegin; 2992 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2993 mat->factor = 0; 2994 if (!mat->ops.setunfactored) PetscFunctionReturn(0); 2995 ierr = (*mat->ops.setunfactored)(mat); CHKERRQ(ierr); 2996 PetscFunctionReturn(0); 2997 } 2998 2999 #undef __FUNC__ 3000 #define __FUNC__ "MatGetType" 3001 /*@C 3002 MatGetType - Gets the matrix type and name (as a string) from the matrix. 3003 3004 Input Parameter: 3005 . mat - the matrix 3006 3007 Output Parameter: 3008 . type - the matrix type (or use PETSC_NULL) 3009 . name - name of matrix type (or use PETSC_NULL) 3010 3011 .keywords: matrix, get, type, name 3012 @*/ 3013 int MatGetType(Mat mat,MatType *type,char **name) 3014 { 3015 int itype = (int)mat->type; 3016 char *matname[10]; 3017 3018 PetscFunctionBegin; 3019 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3020 3021 if (type) *type = (MatType) mat->type; 3022 if (name) { 3023 /* Note: Be sure that this list corresponds to the enum in mat.h */ 3024 matname[0] = "MATSEQDENSE"; 3025 matname[1] = "MATSEQAIJ"; 3026 matname[2] = "MATMPIAIJ"; 3027 matname[3] = "MATSHELL"; 3028 matname[4] = "MATMPIROWBS"; 3029 matname[5] = "MATSEQBDIAG"; 3030 matname[6] = "MATMPIBDIAG"; 3031 matname[7] = "MATMPIDENSE"; 3032 matname[8] = "MATSEQBAIJ"; 3033 matname[9] = "MATMPIBAIJ"; 3034 3035 if (itype < 0 || itype > 9) *name = "Unknown matrix type"; 3036 else *name = matname[itype]; 3037 } 3038 PetscFunctionReturn(0); 3039 } 3040 3041 /*MC 3042 MatGetArrayF90 - Accesses a matrix array from Fortran90. 3043 3044 Input Parameter: 3045 . x - matrix 3046 3047 Output Parameter: 3048 . xx_v - the Fortran90 pointer to the array 3049 . ierr - error code 3050 3051 Synopsis: 3052 MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 3053 3054 Usage: 3055 $ Scalar, pointer xx_v(:) 3056 $ .... 3057 $ call MatGetArrayF90(x,xx_v,ierr) 3058 $ a = xx_v(3) 3059 $ call MatRestoreArrayF90(x,xx_v,ierr) 3060 3061 Notes: 3062 Currently only supported using the NAG F90 compiler. 3063 3064 .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray() 3065 3066 .keywords: matrix, array, f90 3067 M*/ 3068 3069 /*MC 3070 MatRestoreArrayF90 - Restores a matrix array that has been 3071 accessed with MatGetArrayF90(). 3072 3073 Input Parameters: 3074 . x - matrix 3075 . xx_v - the Fortran90 pointer to the array 3076 3077 Output Parameter: 3078 . ierr - error code 3079 3080 Synopsis: 3081 MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 3082 3083 Example of Usage: 3084 $ Scalar, pointer xx_v(:) 3085 $ .... 3086 $ call MatGetArrayF90(x,xx_v,ierr) 3087 $ a = xx_v(3) 3088 $ call MatRestoreArrayF90(x,xx_v,ierr) 3089 3090 Notes: 3091 Currently only supported using the NAG F90 compiler. 3092 3093 .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray() 3094 3095 .keywords: matrix, array, f90 3096 M*/ 3097 3098 3099 #undef __FUNC__ 3100 #define __FUNC__ "MatGetSubMatrix" 3101 /*@ 3102 MatGetSubMatrix - Gets a single submatrix on the same number of processors 3103 as the original matrix. 3104 3105 Input Parameters: 3106 . mat - the original matrix 3107 . isrow - rows this processor should obtain 3108 . iscol - columns for all processors you wish kept 3109 . csize - number of columns "local" to this processor (does nothing for sequential 3110 matrices). This should match the result from VecGetLocalSize() if you 3111 plan to use the matrix in a A*x 3112 . cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 3113 3114 Output Parameters: 3115 . newmat - the new submatrix, of the same type as the old 3116 3117 .seealso: MatGetSubMatrices() 3118 3119 @*/ 3120 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatGetSubMatrixCall cll,Mat *newmat) 3121 { 3122 int ierr, size; 3123 Mat *local; 3124 3125 PetscFunctionBegin; 3126 MPI_Comm_size(mat->comm,&size); 3127 3128 /* if original matrix is on just one processor then use submatrix generated */ 3129 if (size == 1 && cll == MAT_REUSE_MATRIX) { 3130 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr); 3131 PetscFunctionReturn(0); 3132 } else if (size == 1) { 3133 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr); 3134 *newmat = *local; 3135 PetscFree(local); 3136 PetscFunctionReturn(0); 3137 } 3138 3139 if (!mat->ops.getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented"); 3140 ierr = (*mat->ops.getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr); 3141 PetscFunctionReturn(0); 3142 } 3143 3144 3145 3146