1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: matrix.c,v 1.276 1998/01/28 20:32:08 bsmith Exp curfman $"; 3 #endif 4 5 /* 6 This is where the abstract matrix operations are defined 7 */ 8 9 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 10 #include "src/vec/vecimpl.h" 11 #include "pinclude/pviewer.h" 12 13 #undef __FUNC__ 14 #define __FUNC__ "MatGetRow" 15 /*@C 16 MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow() 17 for each row that you get to ensure that your application does 18 not bleed memory. 19 20 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 Note: probably really only in-place when level is zero. 978 .keywords: matrix, factor, ILU, in-place 979 980 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 981 @*/ 982 int MatILUFactor(Mat mat,IS row,IS col,double f,int level) 983 { 984 int ierr; 985 986 PetscFunctionBegin; 987 PetscValidHeaderSpecific(mat,MAT_COOKIE); 988 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square"); 989 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 990 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 991 if (!mat->ops.ilufactor) SETERRQ(PETSC_ERR_SUP,0,""); 992 993 PLogEventBegin(MAT_ILUFactor,mat,row,col,0); 994 ierr = (*mat->ops.ilufactor)(mat,row,col,f,level); CHKERRQ(ierr); 995 PLogEventEnd(MAT_ILUFactor,mat,row,col,0); 996 PetscFunctionReturn(0); 997 } 998 999 #undef __FUNC__ 1000 #define __FUNC__ "MatLUFactorSymbolic" 1001 /*@ 1002 MatLUFactorSymbolic - Performs symbolic LU factorization of matrix. 1003 Call this routine before calling MatLUFactorNumeric(). 1004 1005 Input Parameters: 1006 . mat - the matrix 1007 . row, col - row and column permutations 1008 . f - expected fill as ratio of the original number of nonzeros, 1009 for example 3.0; choosing this parameter well can result in 1010 more efficient use of time and space. Run with the option -log_info 1011 to determine an optimal value to use 1012 1013 Output Parameter: 1014 . fact - new matrix that has been symbolically factored 1015 1016 Notes: 1017 See the file $(PETSC_DIR)/Performace for additional information about 1018 choosing the fill factor for better efficiency. 1019 1020 .keywords: matrix, factor, LU, symbolic, fill 1021 1022 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor() 1023 @*/ 1024 int MatLUFactorSymbolic(Mat mat,IS row,IS col,double f,Mat *fact) 1025 { 1026 int ierr; 1027 1028 PetscFunctionBegin; 1029 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1030 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square"); 1031 PetscValidPointer(fact); 1032 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1033 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1034 if (!mat->ops.lufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,""); 1035 1036 PLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0); 1037 ierr = (*mat->ops.lufactorsymbolic)(mat,row,col,f,fact); CHKERRQ(ierr); 1038 PLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0); 1039 PetscFunctionReturn(0); 1040 } 1041 1042 #undef __FUNC__ 1043 #define __FUNC__ "MatLUFactorNumeric" 1044 /*@ 1045 MatLUFactorNumeric - Performs numeric LU factorization of a matrix. 1046 Call this routine after first calling MatLUFactorSymbolic(). 1047 1048 Input Parameters: 1049 . mat - the matrix 1050 . row, col - row and column permutations 1051 1052 Output Parameters: 1053 . fact - symbolically factored matrix that must have been generated 1054 by MatLUFactorSymbolic() 1055 1056 Notes: 1057 See MatLUFactor() for in-place factorization. See 1058 MatCholeskyFactorNumeric() for the symmetric, positive definite case. 1059 1060 .keywords: matrix, factor, LU, numeric 1061 1062 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor() 1063 @*/ 1064 int MatLUFactorNumeric(Mat mat,Mat *fact) 1065 { 1066 int ierr,flg; 1067 1068 PetscFunctionBegin; 1069 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1070 PetscValidPointer(fact); 1071 PetscValidHeaderSpecific(*fact,MAT_COOKIE); 1072 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1073 if (mat->M != (*fact)->M || mat->N != (*fact)->N) 1074 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dim"); 1075 if (!(*fact)->ops.lufactornumeric) SETERRQ(PETSC_ERR_SUP,0,""); 1076 1077 PLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0); 1078 ierr = (*(*fact)->ops.lufactornumeric)(mat,fact); CHKERRQ(ierr); 1079 PLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0); 1080 ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg); CHKERRQ(ierr); 1081 if (flg) { 1082 ierr = OptionsHasName(0,"-mat_view_contour",&flg); CHKERRQ(ierr); 1083 if (flg) { 1084 ViewerPushFormat(VIEWER_DRAWX_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr); 1085 } 1086 ierr = MatView(*fact,VIEWER_DRAWX_(mat->comm)); CHKERRQ(ierr); 1087 ierr = ViewerFlush(VIEWER_DRAWX_(mat->comm)); CHKERRQ(ierr); 1088 if (flg) { 1089 ViewerPopFormat(VIEWER_DRAWX_(mat->comm));CHKERRQ(ierr); 1090 } 1091 } 1092 PetscFunctionReturn(0); 1093 } 1094 1095 #undef __FUNC__ 1096 #define __FUNC__ "MatCholeskyFactor" 1097 /*@ 1098 MatCholeskyFactor - Performs in-place Cholesky factorization of a 1099 symmetric matrix. 1100 1101 Input Parameters: 1102 . mat - the matrix 1103 . perm - row and column permutations 1104 . f - expected fill as ratio of original fill 1105 1106 Notes: 1107 See MatLUFactor() for the nonsymmetric case. See also 1108 MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric(). 1109 1110 .keywords: matrix, factor, in-place, Cholesky 1111 1112 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric() 1113 @*/ 1114 int MatCholeskyFactor(Mat mat,IS perm,double f) 1115 { 1116 int ierr; 1117 1118 PetscFunctionBegin; 1119 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1120 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square"); 1121 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1122 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1123 if (!mat->ops.choleskyfactor) SETERRQ(PETSC_ERR_SUP,0,""); 1124 1125 PLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0); 1126 ierr = (*mat->ops.choleskyfactor)(mat,perm,f); CHKERRQ(ierr); 1127 PLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0); 1128 PetscFunctionReturn(0); 1129 } 1130 1131 #undef __FUNC__ 1132 #define __FUNC__ "MatCholeskyFactorSymbolic" 1133 /*@ 1134 MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization 1135 of a symmetric matrix. 1136 1137 Input Parameters: 1138 . mat - the matrix 1139 . perm - row and column permutations 1140 . f - expected fill as ratio of original 1141 1142 Output Parameter: 1143 . fact - the factored matrix 1144 1145 Notes: 1146 See MatLUFactorSymbolic() for the nonsymmetric case. See also 1147 MatCholeskyFactor() and MatCholeskyFactorNumeric(). 1148 1149 .keywords: matrix, factor, factorization, symbolic, Cholesky 1150 1151 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric() 1152 @*/ 1153 int MatCholeskyFactorSymbolic(Mat mat,IS perm,double f,Mat *fact) 1154 { 1155 int ierr; 1156 1157 PetscFunctionBegin; 1158 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1159 PetscValidPointer(fact); 1160 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,0,"matrix must be square"); 1161 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1162 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1163 if (!mat->ops.choleskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,""); 1164 1165 PLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0); 1166 ierr = (*mat->ops.choleskyfactorsymbolic)(mat,perm,f,fact); CHKERRQ(ierr); 1167 PLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0); 1168 PetscFunctionReturn(0); 1169 } 1170 1171 #undef __FUNC__ 1172 #define __FUNC__ "MatCholeskyFactorNumeric" 1173 /*@ 1174 MatCholeskyFactorNumeric - Performs numeric Cholesky factorization 1175 of a symmetric matrix. Call this routine after first calling 1176 MatCholeskyFactorSymbolic(). 1177 1178 Input Parameter: 1179 . mat - the initial matrix 1180 1181 Output Parameter: 1182 . fact - the factored matrix 1183 1184 .keywords: matrix, factor, numeric, Cholesky 1185 1186 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric() 1187 @*/ 1188 int MatCholeskyFactorNumeric(Mat mat,Mat *fact) 1189 { 1190 int ierr; 1191 1192 PetscFunctionBegin; 1193 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1194 PetscValidPointer(fact); 1195 if (!mat->ops.choleskyfactornumeric) SETERRQ(PETSC_ERR_SUP,0,""); 1196 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1197 if (mat->M != (*fact)->M || mat->N != (*fact)->N) 1198 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dim"); 1199 1200 PLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0); 1201 ierr = (*mat->ops.choleskyfactornumeric)(mat,fact); CHKERRQ(ierr); 1202 PLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0); 1203 PetscFunctionReturn(0); 1204 } 1205 1206 /* ----------------------------------------------------------------*/ 1207 #undef __FUNC__ 1208 #define __FUNC__ "MatSolve" 1209 /*@ 1210 MatSolve - Solves A x = b, given a factored matrix. 1211 1212 Input Parameters: 1213 . mat - the factored matrix 1214 . b - the right-hand-side vector 1215 1216 Output Parameter: 1217 . x - the result vector 1218 1219 Notes: 1220 The vectors b and x cannot be the same. I.e., one cannot 1221 call MatSolve(A,x,x). 1222 1223 .keywords: matrix, linear system, solve, LU, Cholesky, triangular solve 1224 1225 .seealso: MatSolveAdd(), MatSolveTrans(), MatSolveTransAdd() 1226 @*/ 1227 int MatSolve(Mat mat,Vec b,Vec x) 1228 { 1229 int ierr; 1230 1231 PetscFunctionBegin; 1232 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1233 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1234 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1235 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1236 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1237 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1238 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1239 if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0); 1240 1241 if (!mat->ops.solve) SETERRQ(PETSC_ERR_SUP,0,""); 1242 PLogEventBegin(MAT_Solve,mat,b,x,0); 1243 ierr = (*mat->ops.solve)(mat,b,x); CHKERRQ(ierr); 1244 PLogEventEnd(MAT_Solve,mat,b,x,0); 1245 PetscFunctionReturn(0); 1246 } 1247 1248 #undef __FUNC__ 1249 #define __FUNC__ "MatForwardSolve" 1250 /* @ 1251 MatForwardSolve - Solves L x = b, given a factored matrix, A = LU. 1252 1253 Input Parameters: 1254 . mat - the factored matrix 1255 . b - the right-hand-side vector 1256 1257 Output Parameter: 1258 . x - the result vector 1259 1260 Notes: 1261 MatSolve() should be used for most applications, as it performs 1262 a forward solve followed by a backward solve. 1263 1264 The vectors b and x cannot be the same. I.e., one cannot 1265 call MatForwardSolve(A,x,x). 1266 1267 .keywords: matrix, forward, LU, Cholesky, triangular solve 1268 1269 .seealso: MatSolve(), MatBackwardSolve() 1270 @ */ 1271 int MatForwardSolve(Mat mat,Vec b,Vec x) 1272 { 1273 int ierr; 1274 1275 PetscFunctionBegin; 1276 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1277 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1278 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1279 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1280 if (!mat->ops.forwardsolve) SETERRQ(PETSC_ERR_SUP,0,""); 1281 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1282 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1283 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1284 1285 PLogEventBegin(MAT_ForwardSolve,mat,b,x,0); 1286 ierr = (*mat->ops.forwardsolve)(mat,b,x); CHKERRQ(ierr); 1287 PLogEventEnd(MAT_ForwardSolve,mat,b,x,0); 1288 PetscFunctionReturn(0); 1289 } 1290 1291 #undef __FUNC__ 1292 #define __FUNC__ "MatBackwardSolve" 1293 /* @ 1294 MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU. 1295 1296 Input Parameters: 1297 . mat - the factored matrix 1298 . b - the right-hand-side vector 1299 1300 Output Parameter: 1301 . x - the result vector 1302 1303 Notes: 1304 MatSolve() should be used for most applications, as it performs 1305 a forward solve followed by a backward solve. 1306 1307 The vectors b and x cannot be the same. I.e., one cannot 1308 call MatBackwardSolve(A,x,x). 1309 1310 .keywords: matrix, backward, LU, Cholesky, triangular solve 1311 1312 .seealso: MatSolve(), MatForwardSolve() 1313 @ */ 1314 int MatBackwardSolve(Mat mat,Vec b,Vec x) 1315 { 1316 int ierr; 1317 1318 PetscFunctionBegin; 1319 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1320 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1321 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1322 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1323 if (!mat->ops.backwardsolve) SETERRQ(PETSC_ERR_SUP,0,""); 1324 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1325 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1326 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1327 1328 PLogEventBegin(MAT_BackwardSolve,mat,b,x,0); 1329 ierr = (*mat->ops.backwardsolve)(mat,b,x); CHKERRQ(ierr); 1330 PLogEventEnd(MAT_BackwardSolve,mat,b,x,0); 1331 PetscFunctionReturn(0); 1332 } 1333 1334 #undef __FUNC__ 1335 #define __FUNC__ "MatSolveAdd" 1336 /*@ 1337 MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix. 1338 1339 Input Parameters: 1340 . mat - the factored matrix 1341 . b - the right-hand-side vector 1342 . y - the vector to be added to 1343 1344 Output Parameter: 1345 . x - the result vector 1346 1347 Notes: 1348 The vectors b and x cannot be the same. I.e., one cannot 1349 call MatSolveAdd(A,x,y,x). 1350 1351 .keywords: matrix, linear system, solve, LU, Cholesky, add 1352 1353 .seealso: MatSolve(), MatSolveTrans(), MatSolveTransAdd() 1354 @*/ 1355 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x) 1356 { 1357 Scalar one = 1.0; 1358 Vec tmp; 1359 int ierr; 1360 1361 PetscFunctionBegin; 1362 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1363 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1364 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1365 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1366 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1367 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1368 if (mat->M != y->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim"); 1369 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1370 if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim"); 1371 1372 PLogEventBegin(MAT_SolveAdd,mat,b,x,y); 1373 if (mat->ops.solveadd) { 1374 ierr = (*mat->ops.solveadd)(mat,b,y,x); CHKERRQ(ierr); 1375 } 1376 else { 1377 /* do the solve then the add manually */ 1378 if (x != y) { 1379 ierr = MatSolve(mat,b,x); CHKERRQ(ierr); 1380 ierr = VecAXPY(&one,y,x); CHKERRQ(ierr); 1381 } else { 1382 ierr = VecDuplicate(x,&tmp); CHKERRQ(ierr); 1383 PLogObjectParent(mat,tmp); 1384 ierr = VecCopy(x,tmp); CHKERRQ(ierr); 1385 ierr = MatSolve(mat,b,x); CHKERRQ(ierr); 1386 ierr = VecAXPY(&one,tmp,x); CHKERRQ(ierr); 1387 ierr = VecDestroy(tmp); CHKERRQ(ierr); 1388 } 1389 } 1390 PLogEventEnd(MAT_SolveAdd,mat,b,x,y); 1391 PetscFunctionReturn(0); 1392 } 1393 1394 #undef __FUNC__ 1395 #define __FUNC__ "MatSolveTrans" 1396 /*@ 1397 MatSolveTrans - Solves A' x = b, given a factored matrix. 1398 1399 Input Parameters: 1400 . mat - the factored matrix 1401 . b - the right-hand-side vector 1402 1403 Output Parameter: 1404 . x - the result vector 1405 1406 Notes: 1407 The vectors b and x cannot be the same. I.e., one cannot 1408 call MatSolveTrans(A,x,x). 1409 1410 .keywords: matrix, linear system, solve, LU, Cholesky, transpose 1411 1412 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransAdd() 1413 @*/ 1414 int MatSolveTrans(Mat mat,Vec b,Vec x) 1415 { 1416 int ierr; 1417 1418 PetscFunctionBegin; 1419 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1420 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1421 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1422 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1423 if (!mat->ops.solvetrans) SETERRQ(PETSC_ERR_SUP,0,""); 1424 if (mat->M != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1425 if (mat->N != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1426 1427 PLogEventBegin(MAT_SolveTrans,mat,b,x,0); 1428 ierr = (*mat->ops.solvetrans)(mat,b,x); CHKERRQ(ierr); 1429 PLogEventEnd(MAT_SolveTrans,mat,b,x,0); 1430 PetscFunctionReturn(0); 1431 } 1432 1433 #undef __FUNC__ 1434 #define __FUNC__ "MatSolveTransAdd" 1435 /*@ 1436 MatSolveTransAdd - Computes x = y + inv(trans(A)) b, given a 1437 factored matrix. 1438 1439 Input Parameters: 1440 . mat - the factored matrix 1441 . b - the right-hand-side vector 1442 . y - the vector to be added to 1443 1444 Output Parameter: 1445 . x - the result vector 1446 1447 Notes: 1448 The vectors b and x cannot be the same. I.e., one cannot 1449 call MatSolveTransAdd(A,x,y,x). 1450 1451 .keywords: matrix, linear system, solve, LU, Cholesky, transpose, add 1452 1453 .seealso: MatSolve(), MatSolveAdd(), MatSolveTrans() 1454 @*/ 1455 int MatSolveTransAdd(Mat mat,Vec b,Vec y,Vec x) 1456 { 1457 Scalar one = 1.0; 1458 int ierr; 1459 Vec tmp; 1460 1461 PetscFunctionBegin; 1462 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1463 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1464 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1465 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1466 if (mat->M != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1467 if (mat->N != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1468 if (mat->N != y->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim"); 1469 if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim"); 1470 1471 PLogEventBegin(MAT_SolveTransAdd,mat,b,x,y); 1472 if (mat->ops.solvetransadd) { 1473 ierr = (*mat->ops.solvetransadd)(mat,b,y,x); CHKERRQ(ierr); 1474 } 1475 else { 1476 /* do the solve then the add manually */ 1477 if (x != y) { 1478 ierr = MatSolveTrans(mat,b,x); CHKERRQ(ierr); 1479 ierr = VecAXPY(&one,y,x); CHKERRQ(ierr); 1480 } 1481 else { 1482 ierr = VecDuplicate(x,&tmp); CHKERRQ(ierr); 1483 PLogObjectParent(mat,tmp); 1484 ierr = VecCopy(x,tmp); CHKERRQ(ierr); 1485 ierr = MatSolveTrans(mat,b,x); CHKERRQ(ierr); 1486 ierr = VecAXPY(&one,tmp,x); CHKERRQ(ierr); 1487 ierr = VecDestroy(tmp); CHKERRQ(ierr); 1488 } 1489 } 1490 PLogEventEnd(MAT_SolveTransAdd,mat,b,x,y); 1491 PetscFunctionReturn(0); 1492 } 1493 /* ----------------------------------------------------------------*/ 1494 1495 #undef __FUNC__ 1496 #define __FUNC__ "MatRelax" 1497 /*@ 1498 MatRelax - Computes one relaxation sweep. 1499 1500 Input Parameters: 1501 . mat - the matrix 1502 . b - the right hand side 1503 . omega - the relaxation factor 1504 . flag - flag indicating the type of SOR, one of 1505 $ SOR_FORWARD_SWEEP 1506 $ SOR_BACKWARD_SWEEP 1507 $ SOR_SYMMETRIC_SWEEP (SSOR method) 1508 $ SOR_LOCAL_FORWARD_SWEEP 1509 $ SOR_LOCAL_BACKWARD_SWEEP 1510 $ SOR_LOCAL_SYMMETRIC_SWEEP (local SSOR) 1511 $ SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies 1512 $ upper/lower triangular part of matrix to 1513 $ vector (with omega) 1514 $ SOR_ZERO_INITIAL_GUESS - zero initial guess 1515 . shift - diagonal shift 1516 . its - the number of iterations 1517 1518 Output Parameters: 1519 . x - the solution (can contain an initial guess) 1520 1521 Notes: 1522 SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and 1523 SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings 1524 on each processor. 1525 1526 Application programmers will not generally use MatRelax() directly, 1527 but instead will employ the SLES/PC interface. 1528 1529 Notes for Advanced Users: 1530 The flags are implemented as bitwise inclusive or operations. 1531 For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP) 1532 to specify a zero initial guess for SSOR. 1533 1534 .keywords: matrix, relax, relaxation, sweep 1535 @*/ 1536 int MatRelax(Mat mat,Vec b,double omega,MatSORType flag,double shift, 1537 int its,Vec x) 1538 { 1539 int ierr; 1540 1541 PetscFunctionBegin; 1542 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1543 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1544 if (!mat->ops.relax) SETERRQ(PETSC_ERR_SUP,0,""); 1545 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1546 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1547 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1548 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1549 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1550 1551 PLogEventBegin(MAT_Relax,mat,b,x,0); 1552 ierr =(*mat->ops.relax)(mat,b,omega,flag,shift,its,x); CHKERRQ(ierr); 1553 PLogEventEnd(MAT_Relax,mat,b,x,0); 1554 PetscFunctionReturn(0); 1555 } 1556 1557 #undef __FUNC__ 1558 #define __FUNC__ "MatCopy_Basic" 1559 /* 1560 Default matrix copy routine. 1561 */ 1562 int MatCopy_Basic(Mat A,Mat B) 1563 { 1564 int ierr,i,rstart,rend,nz,*cwork; 1565 Scalar *vwork; 1566 1567 PetscFunctionBegin; 1568 ierr = MatZeroEntries(B); CHKERRQ(ierr); 1569 ierr = MatGetOwnershipRange(A,&rstart,&rend); CHKERRQ(ierr); 1570 for (i=rstart; i<rend; i++) { 1571 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1572 ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 1573 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1574 } 1575 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1576 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1577 PetscFunctionReturn(0); 1578 } 1579 1580 #undef __FUNC__ 1581 #define __FUNC__ "MatCopy" 1582 /*@C 1583 MatCopy - Copys a matrix to another matrix. 1584 1585 Input Parameters: 1586 . A - the matrix 1587 1588 Output Parameter: 1589 . B - where the copy is put 1590 1591 Notes: 1592 MatCopy() copies the matrix entries of a matrix to another existing 1593 matrix (after first zeroing the second matrix). A related routine is 1594 MatConvert(), which first creates a new matrix and then copies the data. 1595 1596 .keywords: matrix, copy, convert 1597 1598 .seealso: MatConvert() 1599 @*/ 1600 int MatCopy(Mat A,Mat B) 1601 { 1602 int ierr; 1603 1604 PetscFunctionBegin; 1605 PetscValidHeaderSpecific(A,MAT_COOKIE); PetscValidHeaderSpecific(B,MAT_COOKIE); 1606 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1607 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1608 if (A->M != B->M || A->N != B->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim"); 1609 1610 PLogEventBegin(MAT_Copy,A,B,0,0); 1611 if (A->ops.copy) { 1612 ierr = (*A->ops.copy)(A,B); CHKERRQ(ierr); 1613 } 1614 else { /* generic conversion */ 1615 ierr = MatCopy_Basic(A,B); CHKERRQ(ierr); 1616 } 1617 PLogEventEnd(MAT_Copy,A,B,0,0); 1618 PetscFunctionReturn(0); 1619 } 1620 1621 static int MatConvertersSet = 0; 1622 static int (*MatConverters[MAX_MATRIX_TYPES][MAX_MATRIX_TYPES])(Mat,MatType,Mat*) = 1623 {{0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1624 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 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 1630 #undef __FUNC__ 1631 #define __FUNC__ "MatConvertRegister" 1632 /*@C 1633 MatConvertRegister - Allows one to register a routine that converts between 1634 two matrix types. 1635 1636 Input Parameters: 1637 . intype - the type of matrix (defined in include/mat.h), for example, MATSEQAIJ. 1638 . outtype - new matrix type, or MATSAME 1639 1640 .seealso: MatConvertRegisterAll() 1641 @*/ 1642 int MatConvertRegister(MatType intype,MatType outtype,int (*converter)(Mat,MatType,Mat*)) 1643 { 1644 PetscFunctionBegin; 1645 MatConverters[intype][outtype] = converter; 1646 MatConvertersSet = 1; 1647 PetscFunctionReturn(0); 1648 } 1649 1650 #undef __FUNC__ 1651 #define __FUNC__ "MatConvert" 1652 /*@C 1653 MatConvert - Converts a matrix to another matrix, either of the same 1654 or different type. 1655 1656 Input Parameters: 1657 . mat - the matrix 1658 . newtype - new matrix type. Use MATSAME to create a new matrix of the 1659 same type as the original matrix. 1660 1661 Output Parameter: 1662 . M - pointer to place new matrix 1663 1664 Notes: 1665 MatConvert() first creates a new matrix and then copies the data from 1666 the first matrix. A related routine is MatCopy(), which copies the matrix 1667 entries of one matrix to another already existing matrix context. 1668 1669 .keywords: matrix, copy, convert 1670 1671 .seealso: MatCopy(), MatDuplicate() 1672 @*/ 1673 int MatConvert(Mat mat,MatType newtype,Mat *M) 1674 { 1675 int ierr; 1676 1677 PetscFunctionBegin; 1678 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1679 PetscValidPointer(M); 1680 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1681 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1682 1683 if (newtype > MAX_MATRIX_TYPES || newtype < -1) { 1684 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Not a valid matrix type"); 1685 } 1686 *M = 0; 1687 1688 if (!MatConvertersSet) { 1689 ierr = MatLoadRegisterAll(); CHKERRQ(ierr); 1690 } 1691 1692 PLogEventBegin(MAT_Convert,mat,0,0,0); 1693 if ((newtype == mat->type || newtype == MATSAME) && mat->ops.convertsametype) { 1694 ierr = (*mat->ops.convertsametype)(mat,M,COPY_VALUES); CHKERRQ(ierr); 1695 } else { 1696 if (!MatConvertersSet) { 1697 ierr = MatConvertRegisterAll(); CHKERRQ(ierr); 1698 } 1699 if (!MatConverters[mat->type][newtype]) { 1700 SETERRQ(PETSC_ERR_ARG_WRONG,1,"Invalid matrix type, or matrix converter not registered"); 1701 } 1702 ierr = (*MatConverters[mat->type][newtype])(mat,newtype,M); CHKERRQ(ierr); 1703 } 1704 PLogEventEnd(MAT_Convert,mat,0,0,0); 1705 PetscFunctionReturn(0); 1706 } 1707 1708 #undef __FUNC__ 1709 #define __FUNC__ "MatDuplicate" 1710 /*@C 1711 MatDuplicate - Duplicates a matrix including the non-zero structure, but 1712 does not copy over the values. 1713 1714 Input Parameters: 1715 . mat - the matrix 1716 1717 Output Parameter: 1718 . M - pointer to place new matrix 1719 1720 .keywords: matrix, copy, convert, duplicate 1721 1722 .seealso: MatCopy(), MatDuplicate(), MatConvert() 1723 @*/ 1724 int MatDuplicate(Mat mat,Mat *M) 1725 { 1726 int ierr; 1727 1728 PetscFunctionBegin; 1729 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1730 PetscValidPointer(M); 1731 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1732 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1733 1734 *M = 0; 1735 PLogEventBegin(MAT_Convert,mat,0,0,0); 1736 if (!mat->ops.convertsametype) { 1737 SETERRQ(PETSC_ERR_SUP,1,"Not written for this matrix type"); 1738 } 1739 ierr = (*mat->ops.convertsametype)(mat,M,DO_NOT_COPY_VALUES); CHKERRQ(ierr); 1740 PLogEventEnd(MAT_Convert,mat,0,0,0); 1741 PetscFunctionReturn(0); 1742 } 1743 1744 #undef __FUNC__ 1745 #define __FUNC__ "MatGetDiagonal" 1746 /*@ 1747 MatGetDiagonal - Gets the diagonal of a matrix. 1748 1749 Input Parameters: 1750 . mat - the matrix 1751 . v - the vector for storing the diagonal 1752 1753 Output Parameter: 1754 . v - the diagonal of the matrix 1755 1756 Notes: For the SeqAIJ matrix format, this routine may also be called 1757 on a LU factored matrix; in that case it routines the reciprocal of 1758 the diagonal entries in U. It returns the entries permuted by the 1759 row and column permutation used during the symbolic factorization. 1760 1761 .keywords: matrix, get, diagonal 1762 @*/ 1763 int MatGetDiagonal(Mat mat,Vec v) 1764 { 1765 int ierr; 1766 1767 PetscFunctionBegin; 1768 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v,VEC_COOKIE); 1769 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1770 /* 1771 The error checking for a factored matrix is handled inside 1772 each matrix type, since MatGetDiagonal() is supported by 1773 factored AIJ matrices 1774 */ 1775 /* if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); */ 1776 if (!mat->ops.getdiagonal) SETERRQ(PETSC_ERR_SUP,0,""); 1777 ierr = (*mat->ops.getdiagonal)(mat,v);CHKERRQ(ierr); 1778 PetscFunctionReturn(0); 1779 } 1780 1781 #undef __FUNC__ 1782 #define __FUNC__ "MatTranspose" 1783 /*@C 1784 MatTranspose - Computes an in-place or out-of-place transpose of a matrix. 1785 1786 Input Parameter: 1787 . mat - the matrix to transpose 1788 1789 Output Parameters: 1790 . B - the transpose (or pass in PETSC_NULL for an in-place transpose) 1791 1792 .keywords: matrix, transpose 1793 1794 .seealso: MatMultTrans(), MatMultTransAdd() 1795 @*/ 1796 int MatTranspose(Mat mat,Mat *B) 1797 { 1798 int ierr; 1799 1800 PetscFunctionBegin; 1801 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1802 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1803 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1804 if (!mat->ops.transpose) SETERRQ(PETSC_ERR_SUP,0,""); 1805 ierr = (*mat->ops.transpose)(mat,B);CHKERRQ(ierr); 1806 PetscFunctionReturn(0); 1807 } 1808 1809 #undef __FUNC__ 1810 #define __FUNC__ "MatPermute" 1811 /*@C 1812 MatPermute - Creates a new matrix with rows and columns permuted from the 1813 original. 1814 1815 Input Parameter: 1816 . mat - the matrix to permute 1817 . row - row permutation 1818 . col - column permutation 1819 1820 Output Parameters: 1821 . B - the permuted matrix 1822 1823 .keywords: matrix, transpose 1824 1825 .seealso: MatGetReordering() 1826 @*/ 1827 int MatPermute(Mat mat,IS row,IS col,Mat *B) 1828 { 1829 int ierr; 1830 1831 PetscFunctionBegin; 1832 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1833 PetscValidHeaderSpecific(row,IS_COOKIE); 1834 PetscValidHeaderSpecific(col,IS_COOKIE); 1835 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1836 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1837 if (!mat->ops.permute) SETERRQ(PETSC_ERR_SUP,0,""); 1838 ierr = (*mat->ops.permute)(mat,row,col,B);CHKERRQ(ierr); 1839 PetscFunctionReturn(0); 1840 } 1841 1842 #undef __FUNC__ 1843 #define __FUNC__ "MatEqual" 1844 /*@ 1845 MatEqual - Compares two matrices. 1846 1847 Input Parameters: 1848 . A - the first matrix 1849 . B - the second matrix 1850 1851 Output Parameter: 1852 . flg : PETSC_TRUE if the matrices are equal; 1853 PETSC_FALSE otherwise. 1854 1855 .keywords: matrix, equal, equivalent 1856 @*/ 1857 int MatEqual(Mat A,Mat B,PetscTruth *flg) 1858 { 1859 int ierr; 1860 1861 PetscFunctionBegin; 1862 PetscValidHeaderSpecific(A,MAT_COOKIE); PetscValidHeaderSpecific(B,MAT_COOKIE); 1863 PetscValidIntPointer(flg); 1864 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1865 if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1866 if (A->M != B->M || A->N != B->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim"); 1867 if (!A->ops.equal) SETERRQ(PETSC_ERR_SUP,0,""); 1868 ierr = (*A->ops.equal)(A,B,flg);CHKERRQ(ierr); 1869 PetscFunctionReturn(0); 1870 } 1871 1872 #undef __FUNC__ 1873 #define __FUNC__ "MatDiagonalScale" 1874 /*@ 1875 MatDiagonalScale - Scales a matrix on the left and right by diagonal 1876 matrices that are stored as vectors. Either of the two scaling 1877 matrices can be PETSC_NULL. 1878 1879 Input Parameters: 1880 . mat - the matrix to be scaled 1881 . l - the left scaling vector (or PETSC_NULL) 1882 . r - the right scaling vector (or PETSC_NULL) 1883 1884 Notes: 1885 MatDiagonalScale() computes A <- LAR, where 1886 $ L = a diagonal matrix 1887 $ R = a diagonal matrix 1888 1889 .keywords: matrix, diagonal, scale 1890 1891 .seealso: MatDiagonalScale() 1892 @*/ 1893 int MatDiagonalScale(Mat mat,Vec l,Vec r) 1894 { 1895 int ierr; 1896 1897 PetscFunctionBegin; 1898 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1899 if (!mat->ops.diagonalscale) SETERRQ(PETSC_ERR_SUP,0,""); 1900 if (l) PetscValidHeaderSpecific(l,VEC_COOKIE); 1901 if (r) PetscValidHeaderSpecific(r,VEC_COOKIE); 1902 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1903 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1904 1905 PLogEventBegin(MAT_Scale,mat,0,0,0); 1906 ierr = (*mat->ops.diagonalscale)(mat,l,r); CHKERRQ(ierr); 1907 PLogEventEnd(MAT_Scale,mat,0,0,0); 1908 PetscFunctionReturn(0); 1909 } 1910 1911 #undef __FUNC__ 1912 #define __FUNC__ "MatScale" 1913 /*@ 1914 MatScale - Scales all elements of a matrix by a given number. 1915 1916 Input Parameters: 1917 . mat - the matrix to be scaled 1918 . a - the scaling value 1919 1920 Output Parameter: 1921 . mat - the scaled matrix 1922 1923 .keywords: matrix, scale 1924 1925 .seealso: MatDiagonalScale() 1926 @*/ 1927 int MatScale(Scalar *a,Mat mat) 1928 { 1929 int ierr; 1930 1931 PetscFunctionBegin; 1932 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1933 PetscValidScalarPointer(a); 1934 if (!mat->ops.scale) SETERRQ(PETSC_ERR_SUP,0,""); 1935 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1936 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1937 1938 PLogEventBegin(MAT_Scale,mat,0,0,0); 1939 ierr = (*mat->ops.scale)(a,mat); CHKERRQ(ierr); 1940 PLogEventEnd(MAT_Scale,mat,0,0,0); 1941 PetscFunctionReturn(0); 1942 } 1943 1944 #undef __FUNC__ 1945 #define __FUNC__ "MatNorm" 1946 /*@ 1947 MatNorm - Calculates various norms of a matrix. 1948 1949 Input Parameters: 1950 . mat - the matrix 1951 . type - the type of norm, NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY 1952 1953 Output Parameters: 1954 . norm - the resulting norm 1955 1956 .keywords: matrix, norm, Frobenius 1957 @*/ 1958 int MatNorm(Mat mat,NormType type,double *norm) 1959 { 1960 int ierr; 1961 1962 PetscFunctionBegin; 1963 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1964 PetscValidScalarPointer(norm); 1965 1966 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1967 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1968 if (!mat->ops.norm) SETERRQ(PETSC_ERR_SUP,0,"Not for this matrix type"); 1969 ierr = (*mat->ops.norm)(mat,type,norm);CHKERRQ(ierr); 1970 PetscFunctionReturn(0); 1971 } 1972 1973 /* 1974 This variable is used to prevent counting of MatAssemblyBegin() that 1975 are called from within a MatAssemblyEnd(). 1976 */ 1977 static int MatAssemblyEnd_InUse = 0; 1978 #undef __FUNC__ 1979 #define __FUNC__ "MatAssemblyBegin" 1980 /*@ 1981 MatAssemblyBegin - Begins assembling the matrix. This routine should 1982 be called after completing all calls to MatSetValues(). 1983 1984 Input Parameters: 1985 . mat - the matrix 1986 . type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 1987 1988 Notes: 1989 MatSetValues() generally caches the values. The matrix is ready to 1990 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 1991 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 1992 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 1993 using the matrix. 1994 1995 .keywords: matrix, assembly, assemble, begin 1996 1997 .seealso: MatAssemblyEnd(), MatSetValues() 1998 @*/ 1999 int MatAssemblyBegin(Mat mat,MatAssemblyType type) 2000 { 2001 int ierr; 2002 2003 PetscFunctionBegin; 2004 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2005 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix.\n did you forget to call MatSetUnfactored()?"); 2006 if (mat->assembled) { 2007 mat->was_assembled = PETSC_TRUE; 2008 mat->assembled = PETSC_FALSE; 2009 } 2010 if (!MatAssemblyEnd_InUse) { 2011 PLogEventBegin(MAT_AssemblyBegin,mat,0,0,0); 2012 if (mat->ops.assemblybegin){ierr = (*mat->ops.assemblybegin)(mat,type);CHKERRQ(ierr);} 2013 PLogEventEnd(MAT_AssemblyBegin,mat,0,0,0); 2014 } else { 2015 if (mat->ops.assemblybegin){ierr = (*mat->ops.assemblybegin)(mat,type);CHKERRQ(ierr);} 2016 } 2017 PetscFunctionReturn(0); 2018 } 2019 2020 2021 #undef __FUNC__ 2022 #define __FUNC__ "MatView_Private" 2023 /* 2024 Processes command line options to determine if/how a matrix 2025 is to be viewed. Called by MatAssemblyEnd() and MatLoad(). 2026 */ 2027 int MatView_Private(Mat mat) 2028 { 2029 int ierr,flg; 2030 2031 PetscFunctionBegin; 2032 ierr = OptionsHasName(PETSC_NULL,"-mat_view_info",&flg); CHKERRQ(ierr); 2033 if (flg) { 2034 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO,0);CHKERRQ(ierr); 2035 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2036 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2037 } 2038 ierr = OptionsHasName(PETSC_NULL,"-mat_view_info_detailed",&flg);CHKERRQ(ierr); 2039 if (flg) { 2040 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO_LONG,0);CHKERRQ(ierr); 2041 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2042 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2043 } 2044 ierr = OptionsHasName(PETSC_NULL,"-mat_view",&flg); CHKERRQ(ierr); 2045 if (flg) { 2046 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2047 } 2048 ierr = OptionsHasName(PETSC_NULL,"-mat_view_matlab",&flg); CHKERRQ(ierr); 2049 if (flg) { 2050 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_MATLAB,"M");CHKERRQ(ierr); 2051 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2052 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2053 } 2054 ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg); CHKERRQ(ierr); 2055 if (flg) { 2056 ierr = OptionsHasName(0,"-mat_view_contour",&flg); CHKERRQ(ierr); 2057 if (flg) { 2058 ViewerPushFormat(VIEWER_DRAWX_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr); 2059 } 2060 ierr = MatView(mat,VIEWER_DRAWX_(mat->comm)); CHKERRQ(ierr); 2061 ierr = ViewerFlush(VIEWER_DRAWX_(mat->comm)); CHKERRQ(ierr); 2062 if (flg) { 2063 ViewerPopFormat(VIEWER_DRAWX_(mat->comm));CHKERRQ(ierr); 2064 } 2065 } 2066 PetscFunctionReturn(0); 2067 } 2068 2069 #undef __FUNC__ 2070 #define __FUNC__ "MatAssemblyEnd" 2071 /*@ 2072 MatAssemblyEnd - Completes assembling the matrix. This routine should 2073 be called after MatAssemblyBegin(). 2074 2075 Input Parameters: 2076 . mat - the matrix 2077 . type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 2078 2079 Options Database Keys: 2080 $ -mat_view_info : Prints info on matrix at 2081 $ conclusion of MatEndAssembly() 2082 $ -mat_view_info_detailed: Prints more detailed info. 2083 $ -mat_view : Prints matrix in ASCII format. 2084 $ -mat_view_matlab : Prints matrix in Matlab format. 2085 $ -mat_view_draw : Draws nonzero structure of matrix, 2086 $ using MatView() and DrawOpenX(). 2087 $ -display <name> : Set display name (default is host) 2088 $ -draw_pause <sec> : Set number of seconds to pause after display 2089 2090 Notes: 2091 MatSetValues() generally caches the values. The matrix is ready to 2092 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 2093 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 2094 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 2095 using the matrix. 2096 2097 .keywords: matrix, assembly, assemble, end 2098 2099 .seealso: MatAssemblyBegin(), MatSetValues(), DrawOpenX(), MatView() 2100 @*/ 2101 int MatAssemblyEnd(Mat mat,MatAssemblyType type) 2102 { 2103 int ierr; 2104 static int inassm = 0; 2105 2106 PetscFunctionBegin; 2107 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2108 2109 inassm++; 2110 MatAssemblyEnd_InUse++; 2111 PLogEventBegin(MAT_AssemblyEnd,mat,0,0,0); 2112 if (mat->ops.assemblyend) { 2113 ierr = (*mat->ops.assemblyend)(mat,type); CHKERRQ(ierr); 2114 } 2115 2116 /* Flush assembly is not a true assembly */ 2117 if (type != MAT_FLUSH_ASSEMBLY) { 2118 mat->assembled = PETSC_TRUE; mat->num_ass++; 2119 } 2120 mat->insertmode = NOT_SET_VALUES; 2121 PLogEventEnd(MAT_AssemblyEnd,mat,0,0,0); 2122 MatAssemblyEnd_InUse--; 2123 2124 if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) { 2125 ierr = MatView_Private(mat); CHKERRQ(ierr); 2126 } 2127 inassm--; 2128 PetscFunctionReturn(0); 2129 } 2130 2131 2132 #undef __FUNC__ 2133 #define __FUNC__ "MatCompress" 2134 /*@ 2135 MatCompress - Tries to store the matrix in as little space as 2136 possible. May fail if memory is already fully used, since it 2137 tries to allocate new space. 2138 2139 Input Parameters: 2140 . mat - the matrix 2141 2142 .keywords: matrix, compress 2143 @*/ 2144 int MatCompress(Mat mat) 2145 { 2146 int ierr; 2147 2148 PetscFunctionBegin; 2149 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2150 if (mat->ops.compress) {ierr = (*mat->ops.compress)(mat);CHKERRQ(ierr);} 2151 PetscFunctionReturn(0); 2152 } 2153 2154 #undef __FUNC__ 2155 #define __FUNC__ "MatSetOption" 2156 /*@ 2157 MatSetOption - Sets a parameter option for a matrix. Some options 2158 may be specific to certain storage formats. Some options 2159 determine how values will be inserted (or added). Sorted, 2160 row-oriented input will generally assemble the fastest. The default 2161 is row-oriented, nonsorted input. 2162 2163 Input Parameters: 2164 . mat - the matrix 2165 . option - the option, one of those listed below (and possibly others), 2166 e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERROR 2167 2168 Options Describing Matrix Structure: 2169 $ MAT_SYMMETRIC - symmetric in terms of both structure and value 2170 $ MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure 2171 2172 Options For Use with MatSetValues(): 2173 Insert a logically dense subblock, which can be 2174 2175 $ MAT_ROW_ORIENTED - row-oriented 2176 $ MAT_COLUMN_ORIENTED - column-oriented 2177 $ MAT_ROWS_SORTED - sorted by row 2178 $ MAT_ROWS_UNSORTED - not sorted by row 2179 $ MAT_COLUMNS_SORTED - sorted by column 2180 $ MAT_COLUMNS_UNSORTED - not sorted by column 2181 2182 Not these options reflect the data you pass in with MatSetValues(); it has 2183 nothing to do with how the data is stored internally in the matrix 2184 data structure. 2185 2186 When (re)assembling a matrix, we can restrict the input for 2187 efficiency/debugging purposes. 2188 2189 $ MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be 2190 allowed if they generate a new nonzero 2191 $ MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed 2192 $ MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if 2193 they generate a nonzero in a new diagonal (for block diagonal format only) 2194 $ MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only) 2195 $ MAT_IGNORE_OFF_PROC_ENTRIES - drop off-processor entries 2196 $ MAT_NEW_NONZERO_LOCATION_ERROR - generate error for new matrix entry 2197 2198 Notes: 2199 Some options are relevant only for particular matrix types and 2200 are thus ignored by others. Other options are not supported by 2201 certain matrix types and will generate an error message if set. 2202 2203 If using a Fortran 77 module to compute a matrix, one may need to 2204 use the column-oriented option (or convert to the row-oriented 2205 format). 2206 2207 MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion 2208 that would generate a new entry in the nonzero structure is instead 2209 ignored. Thus, if memory has not alredy been allocated for this particular 2210 data, then the insertion is ignored. For dense matrices, in which 2211 the entire array is allocated, no entries are ever ignored. 2212 2213 MAT_NEW_NONZERO_LOCATION_ERROR indicates that any add or insertion 2214 that would generate a new entry in the nonzero structure instead produces 2215 an error. (Currently supported for AIJ and BAIJ formats only.) 2216 This is a useful flag when using SAME_NONZERO_PATTERN in calling 2217 SLESSetOperators() to ensure that the nonzero pattern truely does 2218 remain unchanged. 2219 2220 MAT_NEW_NONZERO_ALLOCATION_ERROR indicates that any add or insertion 2221 that would generate a new entry that has not been preallocated will 2222 instead produce an error. (Currently supported for AIJ and BAIJ formats 2223 only.) This is a useful flag when debugging matrix memory preallocation. 2224 2225 MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for 2226 other processors should be dropped, rather than stashed. 2227 This is useful if you know that the "owning" processor is also 2228 always generating the correct matrix entries, so that PETSc need 2229 not transfer duplicate entries generated on another processor. 2230 2231 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero 2232 @*/ 2233 int MatSetOption(Mat mat,MatOption op) 2234 { 2235 int ierr; 2236 2237 PetscFunctionBegin; 2238 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2239 if (mat->ops.setoption) {ierr = (*mat->ops.setoption)(mat,op);CHKERRQ(ierr);} 2240 PetscFunctionReturn(0); 2241 } 2242 2243 #undef __FUNC__ 2244 #define __FUNC__ "MatZeroEntries" 2245 /*@ 2246 MatZeroEntries - Zeros all entries of a matrix. For sparse matrices 2247 this routine retains the old nonzero structure. 2248 2249 Input Parameters: 2250 . mat - the matrix 2251 2252 .keywords: matrix, zero, entries 2253 2254 .seealso: MatZeroRows() 2255 @*/ 2256 int MatZeroEntries(Mat mat) 2257 { 2258 int ierr; 2259 2260 PetscFunctionBegin; 2261 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2262 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2263 if (!mat->ops.zeroentries) SETERRQ(PETSC_ERR_SUP,0,""); 2264 2265 PLogEventBegin(MAT_ZeroEntries,mat,0,0,0); 2266 ierr = (*mat->ops.zeroentries)(mat); CHKERRQ(ierr); 2267 PLogEventEnd(MAT_ZeroEntries,mat,0,0,0); 2268 PetscFunctionReturn(0); 2269 } 2270 2271 #undef __FUNC__ 2272 #define __FUNC__ "MatZeroRows" 2273 /*@ 2274 MatZeroRows - Zeros all entries (except possibly the main diagonal) 2275 of a set of rows of a matrix. 2276 2277 Input Parameters: 2278 . mat - the matrix 2279 . is - index set of rows to remove 2280 . diag - pointer to value put in all diagonals of eliminated rows. 2281 Note that diag is not a pointer to an array, but merely a 2282 pointer to a single value. 2283 2284 Notes: 2285 For the AIJ matrix formats this removes the old nonzero structure, 2286 but does not release memory. For the dense and block diagonal 2287 formats this does not alter the nonzero structure. 2288 2289 The user can set a value in the diagonal entry (or for the AIJ and 2290 row formats can optionally remove the main diagonal entry from the 2291 nonzero structure as well, by passing a null pointer as the final 2292 argument). 2293 2294 For the parallel case, all processes that share the matrix (i.e., 2295 those in the communicator used for matrix creation) MUST call this 2296 routine, regardless of whether any rows being zeroed are owned by 2297 them. 2298 2299 .keywords: matrix, zero, rows, boundary conditions 2300 2301 .seealso: MatZeroEntries(), 2302 @*/ 2303 int MatZeroRows(Mat mat,IS is, Scalar *diag) 2304 { 2305 int ierr; 2306 2307 PetscFunctionBegin; 2308 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2309 PetscValidHeaderSpecific(is,IS_COOKIE); 2310 if (diag) PetscValidScalarPointer(diag); 2311 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2312 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2313 if (!mat->ops.zerorows) SETERRQ(PETSC_ERR_SUP,0,""); 2314 2315 ierr = (*mat->ops.zerorows)(mat,is,diag); CHKERRQ(ierr); 2316 ierr = MatView_Private(mat); CHKERRQ(ierr); 2317 PetscFunctionReturn(0); 2318 } 2319 2320 #undef __FUNC__ 2321 #define __FUNC__ "MatZeroRowsLocal" 2322 /*@ 2323 MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal) 2324 of a set of rows of a matrix; using local numbering of rows. 2325 2326 Input Parameters: 2327 . mat - the matrix 2328 . is - index set of rows to remove 2329 . diag - pointer to value put in all diagonals of eliminated rows. 2330 Note that diag is not a pointer to an array, but merely a 2331 pointer to a single value. 2332 2333 Notes: 2334 For the AIJ matrix formats this removes the old nonzero structure, 2335 but does not release memory. For the dense and block diagonal 2336 formats this does not alter the nonzero structure. 2337 2338 The user can set a value in the diagonal entry (or for the AIJ and 2339 row formats can optionally remove the main diagonal entry from the 2340 nonzero structure as well, by passing a null pointer as the final 2341 argument). 2342 2343 .keywords: matrix, zero, rows, boundary conditions 2344 2345 .seealso: MatZeroEntries(), 2346 @*/ 2347 int MatZeroRowsLocal(Mat mat,IS is, Scalar *diag) 2348 { 2349 int ierr; 2350 IS newis; 2351 2352 PetscFunctionBegin; 2353 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2354 PetscValidHeaderSpecific(is,IS_COOKIE); 2355 if (diag) PetscValidScalarPointer(diag); 2356 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2357 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2358 if (!mat->ops.zerorows) SETERRQ(PETSC_ERR_SUP,0,""); 2359 2360 ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis); CHKERRQ(ierr); 2361 ierr = (*mat->ops.zerorows)(mat,newis,diag); CHKERRQ(ierr); 2362 ierr = ISDestroy(newis); 2363 PetscFunctionReturn(0); 2364 } 2365 2366 #undef __FUNC__ 2367 #define __FUNC__ "MatGetSize" 2368 /*@ 2369 MatGetSize - Returns the numbers of rows and columns in a matrix. 2370 2371 Input Parameter: 2372 . mat - the matrix 2373 2374 Output Parameters: 2375 . m - the number of global rows 2376 . n - the number of global columns 2377 2378 .keywords: matrix, dimension, size, rows, columns, global, get 2379 2380 .seealso: MatGetLocalSize() 2381 @*/ 2382 int MatGetSize(Mat mat,int *m,int* n) 2383 { 2384 int ierr; 2385 2386 PetscFunctionBegin; 2387 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2388 ierr = (*mat->ops.getsize)(mat,m,n);CHKERRQ(ierr); 2389 PetscFunctionReturn(0); 2390 } 2391 2392 #undef __FUNC__ 2393 #define __FUNC__ "MatGetLocalSize" 2394 /*@ 2395 MatGetLocalSize - Returns the number of rows and columns in a matrix 2396 stored locally. This information may be implementation dependent, so 2397 use with care. 2398 2399 Input Parameters: 2400 . mat - the matrix 2401 2402 Output Parameters: 2403 . m - the number of local rows 2404 . n - the number of local columns 2405 2406 .keywords: matrix, dimension, size, local, rows, columns, get 2407 2408 .seealso: MatGetSize() 2409 @*/ 2410 int MatGetLocalSize(Mat mat,int *m,int* n) 2411 { 2412 int ierr; 2413 2414 PetscFunctionBegin; 2415 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2416 ierr = (*mat->ops.getlocalsize)(mat,m,n);CHKERRQ(ierr); 2417 PetscFunctionReturn(0); 2418 } 2419 2420 #undef __FUNC__ 2421 #define __FUNC__ "MatGetOwnershipRange" 2422 /*@ 2423 MatGetOwnershipRange - Returns the range of matrix rows owned by 2424 this processor, assuming that the matrix is laid out with the first 2425 n1 rows on the first processor, the next n2 rows on the second, etc. 2426 For certain parallel layouts this range may not be well defined. 2427 2428 Input Parameters: 2429 . mat - the matrix 2430 2431 Output Parameters: 2432 . m - the global index of the first local row 2433 . n - one more than the global index of the last local row 2434 2435 .keywords: matrix, get, range, ownership 2436 @*/ 2437 int MatGetOwnershipRange(Mat mat,int *m,int* n) 2438 { 2439 int ierr; 2440 2441 PetscFunctionBegin; 2442 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2443 PetscValidIntPointer(m); 2444 PetscValidIntPointer(n); 2445 if (!mat->ops.getownershiprange) SETERRQ(PETSC_ERR_SUP,0,""); 2446 ierr = (*mat->ops.getownershiprange)(mat,m,n);CHKERRQ(ierr); 2447 PetscFunctionReturn(0); 2448 } 2449 2450 #undef __FUNC__ 2451 #define __FUNC__ "MatILUFactorSymbolic" 2452 /*@ 2453 MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix. 2454 Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric() 2455 to complete the factorization. 2456 2457 Input Parameters: 2458 . mat - the matrix 2459 . row - row permutation 2460 . column - column permutation 2461 . fill - number of levels of fill 2462 . f - expected fill as ratio of the original number of nonzeros, 2463 for example 3.0; choosing this parameter well can result in 2464 more efficient use of time and space. Run your code with -log_info 2465 to determine an optimal value to use. 2466 2467 Output Parameters: 2468 . fact - new matrix that has been symbolically factored 2469 2470 Notes: 2471 See the file $(PETSC_DIR)/Performace for additional information about 2472 choosing the fill factor for better efficiency. 2473 2474 .keywords: matrix, factor, incomplete, ILU, symbolic, fill 2475 2476 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric() 2477 @*/ 2478 int MatILUFactorSymbolic(Mat mat,IS row,IS col,double f,int fill,Mat *fact) 2479 { 2480 int ierr; 2481 2482 PetscFunctionBegin; 2483 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2484 PetscValidPointer(fact); 2485 if (fill < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative"); 2486 if (!mat->ops.ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support parallel ILU"); 2487 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2488 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2489 2490 PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0); 2491 ierr = (*mat->ops.ilufactorsymbolic)(mat,row,col,f,fill,fact); CHKERRQ(ierr); 2492 PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0); 2493 PetscFunctionReturn(0); 2494 } 2495 2496 #undef __FUNC__ 2497 #define __FUNC__ "MatIncompleteCholeskyFactorSymbolic" 2498 /*@ 2499 MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete 2500 Cholesky factorization for a symmetric matrix. Use 2501 MatCholeskyFactorNumeric() to complete the factorization. 2502 2503 Input Parameters: 2504 . mat - the matrix 2505 . perm - row and column permutation 2506 . fill - levels of fill 2507 . f - expected fill as ratio of original fill 2508 2509 Output Parameter: 2510 . fact - the factored matrix 2511 2512 Note: Currently only no-fill factorization is supported. 2513 2514 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill 2515 2516 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor() 2517 @*/ 2518 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,double f,int fill,Mat *fact) 2519 { 2520 int ierr; 2521 2522 PetscFunctionBegin; 2523 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2524 PetscValidPointer(fact); 2525 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2526 if (fill < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative"); 2527 if (!mat->ops.incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel"); 2528 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2529 2530 PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0); 2531 ierr = (*mat->ops.incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr); 2532 PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0); 2533 PetscFunctionReturn(0); 2534 } 2535 2536 #undef __FUNC__ 2537 #define __FUNC__ "MatGetArray" 2538 /*@C 2539 MatGetArray - Returns a pointer to the element values in the matrix. 2540 This routine is implementation dependent, and may not even work for 2541 certain matrix types. You MUST call MatRestoreArray() when you no 2542 longer need to access the array. 2543 2544 Input Parameter: 2545 . mat - the matrix 2546 2547 Output Parameter: 2548 . v - the location of the values 2549 2550 Fortran Note: 2551 The Fortran interface is slightly different from that given below. 2552 See the Fortran chapter of the users manual and 2553 petsc/src/mat/examples for details. 2554 2555 .keywords: matrix, array, elements, values 2556 2557 .seealso: MatRestoreArray() 2558 @*/ 2559 int MatGetArray(Mat mat,Scalar **v) 2560 { 2561 int ierr; 2562 2563 PetscFunctionBegin; 2564 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2565 PetscValidPointer(v); 2566 if (!mat->ops.getarray) SETERRQ(PETSC_ERR_SUP,0,""); 2567 ierr = (*mat->ops.getarray)(mat,v);CHKERRQ(ierr); 2568 PetscFunctionReturn(0); 2569 } 2570 2571 #undef __FUNC__ 2572 #define __FUNC__ "MatRestoreArray" 2573 /*@C 2574 MatRestoreArray - Restores the matrix after MatGetArray has been called. 2575 2576 Input Parameter: 2577 . mat - the matrix 2578 . v - the location of the values 2579 2580 Fortran Note: 2581 The Fortran interface is slightly different from that given below. 2582 See the users manual and petsc/src/mat/examples for details. 2583 2584 .keywords: matrix, array, elements, values, restore 2585 2586 .seealso: MatGetArray() 2587 @*/ 2588 int MatRestoreArray(Mat mat,Scalar **v) 2589 { 2590 int ierr; 2591 2592 PetscFunctionBegin; 2593 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2594 PetscValidPointer(v); 2595 if (!mat->ops.restorearray) SETERRQ(PETSC_ERR_SUP,0,""); 2596 ierr = (*mat->ops.restorearray)(mat,v);CHKERRQ(ierr); 2597 PetscFunctionReturn(0); 2598 } 2599 2600 #undef __FUNC__ 2601 #define __FUNC__ "MatGetSubMatrices" 2602 /*@C 2603 MatGetSubMatrices - Extracts several submatrices from a matrix. If submat 2604 points to an array of valid matrices, they may be reused to store the new 2605 submatrices. 2606 2607 Input Parameters: 2608 . mat - the matrix 2609 . n - the number of submatrixes to be extracted 2610 . irow, icol - index sets of rows and columns to extract 2611 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 2612 2613 Output Parameter: 2614 . submat - the array of submatrices 2615 2616 Notes: 2617 MatGetSubMatrices() can extract only sequential submatrices 2618 (from both sequential and parallel matrices). Use MatGetSubMatrix() 2619 to extract a parallel submatrix. 2620 2621 When extracting submatrices from a parallel matrix, each processor can 2622 form a different submatrix by setting the rows and columns of its 2623 individual index sets according to the local submatrix desired. 2624 2625 When finished using the submatrices, the user should destroy 2626 them with MatDestroySubMatrices(). 2627 2628 .keywords: matrix, get, submatrix, submatrices 2629 2630 .seealso: MatDestroyMatrices(), MatGetSubMatrix() 2631 @*/ 2632 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **submat) 2633 { 2634 int ierr; 2635 2636 PetscFunctionBegin; 2637 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2638 if (!mat->ops.getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,""); 2639 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2640 2641 PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0); 2642 ierr = (*mat->ops.getsubmatrices)(mat,n,irow,icol,scall,submat); CHKERRQ(ierr); 2643 PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0); 2644 2645 PetscFunctionReturn(0); 2646 } 2647 2648 #undef __FUNC__ 2649 #define __FUNC__ "MatDestroyMatrices" 2650 /*@C 2651 MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices(). 2652 2653 Input Parameters: 2654 . n - the number of local matrices 2655 . mat - the matrices 2656 2657 .keywords: matrix, destroy, submatrix, submatrices 2658 2659 .seealso: MatGetSubMatrices() 2660 @*/ 2661 int MatDestroyMatrices(int n,Mat **mat) 2662 { 2663 int ierr,i; 2664 2665 PetscFunctionBegin; 2666 if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices"); 2667 PetscValidPointer(mat); 2668 for ( i=0; i<n; i++ ) { 2669 ierr = MatDestroy((*mat)[i]); CHKERRQ(ierr); 2670 } 2671 if (n) PetscFree(*mat); 2672 PetscFunctionReturn(0); 2673 } 2674 2675 #undef __FUNC__ 2676 #define __FUNC__ "MatIncreaseOverlap" 2677 /*@ 2678 MatIncreaseOverlap - Given a set of submatrices indicated by index sets, 2679 replaces the index sets by larger ones that represent submatrices with 2680 additional overlap. 2681 2682 Input Parameters: 2683 . mat - the matrix 2684 . n - the number of index sets 2685 . is - the array of pointers to index sets 2686 . ov - the additional overlap requested 2687 2688 .keywords: matrix, overlap, Schwarz 2689 2690 .seealso: MatGetSubMatrices() 2691 @*/ 2692 int MatIncreaseOverlap(Mat mat,int n, IS *is,int ov) 2693 { 2694 int ierr; 2695 2696 PetscFunctionBegin; 2697 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2698 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2699 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2700 2701 if (ov == 0) PetscFunctionReturn(0); 2702 if (!mat->ops.increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,""); 2703 PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0); 2704 ierr = (*mat->ops.increaseoverlap)(mat,n,is,ov); CHKERRQ(ierr); 2705 PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0); 2706 PetscFunctionReturn(0); 2707 } 2708 2709 #undef __FUNC__ 2710 #define __FUNC__ "MatPrintHelp" 2711 /*@ 2712 MatPrintHelp - Prints all the options for the matrix. 2713 2714 Input Parameter: 2715 . mat - the matrix 2716 2717 Options Database Keys: 2718 $ -help, -h 2719 2720 .keywords: mat, help 2721 2722 .seealso: MatCreate(), MatCreateXXX() 2723 @*/ 2724 int MatPrintHelp(Mat mat) 2725 { 2726 static int called = 0; 2727 int ierr; 2728 MPI_Comm comm; 2729 2730 PetscFunctionBegin; 2731 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2732 2733 comm = mat->comm; 2734 if (!called) { 2735 (*PetscHelpPrintf)(comm,"General matrix options:\n"); 2736 (*PetscHelpPrintf)(comm," -mat_view_info: view basic matrix info during MatAssemblyEnd()\n"); 2737 (*PetscHelpPrintf)(comm," -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n"); 2738 (*PetscHelpPrintf)(comm," -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n"); 2739 (*PetscHelpPrintf)(comm," -draw_pause <sec>: set seconds of display pause\n"); 2740 (*PetscHelpPrintf)(comm," -display <name>: set alternate display\n"); 2741 called = 1; 2742 } 2743 if (mat->ops.printhelp) { 2744 ierr = (*mat->ops.printhelp)(mat); CHKERRQ(ierr); 2745 } 2746 PetscFunctionReturn(0); 2747 } 2748 2749 #undef __FUNC__ 2750 #define __FUNC__ "MatGetBlockSize" 2751 /*@ 2752 MatGetBlockSize - Returns the matrix block size; useful especially for the 2753 block row and block diagonal formats. 2754 2755 Input Parameter: 2756 . mat - the matrix 2757 2758 Output Parameter: 2759 . bs - block size 2760 2761 Notes: 2762 $ block diagonal formats: MATSEQBDIAG, MATMPIBDIAG 2763 $ block row formats: MATSEQBAIJ, MATMPIBAIJ 2764 2765 .keywords: matrix, get, block, size 2766 2767 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag() 2768 @*/ 2769 int MatGetBlockSize(Mat mat,int *bs) 2770 { 2771 int ierr; 2772 2773 PetscFunctionBegin; 2774 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2775 PetscValidIntPointer(bs); 2776 if (!mat->ops.getblocksize) SETERRQ(PETSC_ERR_SUP,0,""); 2777 ierr = (*mat->ops.getblocksize)(mat,bs);CHKERRQ(ierr); 2778 PetscFunctionReturn(0); 2779 } 2780 2781 #undef __FUNC__ 2782 #define __FUNC__ "MatGetRowIJ" 2783 /*@C 2784 MatGetRowIJ - Returns the compress row storage i and j indices for sequential matrices. 2785 EXPERTS ONLY. 2786 2787 Input Parameters: 2788 . mat - the matrix 2789 . shift - 0 or 1 indicating we want the indices starting at 0 or 1 2790 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2791 symmetrized 2792 2793 Output Parameters: 2794 . n - number of rows and columns in the (possibly compressed) matrix 2795 . ia - the row indices 2796 . ja - the column indices 2797 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2798 @*/ 2799 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2800 { 2801 int ierr; 2802 2803 PetscFunctionBegin; 2804 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2805 if (ia) PetscValidIntPointer(ia); 2806 if (ja) PetscValidIntPointer(ja); 2807 PetscValidIntPointer(done); 2808 if (!mat->ops.getrowij) *done = PETSC_FALSE; 2809 else { 2810 *done = PETSC_TRUE; 2811 ierr = (*mat->ops.getrowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2812 } 2813 PetscFunctionReturn(0); 2814 } 2815 2816 #undef __FUNC__ 2817 #define __FUNC__ "MatGetColumnIJ" 2818 /*@C 2819 MatGetColumnIJ - Returns the compress Column storage i and j indices for sequential matrices. 2820 EXPERTS ONLY. 2821 2822 Input Parameters: 2823 . mat - the matrix 2824 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2825 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2826 symmetrized 2827 2828 Output Parameters: 2829 . n - number of Columns and columns in the (possibly compressed) matrix 2830 . ia - the Column indices 2831 . ja - the column indices 2832 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2833 @*/ 2834 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2835 { 2836 int ierr; 2837 2838 PetscFunctionBegin; 2839 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2840 if (ia) PetscValidIntPointer(ia); 2841 if (ja) PetscValidIntPointer(ja); 2842 PetscValidIntPointer(done); 2843 2844 if (!mat->ops.getcolumnij) *done = PETSC_FALSE; 2845 else { 2846 *done = PETSC_TRUE; 2847 ierr = (*mat->ops.getcolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2848 } 2849 PetscFunctionReturn(0); 2850 } 2851 2852 #undef __FUNC__ 2853 #define __FUNC__ "MatRestoreRowIJ" 2854 /*@C 2855 MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with 2856 MatGetRowIJ(). EXPERTS ONLY. 2857 2858 Input Parameters: 2859 . mat - the matrix 2860 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2861 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2862 symmetrized 2863 2864 Output Parameters: 2865 . n - size of (possibly compressed) matrix 2866 . ia - the row indices 2867 . ja - the column indices 2868 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2869 @*/ 2870 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2871 { 2872 int ierr; 2873 2874 PetscFunctionBegin; 2875 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2876 if (ia) PetscValidIntPointer(ia); 2877 if (ja) PetscValidIntPointer(ja); 2878 PetscValidIntPointer(done); 2879 2880 if (!mat->ops.restorerowij) *done = PETSC_FALSE; 2881 else { 2882 *done = PETSC_TRUE; 2883 ierr = (*mat->ops.restorerowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2884 } 2885 PetscFunctionReturn(0); 2886 } 2887 2888 #undef __FUNC__ 2889 #define __FUNC__ "MatRestoreColumnIJ" 2890 /*@C 2891 MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with 2892 MatGetColumnIJ(). EXPERTS ONLY. 2893 2894 Input Parameters: 2895 . mat - the matrix 2896 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2897 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2898 symmetrized 2899 2900 Output Parameters: 2901 . n - size of (possibly compressed) matrix 2902 . ia - the Column indices 2903 . ja - the column indices 2904 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2905 @*/ 2906 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2907 { 2908 int ierr; 2909 2910 PetscFunctionBegin; 2911 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2912 if (ia) PetscValidIntPointer(ia); 2913 if (ja) PetscValidIntPointer(ja); 2914 PetscValidIntPointer(done); 2915 2916 if (!mat->ops.restorecolumnij) *done = PETSC_FALSE; 2917 else { 2918 *done = PETSC_TRUE; 2919 ierr = (*mat->ops.restorecolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2920 } 2921 PetscFunctionReturn(0); 2922 } 2923 2924 #undef __FUNC__ 2925 #define __FUNC__ "MatColoringPatch" 2926 /*@C 2927 MatColoringPatch - EXPERTS ONLY, used inside matrix coloring routines that 2928 use matGetRowIJ() and/or MatGetColumnIJ(). 2929 2930 Input Parameters: 2931 . mat - the matrix 2932 . n - number of colors 2933 . colorarray - array indicating color for each column 2934 2935 Output Parameters: 2936 . iscoloring - coloring generated using colorarray information 2937 2938 @*/ 2939 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring) 2940 { 2941 int ierr; 2942 2943 PetscFunctionBegin; 2944 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2945 PetscValidIntPointer(colorarray); 2946 2947 if (!mat->ops.coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");} 2948 else { 2949 ierr = (*mat->ops.coloringpatch)(mat,n,colorarray,iscoloring); CHKERRQ(ierr); 2950 } 2951 PetscFunctionReturn(0); 2952 } 2953 2954 2955 #undef __FUNC__ 2956 #define __FUNC__ "MatSetUnfactored" 2957 /*@ 2958 MatSetUnfactored - Resets a factored matrix to be treated as unfactored. 2959 2960 Input Paramter: 2961 . mat - the factored matrix to be reset 2962 2963 Notes: 2964 This routine should be used only with factored matrices formed by in-place 2965 factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE 2966 format). This option can save memory, for example, when solving nonlinear 2967 systems with a matrix-free Newton-Krylov method and a matrix-based, in-place 2968 ILU(0) preconditioner. 2969 2970 Note that one can specify in-place ILU(0) factorization by calling 2971 $ PCType(pc,PCILU); 2972 $ PCILUSeUseInPlace(pc); 2973 or by using the options -pc_type ilu -pc_ilu_in_place 2974 2975 In-place factorization ILU(0) can also be used as a local 2976 solver for the blocks within the block Jacobi or additive Schwarz 2977 methods (runtime option: -sub_pc_ilu_in_place). See the discussion 2978 of these preconditioners in the users manual for details on setting 2979 local solver options. 2980 2981 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace() 2982 2983 .keywords: matrix-free, in-place ILU, in-place LU 2984 @*/ 2985 int MatSetUnfactored(Mat mat) 2986 { 2987 int ierr; 2988 2989 PetscFunctionBegin; 2990 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2991 mat->factor = 0; 2992 if (!mat->ops.setunfactored) PetscFunctionReturn(0); 2993 ierr = (*mat->ops.setunfactored)(mat); CHKERRQ(ierr); 2994 PetscFunctionReturn(0); 2995 } 2996 2997 #undef __FUNC__ 2998 #define __FUNC__ "MatGetType" 2999 /*@C 3000 MatGetType - Gets the matrix type and name (as a string) from the matrix. 3001 3002 Input Parameter: 3003 . mat - the matrix 3004 3005 Output Parameter: 3006 . type - the matrix type (or use PETSC_NULL) 3007 . name - name of matrix type (or use PETSC_NULL) 3008 3009 .keywords: matrix, get, type, name 3010 @*/ 3011 int MatGetType(Mat mat,MatType *type,char **name) 3012 { 3013 int itype = (int)mat->type; 3014 char *matname[10]; 3015 3016 PetscFunctionBegin; 3017 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3018 3019 if (type) *type = (MatType) mat->type; 3020 if (name) { 3021 /* Note: Be sure that this list corresponds to the enum in mat.h */ 3022 matname[0] = "MATSEQDENSE"; 3023 matname[1] = "MATSEQAIJ"; 3024 matname[2] = "MATMPIAIJ"; 3025 matname[3] = "MATSHELL"; 3026 matname[4] = "MATMPIROWBS"; 3027 matname[5] = "MATSEQBDIAG"; 3028 matname[6] = "MATMPIBDIAG"; 3029 matname[7] = "MATMPIDENSE"; 3030 matname[8] = "MATSEQBAIJ"; 3031 matname[9] = "MATMPIBAIJ"; 3032 3033 if (itype < 0 || itype > 9) *name = "Unknown matrix type"; 3034 else *name = matname[itype]; 3035 } 3036 PetscFunctionReturn(0); 3037 } 3038 3039 /*MC 3040 MatGetArrayF90 - Accesses a matrix array from Fortran90. 3041 3042 Input Parameter: 3043 . x - matrix 3044 3045 Output Parameter: 3046 . xx_v - the Fortran90 pointer to the array 3047 . ierr - error code 3048 3049 Synopsis: 3050 MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 3051 3052 Usage: 3053 $ Scalar, pointer xx_v(:) 3054 $ .... 3055 $ call MatGetArrayF90(x,xx_v,ierr) 3056 $ a = xx_v(3) 3057 $ call MatRestoreArrayF90(x,xx_v,ierr) 3058 3059 Notes: 3060 Currently only supported using the NAG F90 compiler. 3061 3062 .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray() 3063 3064 .keywords: matrix, array, f90 3065 M*/ 3066 3067 /*MC 3068 MatRestoreArrayF90 - Restores a matrix array that has been 3069 accessed with MatGetArrayF90(). 3070 3071 Input Parameters: 3072 . x - matrix 3073 . xx_v - the Fortran90 pointer to the array 3074 3075 Output Parameter: 3076 . ierr - error code 3077 3078 Synopsis: 3079 MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 3080 3081 Example of Usage: 3082 $ Scalar, pointer xx_v(:) 3083 $ .... 3084 $ call MatGetArrayF90(x,xx_v,ierr) 3085 $ a = xx_v(3) 3086 $ call MatRestoreArrayF90(x,xx_v,ierr) 3087 3088 Notes: 3089 Currently only supported using the NAG F90 compiler. 3090 3091 .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray() 3092 3093 .keywords: matrix, array, f90 3094 M*/ 3095 3096 3097 #undef __FUNC__ 3098 #define __FUNC__ "MatGetSubMatrix" 3099 /*@ 3100 MatGetSubMatrix - Gets a single submatrix on the same number of processors 3101 as the original matrix. 3102 3103 Input Parameters: 3104 . mat - the original matrix 3105 . isrow - rows this processor should obtain 3106 . iscol - columns for all processors you wish kept 3107 . csize - number of columns "local" to this processor (does nothing for sequential 3108 matrices). This should match the result from VecGetLocalSize() if you 3109 plan to use the matrix in a A*x 3110 . cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 3111 3112 Output Parameters: 3113 . newmat - the new submatrix, of the same type as the old 3114 3115 .seealso: MatGetSubMatrices() 3116 3117 @*/ 3118 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatGetSubMatrixCall cll,Mat *newmat) 3119 { 3120 int ierr, size; 3121 Mat *local; 3122 3123 PetscFunctionBegin; 3124 MPI_Comm_size(mat->comm,&size); 3125 3126 /* if original matrix is on just one processor then use submatrix generated */ 3127 if (size == 1 && cll == MAT_REUSE_MATRIX) { 3128 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr); 3129 PetscFunctionReturn(0); 3130 } else if (size == 1) { 3131 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr); 3132 *newmat = *local; 3133 PetscFree(local); 3134 PetscFunctionReturn(0); 3135 } 3136 3137 if (!mat->ops.getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented"); 3138 ierr = (*mat->ops.getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr); 3139 PetscFunctionReturn(0); 3140 } 3141 3142 3143 3144