1 /*$Id: matrix.c,v 1.414 2001/09/28 18:57:28 balay Exp $*/ 2 3 /* 4 This is where the abstract matrix operations are defined 5 */ 6 7 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 8 #include "src/vec/vecimpl.h" 9 10 /* Logging support */ 11 int MAT_COOKIE = 0; 12 int MATSNESMFCTX_COOKIE = 0; 13 int MAT_Mult = 0, MAT_MultMatrixFree = 0, MAT_Mults = 0, MAT_MultConstrained = 0, MAT_MultAdd = 0, MAT_MultTranspose = 0; 14 int MAT_MultTransposeConstrained = 0, MAT_MultTransposeAdd = 0, MAT_Solve = 0, MAT_Solves = 0, MAT_SolveAdd = 0, MAT_SolveTranspose = 0; 15 int MAT_SolveTransposeAdd = 0, MAT_Relax = 0, MAT_ForwardSolve = 0, MAT_BackwardSolve = 0, MAT_LUFactor = 0, MAT_LUFactorSymbolic = 0; 16 int MAT_LUFactorNumeric = 0, MAT_CholeskyFactor = 0, MAT_CholeskyFactorSymbolic = 0, MAT_CholeskyFactorNumeric = 0, MAT_ILUFactor = 0; 17 int MAT_ILUFactorSymbolic = 0, MAT_ICCFactorSymbolic = 0, MAT_Copy = 0, MAT_Convert = 0, MAT_Scale = 0, MAT_AssemblyBegin = 0; 18 int MAT_AssemblyEnd = 0, MAT_SetValues = 0, MAT_GetValues = 0, MAT_GetRow = 0, MAT_GetSubMatrices = 0, MAT_GetColoring = 0, MAT_GetOrdering = 0; 19 int MAT_IncreaseOverlap = 0, MAT_Partitioning = 0, MAT_ZeroEntries = 0, MAT_Load = 0, MAT_View = 0, MAT_AXPY = 0, MAT_FDColoringCreate = 0; 20 int MAT_FDColoringApply = 0,MAT_Transpose = 0,MAT_FDColoringFunction = 0; 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "MatGetRow" 24 /*@C 25 MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow() 26 for each row that you get to ensure that your application does 27 not bleed memory. 28 29 Not Collective 30 31 Input Parameters: 32 + mat - the matrix 33 - row - the row to get 34 35 Output Parameters: 36 + ncols - if not NULL, the number of nonzeros in the row 37 . cols - if not NULL, the column numbers 38 - vals - if not NULL, the values 39 40 Notes: 41 This routine is provided for people who need to have direct access 42 to the structure of a matrix. We hope that we provide enough 43 high-level matrix routines that few users will need it. 44 45 MatGetRow() always returns 0-based column indices, regardless of 46 whether the internal representation is 0-based (default) or 1-based. 47 48 For better efficiency, set cols and/or vals to PETSC_NULL if you do 49 not wish to extract these quantities. 50 51 The user can only examine the values extracted with MatGetRow(); 52 the values cannot be altered. To change the matrix entries, one 53 must use MatSetValues(). 54 55 You can only have one call to MatGetRow() outstanding for a particular 56 matrix at a time, per processor. MatGetRow() can only obtained rows 57 associated with the given processor, it cannot get rows from the 58 other processors; for that we suggest using MatGetSubMatrices(), then 59 MatGetRow() on the submatrix. The row indix passed to MatGetRows() 60 is in the global number of rows. 61 62 Fortran Notes: 63 The calling sequence from Fortran is 64 .vb 65 MatGetRow(matrix,row,ncols,cols,values,ierr) 66 Mat matrix (input) 67 integer row (input) 68 integer ncols (output) 69 integer cols(maxcols) (output) 70 double precision (or double complex) values(maxcols) output 71 .ve 72 where maxcols >= maximum nonzeros in any row of the matrix. 73 74 Caution: 75 Do not try to change the contents of the output arrays (cols and vals). 76 In some cases, this may corrupt the matrix. 77 78 Level: advanced 79 80 Concepts: matrices^row access 81 82 .seealso: MatRestoreRow(), MatSetValues(), MatGetValues(), MatGetSubmatrices(), MatGetDiagonal() 83 @*/ 84 int MatGetRow(Mat mat,int row,int *ncols,int *cols[],PetscScalar *vals[]) 85 { 86 int incols,ierr; 87 88 PetscFunctionBegin; 89 PetscValidHeaderSpecific(mat,MAT_COOKIE); 90 PetscValidType(mat); 91 MatPreallocated(mat); 92 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 93 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 94 if (!mat->ops->getrow) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 95 ierr = PetscLogEventBegin(MAT_GetRow,mat,0,0,0);CHKERRQ(ierr); 96 ierr = (*mat->ops->getrow)(mat,row,&incols,cols,vals);CHKERRQ(ierr); 97 if (ncols) *ncols = incols; 98 ierr = PetscLogEventEnd(MAT_GetRow,mat,0,0,0);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "MatRestoreRow" 104 /*@C 105 MatRestoreRow - Frees any temporary space allocated by MatGetRow(). 106 107 Not Collective 108 109 Input Parameters: 110 + mat - the matrix 111 . row - the row to get 112 . ncols, cols - the number of nonzeros and their columns 113 - vals - if nonzero the column values 114 115 Notes: 116 This routine should be called after you have finished examining the entries. 117 118 Fortran Notes: 119 The calling sequence from Fortran is 120 .vb 121 MatRestoreRow(matrix,row,ncols,cols,values,ierr) 122 Mat matrix (input) 123 integer row (input) 124 integer ncols (output) 125 integer cols(maxcols) (output) 126 double precision (or double complex) values(maxcols) output 127 .ve 128 Where maxcols >= maximum nonzeros in any row of the matrix. 129 130 In Fortran MatRestoreRow() MUST be called after MatGetRow() 131 before another call to MatGetRow() can be made. 132 133 Level: advanced 134 135 .seealso: MatGetRow() 136 @*/ 137 int MatRestoreRow(Mat mat,int row,int *ncols,int *cols[],PetscScalar *vals[]) 138 { 139 int ierr; 140 141 PetscFunctionBegin; 142 PetscValidHeaderSpecific(mat,MAT_COOKIE); 143 PetscValidIntPointer(ncols); 144 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 145 if (!mat->ops->restorerow) PetscFunctionReturn(0); 146 ierr = (*mat->ops->restorerow)(mat,row,ncols,cols,vals);CHKERRQ(ierr); 147 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "MatView" 153 /*@C 154 MatView - Visualizes a matrix object. 155 156 Collective on Mat 157 158 Input Parameters: 159 + mat - the matrix 160 - viewer - visualization context 161 162 Notes: 163 The available visualization contexts include 164 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 165 . PETSC_VIEWER_STDOUT_WORLD - synchronized standard 166 output where only the first processor opens 167 the file. All other processors send their 168 data to the first processor to print. 169 - PETSC_VIEWER_DRAW_WORLD - graphical display of nonzero structure 170 171 The user can open alternative visualization contexts with 172 + PetscViewerASCIIOpen() - Outputs matrix to a specified file 173 . PetscViewerBinaryOpen() - Outputs matrix in binary to a 174 specified file; corresponding input uses MatLoad() 175 . PetscViewerDrawOpen() - Outputs nonzero matrix structure to 176 an X window display 177 - PetscViewerSocketOpen() - Outputs matrix to Socket viewer. 178 Currently only the sequential dense and AIJ 179 matrix types support the Socket viewer. 180 181 The user can call PetscViewerSetFormat() to specify the output 182 format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF, 183 PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include 184 + PETSC_VIEWER_ASCII_DEFAULT - default, prints matrix contents 185 . PETSC_VIEWER_ASCII_MATLAB - prints matrix contents in Matlab format 186 . PETSC_VIEWER_ASCII_DENSE - prints entire matrix including zeros 187 . PETSC_VIEWER_ASCII_COMMON - prints matrix contents, using a sparse 188 format common among all matrix types 189 . PETSC_VIEWER_ASCII_IMPL - prints matrix contents, using an implementation-specific 190 format (which is in many cases the same as the default) 191 . PETSC_VIEWER_ASCII_INFO - prints basic information about the matrix 192 size and structure (not the matrix entries) 193 . PETSC_VIEWER_ASCII_INFO_DETAIL - prints more detailed information about 194 the matrix structure 195 196 Options Database Keys: 197 + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly() 198 . -mat_view_info_detailed - Prints more detailed info 199 . -mat_view - Prints matrix in ASCII format 200 . -mat_view_matlab - Prints matrix in Matlab format 201 . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX(). 202 . -display <name> - Sets display name (default is host) 203 . -draw_pause <sec> - Sets number of seconds to pause after display 204 . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual) 205 . -viewer_socket_machine <machine> 206 . -viewer_socket_port <port> 207 . -mat_view_binary - save matrix to file in binary format 208 - -viewer_binary_filename <name> 209 Level: beginner 210 211 Concepts: matrices^viewing 212 Concepts: matrices^plotting 213 Concepts: matrices^printing 214 215 .seealso: PetscViewerSetFormat(), PetscViewerASCIIOpen(), PetscViewerDrawOpen(), 216 PetscViewerSocketOpen(), PetscViewerBinaryOpen(), MatLoad() 217 @*/ 218 int MatView(Mat mat,PetscViewer viewer) 219 { 220 int ierr,rows,cols; 221 PetscTruth isascii; 222 char *cstr; 223 PetscViewerFormat format; 224 225 PetscFunctionBegin; 226 PetscValidHeaderSpecific(mat,MAT_COOKIE); 227 PetscValidType(mat); 228 MatPreallocated(mat); 229 if (!viewer) viewer = PETSC_VIEWER_STDOUT_(mat->comm); 230 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_COOKIE); 231 PetscCheckSameComm(mat,viewer); 232 if (!mat->assembled) SETERRQ(1,"Must call MatAssemblyBegin/End() before viewing matrix"); 233 234 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 235 if (isascii) { 236 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 237 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 238 if (mat->prefix) { 239 ierr = PetscViewerASCIIPrintf(viewer,"Matrix Object:(%s)\n",mat->prefix);CHKERRQ(ierr); 240 } else { 241 ierr = PetscViewerASCIIPrintf(viewer,"Matrix Object:\n");CHKERRQ(ierr); 242 } 243 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 244 ierr = MatGetType(mat,&cstr);CHKERRQ(ierr); 245 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPrintf(viewer,"type=%s, rows=%d, cols=%d\n",cstr,rows,cols);CHKERRQ(ierr); 247 if (mat->ops->getinfo) { 248 MatInfo info; 249 ierr = MatGetInfo(mat,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr); 250 ierr = PetscViewerASCIIPrintf(viewer,"total: nonzeros=%d, allocated nonzeros=%d\n", 251 (int)info.nz_used,(int)info.nz_allocated);CHKERRQ(ierr); 252 } 253 } 254 } 255 if (mat->ops->view) { 256 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 257 ierr = (*mat->ops->view)(mat,viewer);CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 259 } else if (!isascii) { 260 SETERRQ1(1,"Viewer type %s not supported",((PetscObject)viewer)->type_name); 261 } 262 if (isascii) { 263 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 264 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 265 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 266 } 267 } 268 PetscFunctionReturn(0); 269 } 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "MatScaleSystem" 273 /*@C 274 MatScaleSystem - Scale a vector solution and right hand side to 275 match the scaling of a scaled matrix. 276 277 Collective on Mat 278 279 Input Parameter: 280 + mat - the matrix 281 . x - solution vector (or PETSC_NULL) 282 - b - right hand side vector (or PETSC_NULL) 283 284 285 Notes: 286 For AIJ, BAIJ, and BDiag matrix formats, the matrices are not 287 internally scaled, so this does nothing. For MPIROWBS it 288 permutes and diagonally scales. 289 290 The SLES methods automatically call this routine when required 291 (via PCPreSolve()) so it is rarely used directly. 292 293 Level: Developer 294 295 Concepts: matrices^scaling 296 297 .seealso: MatUseScaledForm(), MatUnScaleSystem() 298 @*/ 299 int MatScaleSystem(Mat mat,Vec x,Vec b) 300 { 301 int ierr; 302 303 PetscFunctionBegin; 304 PetscValidHeaderSpecific(mat,MAT_COOKIE); 305 PetscValidType(mat); 306 MatPreallocated(mat); 307 if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);PetscCheckSameComm(mat,x);} 308 if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);PetscCheckSameComm(mat,b);} 309 310 if (mat->ops->scalesystem) { 311 ierr = (*mat->ops->scalesystem)(mat,x,b);CHKERRQ(ierr); 312 } 313 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 #undef __FUNCT__ 318 #define __FUNCT__ "MatUnScaleSystem" 319 /*@C 320 MatUnScaleSystem - Unscales a vector solution and right hand side to 321 match the original scaling of a scaled matrix. 322 323 Collective on Mat 324 325 Input Parameter: 326 + mat - the matrix 327 . x - solution vector (or PETSC_NULL) 328 - b - right hand side vector (or PETSC_NULL) 329 330 331 Notes: 332 For AIJ, BAIJ, and BDiag matrix formats, the matrices are not 333 internally scaled, so this does nothing. For MPIROWBS it 334 permutes and diagonally scales. 335 336 The SLES methods automatically call this routine when required 337 (via PCPreSolve()) so it is rarely used directly. 338 339 Level: Developer 340 341 .seealso: MatUseScaledForm(), MatScaleSystem() 342 @*/ 343 int MatUnScaleSystem(Mat mat,Vec x,Vec b) 344 { 345 int ierr; 346 347 PetscFunctionBegin; 348 PetscValidHeaderSpecific(mat,MAT_COOKIE); 349 PetscValidType(mat); 350 MatPreallocated(mat); 351 if (x) {PetscValidHeaderSpecific(x,VEC_COOKIE);PetscCheckSameComm(mat,x);} 352 if (b) {PetscValidHeaderSpecific(b,VEC_COOKIE);PetscCheckSameComm(mat,b);} 353 if (mat->ops->unscalesystem) { 354 ierr = (*mat->ops->unscalesystem)(mat,x,b);CHKERRQ(ierr); 355 } 356 PetscFunctionReturn(0); 357 } 358 359 #undef __FUNCT__ 360 #define __FUNCT__ "MatUseScaledForm" 361 /*@C 362 MatUseScaledForm - For matrix storage formats that scale the 363 matrix (for example MPIRowBS matrices are diagonally scaled on 364 assembly) indicates matrix operations (MatMult() etc) are 365 applied using the scaled matrix. 366 367 Collective on Mat 368 369 Input Parameter: 370 + mat - the matrix 371 - scaled - PETSC_TRUE for applying the scaled, PETSC_FALSE for 372 applying the original matrix 373 374 Notes: 375 For scaled matrix formats, applying the original, unscaled matrix 376 will be slightly more expensive 377 378 Level: Developer 379 380 .seealso: MatScaleSystem(), MatUnScaleSystem() 381 @*/ 382 int MatUseScaledForm(Mat mat,PetscTruth scaled) 383 { 384 int ierr; 385 386 PetscFunctionBegin; 387 PetscValidHeaderSpecific(mat,MAT_COOKIE); 388 PetscValidType(mat); 389 MatPreallocated(mat); 390 if (mat->ops->usescaledform) { 391 ierr = (*mat->ops->usescaledform)(mat,scaled);CHKERRQ(ierr); 392 } 393 PetscFunctionReturn(0); 394 } 395 396 #undef __FUNCT__ 397 #define __FUNCT__ "MatDestroy" 398 /*@C 399 MatDestroy - Frees space taken by a matrix. 400 401 Collective on Mat 402 403 Input Parameter: 404 . A - the matrix 405 406 Level: beginner 407 408 @*/ 409 int MatDestroy(Mat A) 410 { 411 int ierr; 412 413 PetscFunctionBegin; 414 PetscValidHeaderSpecific(A,MAT_COOKIE); 415 PetscValidType(A); 416 MatPreallocated(A); 417 if (--A->refct > 0) PetscFunctionReturn(0); 418 419 /* if memory was published with AMS then destroy it */ 420 ierr = PetscObjectDepublish(A);CHKERRQ(ierr); 421 if (A->mapping) { 422 ierr = ISLocalToGlobalMappingDestroy(A->mapping);CHKERRQ(ierr); 423 } 424 if (A->bmapping) { 425 ierr = ISLocalToGlobalMappingDestroy(A->bmapping);CHKERRQ(ierr); 426 } 427 if (A->rmap) { 428 ierr = PetscMapDestroy(A->rmap);CHKERRQ(ierr); 429 } 430 if (A->cmap) { 431 ierr = PetscMapDestroy(A->cmap);CHKERRQ(ierr); 432 } 433 434 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 435 PetscLogObjectDestroy(A); 436 PetscHeaderDestroy(A); 437 PetscFunctionReturn(0); 438 } 439 440 #undef __FUNCT__ 441 #define __FUNCT__ "MatValid" 442 /*@ 443 MatValid - Checks whether a matrix object is valid. 444 445 Collective on Mat 446 447 Input Parameter: 448 . m - the matrix to check 449 450 Output Parameter: 451 flg - flag indicating matrix status, either 452 PETSC_TRUE if matrix is valid, or PETSC_FALSE otherwise. 453 454 Level: developer 455 456 Concepts: matrices^validity 457 @*/ 458 int MatValid(Mat m,PetscTruth *flg) 459 { 460 PetscFunctionBegin; 461 PetscValidIntPointer(flg); 462 if (!m) *flg = PETSC_FALSE; 463 else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE; 464 else *flg = PETSC_TRUE; 465 PetscFunctionReturn(0); 466 } 467 468 #undef __FUNCT__ 469 #define __FUNCT__ "MatSetValues" 470 /*@ 471 MatSetValues - Inserts or adds a block of values into a matrix. 472 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 473 MUST be called after all calls to MatSetValues() have been completed. 474 475 Not Collective 476 477 Input Parameters: 478 + mat - the matrix 479 . v - a logically two-dimensional array of values 480 . m, idxm - the number of rows and their global indices 481 . n, idxn - the number of columns and their global indices 482 - addv - either ADD_VALUES or INSERT_VALUES, where 483 ADD_VALUES adds values to any existing entries, and 484 INSERT_VALUES replaces existing entries with new values 485 486 Notes: 487 By default the values, v, are row-oriented and unsorted. 488 See MatSetOption() for other options. 489 490 Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES 491 options cannot be mixed without intervening calls to the assembly 492 routines. 493 494 MatSetValues() uses 0-based row and column numbers in Fortran 495 as well as in C. 496 497 Negative indices may be passed in idxm and idxn, these rows and columns are 498 simply ignored. This allows easily inserting element stiffness matrices 499 with homogeneous Dirchlet boundary conditions that you don't want represented 500 in the matrix. 501 502 Efficiency Alert: 503 The routine MatSetValuesBlocked() may offer much better efficiency 504 for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ). 505 506 Level: beginner 507 508 Concepts: matrices^putting entries in 509 510 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal() 511 @*/ 512 int MatSetValues(Mat mat,int m,const int idxm[],int n,const int idxn[],const PetscScalar v[],InsertMode addv) 513 { 514 int ierr; 515 516 PetscFunctionBegin; 517 if (!m || !n) PetscFunctionReturn(0); /* no values to insert */ 518 PetscValidHeaderSpecific(mat,MAT_COOKIE); 519 PetscValidType(mat); 520 MatPreallocated(mat); 521 PetscValidIntPointer(idxm); 522 PetscValidIntPointer(idxn); 523 PetscValidScalarPointer(v); 524 if (mat->insertmode == NOT_SET_VALUES) { 525 mat->insertmode = addv; 526 } 527 #if defined(PETSC_USE_BOPT_g) 528 else if (mat->insertmode != addv) { 529 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values"); 530 } 531 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 532 #endif 533 534 if (mat->assembled) { 535 mat->was_assembled = PETSC_TRUE; 536 mat->assembled = PETSC_FALSE; 537 } 538 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 539 if (!mat->ops->setvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 540 ierr = (*mat->ops->setvalues)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr); 541 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 542 PetscFunctionReturn(0); 543 } 544 545 #undef __FUNCT__ 546 #define __FUNCT__ "MatSetValuesStencil" 547 /*@C 548 MatSetValuesStencil - Inserts or adds a block of values into a matrix. 549 Using structured grid indexing 550 551 Not Collective 552 553 Input Parameters: 554 + mat - the matrix 555 . v - a logically two-dimensional array of values 556 . m - number of rows being entered 557 . idxm - grid coordinates (and component number when dof > 1) for matrix rows being entered 558 . n - number of columns being entered 559 . idxn - grid coordinates (and component number when dof > 1) for matrix columns being entered 560 - addv - either ADD_VALUES or INSERT_VALUES, where 561 ADD_VALUES adds values to any existing entries, and 562 INSERT_VALUES replaces existing entries with new values 563 564 Notes: 565 By default the values, v, are row-oriented and unsorted. 566 See MatSetOption() for other options. 567 568 Calls to MatSetValuesStencil() with the INSERT_VALUES and ADD_VALUES 569 options cannot be mixed without intervening calls to the assembly 570 routines. 571 572 The grid coordinates are across the entire grid, not just the local portion 573 574 MatSetValuesStencil() uses 0-based row and column numbers in Fortran 575 as well as in C. 576 577 For setting/accessing vector values via array coordinates you can use the DAVecGetArray() routine 578 579 In order to use this routine you must either obtain the matrix with DAGetMatrix() 580 or call MatSetLocalToGlobalMapping() and MatSetStencil() first. 581 582 The columns and rows in the stencil passed in MUST be contained within the 583 ghost region of the given process as set with DACreateXXX() or MatSetStencil(). For example, 584 if you create a DA with an overlap of one grid level and on a particular process its first 585 local nonghost x logical coordinate is 6 (so its first ghost x logical coordinate is 5) the 586 first i index you can use in your column and row indices in MatSetStencil() is 5. 587 588 In Fortran idxm and idxn should be declared as 589 $ MatStencil idxm(4,m),idxn(4,n) 590 and the values inserted using 591 $ idxm(MatStencil_i,1) = i 592 $ idxm(MatStencil_j,1) = j 593 $ idxm(MatStencil_k,1) = k 594 $ idxm(MatStencil_c,1) = c 595 etc 596 597 Negative indices may be passed in idxm and idxn, these rows and columns are 598 simply ignored. This allows easily inserting element stiffness matrices 599 with homogeneous Dirchlet boundary conditions that you don't want represented 600 in the matrix. 601 602 Inspired by the structured grid interface to the HYPRE package 603 (http://www.llnl.gov/CASC/hypre) 604 605 Efficiency Alert: 606 The routine MatSetValuesBlockedStencil() may offer much better efficiency 607 for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ). 608 609 Level: beginner 610 611 Concepts: matrices^putting entries in 612 613 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal() 614 MatSetValues(), MatSetValuesBlockedStencil(), MatSetStencil(), DAGetMatrix(), DAVecGetArray() 615 @*/ 616 int MatSetValuesStencil(Mat mat,int m,const MatStencil idxm[],int n,const MatStencil idxn[],const PetscScalar v[],InsertMode addv) 617 { 618 int j,i,ierr,jdxm[128],jdxn[256],dim = mat->stencil.dim,*dims = mat->stencil.dims+1,tmp; 619 int *starts = mat->stencil.starts,*dxm = (int*)idxm,*dxn = (int*)idxn,sdim = dim - (1 - (int)mat->stencil.noc); 620 621 PetscFunctionBegin; 622 if (!m || !n) PetscFunctionReturn(0); /* no values to insert */ 623 PetscValidHeaderSpecific(mat,MAT_COOKIE); 624 PetscValidType(mat); 625 PetscValidIntPointer(idxm); 626 PetscValidIntPointer(idxn); 627 PetscValidScalarPointer(v); 628 629 if (m > 128) SETERRQ1(1,"Can only set 128 rows at a time; trying to set %d",m); 630 if (n > 128) SETERRQ1(1,"Can only set 256 columns at a time; trying to set %d",n); 631 632 for (i=0; i<m; i++) { 633 for (j=0; j<3-sdim; j++) dxm++; 634 tmp = *dxm++ - starts[0]; 635 for (j=0; j<dim-1; j++) { 636 tmp = tmp*dims[j] + *dxm++ - starts[j+1]; 637 } 638 if (mat->stencil.noc) dxm++; 639 jdxm[i] = tmp; 640 } 641 for (i=0; i<n; i++) { 642 for (j=0; j<3-sdim; j++) dxn++; 643 tmp = *dxn++ - starts[0]; 644 for (j=0; j<dim-1; j++) { 645 tmp = tmp*dims[j] + *dxn++ - starts[j+1]; 646 } 647 if (mat->stencil.noc) dxn++; 648 jdxn[i] = tmp; 649 } 650 ierr = MatSetValuesLocal(mat,m,jdxm,n,jdxn,v,addv);CHKERRQ(ierr); 651 PetscFunctionReturn(0); 652 } 653 654 #undef __FUNCT__ 655 #define __FUNCT__ "MatSetStencil" 656 /*@ 657 MatSetStencil - Sets the grid information for setting values into a matrix via 658 MatSetValuesStencil() 659 660 Not Collective 661 662 Input Parameters: 663 + mat - the matrix 664 . dim - dimension of the grid 1, 2, or 3 665 . dims - number of grid points in x, y, and z direction, including ghost points on your processor 666 . starts - starting point of ghost nodes on your processor in x, y, and z direction 667 - dof - number of degrees of freedom per node 668 669 670 Inspired by the structured grid interface to the HYPRE package 671 (www.llnl.gov/CASC/hyper) 672 673 For matrices generated with DAGetMatrix() this routine is automatically called and so not needed by the 674 user. 675 676 Level: beginner 677 678 Concepts: matrices^putting entries in 679 680 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal() 681 MatSetValues(), MatSetValuesBlockedStencil(), MatSetValuesStencil() 682 @*/ 683 int MatSetStencil(Mat mat,int dim,const int dims[],const int starts[],int dof) 684 { 685 int i; 686 687 PetscFunctionBegin; 688 PetscValidHeaderSpecific(mat,MAT_COOKIE); 689 PetscValidIntPointer(dims); 690 PetscValidIntPointer(starts); 691 692 mat->stencil.dim = dim + (dof > 1); 693 for (i=0; i<dim; i++) { 694 mat->stencil.dims[i] = dims[dim-i-1]; /* copy the values in backwards */ 695 mat->stencil.starts[i] = starts[dim-i-1]; 696 } 697 mat->stencil.dims[dim] = dof; 698 mat->stencil.starts[dim] = 0; 699 mat->stencil.noc = (PetscTruth)(dof == 1); 700 PetscFunctionReturn(0); 701 } 702 703 #undef __FUNCT__ 704 #define __FUNCT__ "MatSetValuesBlocked" 705 /*@ 706 MatSetValuesBlocked - Inserts or adds a block of values into a matrix. 707 708 Not Collective 709 710 Input Parameters: 711 + mat - the matrix 712 . v - a logically two-dimensional array of values 713 . m, idxm - the number of block rows and their global block indices 714 . n, idxn - the number of block columns and their global block indices 715 - addv - either ADD_VALUES or INSERT_VALUES, where 716 ADD_VALUES adds values to any existing entries, and 717 INSERT_VALUES replaces existing entries with new values 718 719 Notes: 720 By default the values, v, are row-oriented and unsorted. So the layout of 721 v is the same as for MatSetValues(). See MatSetOption() for other options. 722 723 Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES 724 options cannot be mixed without intervening calls to the assembly 725 routines. 726 727 MatSetValuesBlocked() uses 0-based row and column numbers in Fortran 728 as well as in C. 729 730 Negative indices may be passed in idxm and idxn, these rows and columns are 731 simply ignored. This allows easily inserting element stiffness matrices 732 with homogeneous Dirchlet boundary conditions that you don't want represented 733 in the matrix. 734 735 Each time an entry is set within a sparse matrix via MatSetValues(), 736 internal searching must be done to determine where to place the the 737 data in the matrix storage space. By instead inserting blocks of 738 entries via MatSetValuesBlocked(), the overhead of matrix assembly is 739 reduced. 740 741 Restrictions: 742 MatSetValuesBlocked() is currently supported only for the block AIJ 743 matrix format (MATSEQBAIJ and MATMPIBAIJ, which are created via 744 MatCreateSeqBAIJ() and MatCreateMPIBAIJ()). 745 746 Level: intermediate 747 748 Concepts: matrices^putting entries in blocked 749 750 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal() 751 @*/ 752 int MatSetValuesBlocked(Mat mat,int m,const int idxm[],int n,const int idxn[],const PetscScalar v[],InsertMode addv) 753 { 754 int ierr; 755 756 PetscFunctionBegin; 757 if (!m || !n) PetscFunctionReturn(0); /* no values to insert */ 758 PetscValidHeaderSpecific(mat,MAT_COOKIE); 759 PetscValidType(mat); 760 MatPreallocated(mat); 761 PetscValidIntPointer(idxm); 762 PetscValidIntPointer(idxn); 763 PetscValidScalarPointer(v); 764 if (mat->insertmode == NOT_SET_VALUES) { 765 mat->insertmode = addv; 766 } 767 #if defined(PETSC_USE_BOPT_g) 768 else if (mat->insertmode != addv) { 769 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values"); 770 } 771 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 772 #endif 773 774 if (mat->assembled) { 775 mat->was_assembled = PETSC_TRUE; 776 mat->assembled = PETSC_FALSE; 777 } 778 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 779 if (!mat->ops->setvaluesblocked) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 780 ierr = (*mat->ops->setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr); 781 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 782 PetscFunctionReturn(0); 783 } 784 785 #undef __FUNCT__ 786 #define __FUNCT__ "MatGetValues" 787 /*@ 788 MatGetValues - Gets a block of values from a matrix. 789 790 Not Collective; currently only returns a local block 791 792 Input Parameters: 793 + mat - the matrix 794 . v - a logically two-dimensional array for storing the values 795 . m, idxm - the number of rows and their global indices 796 - n, idxn - the number of columns and their global indices 797 798 Notes: 799 The user must allocate space (m*n PetscScalars) for the values, v. 800 The values, v, are then returned in a row-oriented format, 801 analogous to that used by default in MatSetValues(). 802 803 MatGetValues() uses 0-based row and column numbers in 804 Fortran as well as in C. 805 806 MatGetValues() requires that the matrix has been assembled 807 with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to 808 MatSetValues() and MatGetValues() CANNOT be made in succession 809 without intermediate matrix assembly. 810 811 Level: advanced 812 813 Concepts: matrices^accessing values 814 815 .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues() 816 @*/ 817 int MatGetValues(Mat mat,int m,const int idxm[],int n,const int idxn[],PetscScalar v[]) 818 { 819 int ierr; 820 821 PetscFunctionBegin; 822 PetscValidHeaderSpecific(mat,MAT_COOKIE); 823 PetscValidType(mat); 824 MatPreallocated(mat); 825 PetscValidIntPointer(idxm); 826 PetscValidIntPointer(idxn); 827 PetscValidScalarPointer(v); 828 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 829 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 830 if (!mat->ops->getvalues) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 831 832 ierr = PetscLogEventBegin(MAT_GetValues,mat,0,0,0);CHKERRQ(ierr); 833 ierr = (*mat->ops->getvalues)(mat,m,idxm,n,idxn,v);CHKERRQ(ierr); 834 ierr = PetscLogEventEnd(MAT_GetValues,mat,0,0,0);CHKERRQ(ierr); 835 PetscFunctionReturn(0); 836 } 837 838 #undef __FUNCT__ 839 #define __FUNCT__ "MatSetLocalToGlobalMapping" 840 /*@ 841 MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by 842 the routine MatSetValuesLocal() to allow users to insert matrix entries 843 using a local (per-processor) numbering. 844 845 Not Collective 846 847 Input Parameters: 848 + x - the matrix 849 - mapping - mapping created with ISLocalToGlobalMappingCreate() 850 or ISLocalToGlobalMappingCreateIS() 851 852 Level: intermediate 853 854 Concepts: matrices^local to global mapping 855 Concepts: local to global mapping^for matrices 856 857 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal() 858 @*/ 859 int MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping) 860 { 861 int ierr; 862 PetscFunctionBegin; 863 PetscValidHeaderSpecific(x,MAT_COOKIE); 864 PetscValidType(x); 865 MatPreallocated(x); 866 PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE); 867 if (x->mapping) { 868 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix"); 869 } 870 871 if (x->ops->setlocaltoglobalmapping) { 872 ierr = (*x->ops->setlocaltoglobalmapping)(x,mapping);CHKERRQ(ierr); 873 } else { 874 x->mapping = mapping; 875 ierr = PetscObjectReference((PetscObject)mapping);CHKERRQ(ierr); 876 } 877 PetscFunctionReturn(0); 878 } 879 880 #undef __FUNCT__ 881 #define __FUNCT__ "MatSetLocalToGlobalMappingBlock" 882 /*@ 883 MatSetLocalToGlobalMappingBlock - Sets a local-to-global numbering for use 884 by the routine MatSetValuesBlockedLocal() to allow users to insert matrix 885 entries using a local (per-processor) numbering. 886 887 Not Collective 888 889 Input Parameters: 890 + x - the matrix 891 - mapping - mapping created with ISLocalToGlobalMappingCreate() or 892 ISLocalToGlobalMappingCreateIS() 893 894 Level: intermediate 895 896 Concepts: matrices^local to global mapping blocked 897 Concepts: local to global mapping^for matrices, blocked 898 899 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(), 900 MatSetValuesBlocked(), MatSetValuesLocal() 901 @*/ 902 int MatSetLocalToGlobalMappingBlock(Mat x,ISLocalToGlobalMapping mapping) 903 { 904 int ierr; 905 PetscFunctionBegin; 906 PetscValidHeaderSpecific(x,MAT_COOKIE); 907 PetscValidType(x); 908 MatPreallocated(x); 909 PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE); 910 if (x->bmapping) { 911 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for matrix"); 912 } 913 914 x->bmapping = mapping; 915 ierr = PetscObjectReference((PetscObject)mapping);CHKERRQ(ierr); 916 PetscFunctionReturn(0); 917 } 918 919 #undef __FUNCT__ 920 #define __FUNCT__ "MatSetValuesLocal" 921 /*@ 922 MatSetValuesLocal - Inserts or adds values into certain locations of a matrix, 923 using a local ordering of the nodes. 924 925 Not Collective 926 927 Input Parameters: 928 + x - the matrix 929 . nrow, irow - number of rows and their local indices 930 . ncol, icol - number of columns and their local indices 931 . y - a logically two-dimensional array of values 932 - addv - either INSERT_VALUES or ADD_VALUES, where 933 ADD_VALUES adds values to any existing entries, and 934 INSERT_VALUES replaces existing entries with new values 935 936 Notes: 937 Before calling MatSetValuesLocal(), the user must first set the 938 local-to-global mapping by calling MatSetLocalToGlobalMapping(). 939 940 Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES 941 options cannot be mixed without intervening calls to the assembly 942 routines. 943 944 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 945 MUST be called after all calls to MatSetValuesLocal() have been completed. 946 947 Level: intermediate 948 949 Concepts: matrices^putting entries in with local numbering 950 951 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping(), 952 MatSetValueLocal() 953 @*/ 954 int MatSetValuesLocal(Mat mat,int nrow,const int irow[],int ncol,const int icol[],const PetscScalar y[],InsertMode addv) 955 { 956 int ierr,irowm[2048],icolm[2048]; 957 958 PetscFunctionBegin; 959 PetscValidHeaderSpecific(mat,MAT_COOKIE); 960 PetscValidType(mat); 961 MatPreallocated(mat); 962 PetscValidIntPointer(irow); 963 PetscValidIntPointer(icol); 964 PetscValidScalarPointer(y); 965 966 if (mat->insertmode == NOT_SET_VALUES) { 967 mat->insertmode = addv; 968 } 969 #if defined(PETSC_USE_BOPT_g) 970 else if (mat->insertmode != addv) { 971 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values"); 972 } 973 if (!mat->ops->setvalueslocal && (nrow > 2048 || ncol > 2048)) { 974 SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol); 975 } 976 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 977 #endif 978 979 if (mat->assembled) { 980 mat->was_assembled = PETSC_TRUE; 981 mat->assembled = PETSC_FALSE; 982 } 983 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 984 if (!mat->ops->setvalueslocal) { 985 ierr = ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm);CHKERRQ(ierr); 986 ierr = ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);CHKERRQ(ierr); 987 ierr = (*mat->ops->setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr); 988 } else { 989 ierr = (*mat->ops->setvalueslocal)(mat,nrow,irow,ncol,icol,y,addv);CHKERRQ(ierr); 990 } 991 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 992 PetscFunctionReturn(0); 993 } 994 995 #undef __FUNCT__ 996 #define __FUNCT__ "MatSetValuesBlockedLocal" 997 /*@ 998 MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix, 999 using a local ordering of the nodes a block at a time. 1000 1001 Not Collective 1002 1003 Input Parameters: 1004 + x - the matrix 1005 . nrow, irow - number of rows and their local indices 1006 . ncol, icol - number of columns and their local indices 1007 . y - a logically two-dimensional array of values 1008 - addv - either INSERT_VALUES or ADD_VALUES, where 1009 ADD_VALUES adds values to any existing entries, and 1010 INSERT_VALUES replaces existing entries with new values 1011 1012 Notes: 1013 Before calling MatSetValuesBlockedLocal(), the user must first set the 1014 local-to-global mapping by calling MatSetLocalToGlobalMappingBlock(), 1015 where the mapping MUST be set for matrix blocks, not for matrix elements. 1016 1017 Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES 1018 options cannot be mixed without intervening calls to the assembly 1019 routines. 1020 1021 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 1022 MUST be called after all calls to MatSetValuesBlockedLocal() have been completed. 1023 1024 Level: intermediate 1025 1026 Concepts: matrices^putting blocked values in with local numbering 1027 1028 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlock(), MatSetValuesBlocked() 1029 @*/ 1030 int MatSetValuesBlockedLocal(Mat mat,int nrow,const int irow[],int ncol,const int icol[],const PetscScalar y[],InsertMode addv) 1031 { 1032 int ierr,irowm[2048],icolm[2048]; 1033 1034 PetscFunctionBegin; 1035 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1036 PetscValidType(mat); 1037 MatPreallocated(mat); 1038 PetscValidIntPointer(irow); 1039 PetscValidIntPointer(icol); 1040 PetscValidScalarPointer(y); 1041 if (mat->insertmode == NOT_SET_VALUES) { 1042 mat->insertmode = addv; 1043 } 1044 #if defined(PETSC_USE_BOPT_g) 1045 else if (mat->insertmode != addv) { 1046 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Cannot mix add values and insert values"); 1047 } 1048 if (!mat->bmapping) { 1049 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with MatSetLocalToGlobalMappingBlock()"); 1050 } 1051 if (nrow > 2048 || ncol > 2048) { 1052 SETERRQ2(PETSC_ERR_SUP,"Number column/row indices must be <= 2048: are %d %d",nrow,ncol); 1053 } 1054 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1055 #endif 1056 1057 if (mat->assembled) { 1058 mat->was_assembled = PETSC_TRUE; 1059 mat->assembled = PETSC_FALSE; 1060 } 1061 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 1062 ierr = ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm);CHKERRQ(ierr); 1063 ierr = ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm);CHKERRQ(ierr); 1064 ierr = (*mat->ops->setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr); 1065 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 1066 PetscFunctionReturn(0); 1067 } 1068 1069 /* --------------------------------------------------------*/ 1070 #undef __FUNCT__ 1071 #define __FUNCT__ "MatMult" 1072 /*@ 1073 MatMult - Computes the matrix-vector product, y = Ax. 1074 1075 Collective on Mat and Vec 1076 1077 Input Parameters: 1078 + mat - the matrix 1079 - x - the vector to be multiplied 1080 1081 Output Parameters: 1082 . y - the result 1083 1084 Notes: 1085 The vectors x and y cannot be the same. I.e., one cannot 1086 call MatMult(A,y,y). 1087 1088 Level: beginner 1089 1090 Concepts: matrix-vector product 1091 1092 .seealso: MatMultTranspose(), MatMultAdd(), MatMultTransposeAdd() 1093 @*/ 1094 int MatMult(Mat mat,Vec x,Vec y) 1095 { 1096 int ierr; 1097 1098 PetscFunctionBegin; 1099 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1100 PetscValidType(mat); 1101 MatPreallocated(mat); 1102 PetscValidHeaderSpecific(x,VEC_COOKIE); 1103 PetscValidHeaderSpecific(y,VEC_COOKIE); 1104 1105 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1106 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1107 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1108 #ifndef PETSC_HAVE_CONSTRAINTS 1109 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1110 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 1111 if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n); 1112 #endif 1113 1114 if (mat->nullsp) { 1115 ierr = MatNullSpaceRemove(mat->nullsp,x,&x);CHKERRQ(ierr); 1116 } 1117 1118 ierr = PetscLogEventBegin(MAT_Mult,mat,x,y,0);CHKERRQ(ierr); 1119 ierr = (*mat->ops->mult)(mat,x,y);CHKERRQ(ierr); 1120 ierr = PetscLogEventEnd(MAT_Mult,mat,x,y,0);CHKERRQ(ierr); 1121 1122 if (mat->nullsp) { 1123 ierr = MatNullSpaceRemove(mat->nullsp,y,PETSC_NULL);CHKERRQ(ierr); 1124 } 1125 PetscFunctionReturn(0); 1126 } 1127 1128 #undef __FUNCT__ 1129 #define __FUNCT__ "MatMultTranspose" 1130 /*@ 1131 MatMultTranspose - Computes matrix transpose times a vector. 1132 1133 Collective on Mat and Vec 1134 1135 Input Parameters: 1136 + mat - the matrix 1137 - x - the vector to be multilplied 1138 1139 Output Parameters: 1140 . y - the result 1141 1142 Notes: 1143 The vectors x and y cannot be the same. I.e., one cannot 1144 call MatMultTranspose(A,y,y). 1145 1146 Level: beginner 1147 1148 Concepts: matrix vector product^transpose 1149 1150 .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd() 1151 @*/ 1152 int MatMultTranspose(Mat mat,Vec x,Vec y) 1153 { 1154 int ierr; 1155 PetscTruth flg1, flg2; 1156 1157 PetscFunctionBegin; 1158 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1159 PetscValidType(mat); 1160 MatPreallocated(mat); 1161 PetscValidHeaderSpecific(x,VEC_COOKIE); 1162 PetscValidHeaderSpecific(y,VEC_COOKIE); 1163 1164 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1165 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1166 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1167 #ifndef PETSC_HAVE_CONSTRAINTS 1168 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 1169 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N); 1170 #endif 1171 1172 if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP, "Operation not supported"); 1173 ierr = PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);CHKERRQ(ierr); 1174 if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined"); 1175 1176 ierr = PetscTypeCompare((PetscObject)mat,MATSEQSBAIJ,&flg1); 1177 ierr = PetscTypeCompare((PetscObject)mat,MATMPISBAIJ,&flg2); 1178 if (flg1 || flg2) { /* mat is in sbaij format */ 1179 ierr = (*mat->ops->mult)(mat,x,y);CHKERRQ(ierr); 1180 } else { 1181 ierr = (*mat->ops->multtranspose)(mat,x,y);CHKERRQ(ierr); 1182 } 1183 ierr = PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);CHKERRQ(ierr); 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNCT__ 1188 #define __FUNCT__ "MatMultAdd" 1189 /*@ 1190 MatMultAdd - Computes v3 = v2 + A * v1. 1191 1192 Collective on Mat and Vec 1193 1194 Input Parameters: 1195 + mat - the matrix 1196 - v1, v2 - the vectors 1197 1198 Output Parameters: 1199 . v3 - the result 1200 1201 Notes: 1202 The vectors v1 and v3 cannot be the same. I.e., one cannot 1203 call MatMultAdd(A,v1,v2,v1). 1204 1205 Level: beginner 1206 1207 Concepts: matrix vector product^addition 1208 1209 .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd() 1210 @*/ 1211 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3) 1212 { 1213 int ierr; 1214 1215 PetscFunctionBegin; 1216 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1217 PetscValidType(mat); 1218 MatPreallocated(mat); 1219 PetscValidHeaderSpecific(v1,VEC_COOKIE); 1220 PetscValidHeaderSpecific(v2,VEC_COOKIE); 1221 PetscValidHeaderSpecific(v3,VEC_COOKIE); 1222 1223 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1224 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1225 if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N); 1226 if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N); 1227 if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N); 1228 if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n); 1229 if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n); 1230 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors"); 1231 1232 ierr = PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1233 ierr = (*mat->ops->multadd)(mat,v1,v2,v3);CHKERRQ(ierr); 1234 ierr = PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1235 PetscFunctionReturn(0); 1236 } 1237 1238 #undef __FUNCT__ 1239 #define __FUNCT__ "MatMultTransposeAdd" 1240 /*@ 1241 MatMultTransposeAdd - Computes v3 = v2 + A' * v1. 1242 1243 Collective on Mat and Vec 1244 1245 Input Parameters: 1246 + mat - the matrix 1247 - v1, v2 - the vectors 1248 1249 Output Parameters: 1250 . v3 - the result 1251 1252 Notes: 1253 The vectors v1 and v3 cannot be the same. I.e., one cannot 1254 call MatMultTransposeAdd(A,v1,v2,v1). 1255 1256 Level: beginner 1257 1258 Concepts: matrix vector product^transpose and addition 1259 1260 .seealso: MatMultTranspose(), MatMultAdd(), MatMult() 1261 @*/ 1262 int MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3) 1263 { 1264 int ierr; 1265 1266 PetscFunctionBegin; 1267 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1268 PetscValidType(mat); 1269 MatPreallocated(mat); 1270 PetscValidHeaderSpecific(v1,VEC_COOKIE); 1271 PetscValidHeaderSpecific(v2,VEC_COOKIE); 1272 PetscValidHeaderSpecific(v3,VEC_COOKIE); 1273 1274 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1275 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1276 if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1277 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors"); 1278 if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N); 1279 if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N); 1280 if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N); 1281 1282 ierr = PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1283 ierr = (*mat->ops->multtransposeadd)(mat,v1,v2,v3);CHKERRQ(ierr); 1284 ierr = PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1285 PetscFunctionReturn(0); 1286 } 1287 1288 #undef __FUNCT__ 1289 #define __FUNCT__ "MatMultConstrained" 1290 /*@ 1291 MatMultConstrained - The inner multiplication routine for a 1292 constrained matrix P^T A P. 1293 1294 Collective on Mat and Vec 1295 1296 Input Parameters: 1297 + mat - the matrix 1298 - x - the vector to be multilplied 1299 1300 Output Parameters: 1301 . y - the result 1302 1303 Notes: 1304 The vectors x and y cannot be the same. I.e., one cannot 1305 call MatMult(A,y,y). 1306 1307 Level: beginner 1308 1309 .keywords: matrix, multiply, matrix-vector product, constraint 1310 .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd() 1311 @*/ 1312 int MatMultConstrained(Mat mat,Vec x,Vec y) 1313 { 1314 int ierr; 1315 1316 PetscFunctionBegin; 1317 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1318 PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1319 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1320 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1321 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1322 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1323 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 1324 if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n); 1325 1326 ierr = PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1327 ierr = (*mat->ops->multconstrained)(mat,x,y); CHKERRQ(ierr); 1328 ierr = PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1329 1330 PetscFunctionReturn(0); 1331 } 1332 1333 #undef __FUNCT__ 1334 #define __FUNCT__ "MatMultTransposeConstrained" 1335 /*@ 1336 MatMultTransposeConstrained - The inner multiplication routine for a 1337 constrained matrix P^T A^T P. 1338 1339 Collective on Mat and Vec 1340 1341 Input Parameters: 1342 + mat - the matrix 1343 - x - the vector to be multilplied 1344 1345 Output Parameters: 1346 . y - the result 1347 1348 Notes: 1349 The vectors x and y cannot be the same. I.e., one cannot 1350 call MatMult(A,y,y). 1351 1352 Level: beginner 1353 1354 .keywords: matrix, multiply, matrix-vector product, constraint 1355 .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd() 1356 @*/ 1357 int MatMultTransposeConstrained(Mat mat,Vec x,Vec y) 1358 { 1359 int ierr; 1360 1361 PetscFunctionBegin; 1362 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1363 PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1364 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1365 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1366 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1367 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1368 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 1369 1370 ierr = PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1371 ierr = (*mat->ops->multtransposeconstrained)(mat,x,y);CHKERRQ(ierr); 1372 ierr = PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1373 1374 PetscFunctionReturn(0); 1375 } 1376 /* ------------------------------------------------------------*/ 1377 #undef __FUNCT__ 1378 #define __FUNCT__ "MatGetInfo" 1379 /*@C 1380 MatGetInfo - Returns information about matrix storage (number of 1381 nonzeros, memory, etc.). 1382 1383 Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used 1384 as the flag 1385 1386 Input Parameters: 1387 . mat - the matrix 1388 1389 Output Parameters: 1390 + flag - flag indicating the type of parameters to be returned 1391 (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors, 1392 MAT_GLOBAL_SUM - sum over all processors) 1393 - info - matrix information context 1394 1395 Notes: 1396 The MatInfo context contains a variety of matrix data, including 1397 number of nonzeros allocated and used, number of mallocs during 1398 matrix assembly, etc. Additional information for factored matrices 1399 is provided (such as the fill ratio, number of mallocs during 1400 factorization, etc.). Much of this info is printed to STDOUT 1401 when using the runtime options 1402 $ -log_info -mat_view_info 1403 1404 Example for C/C++ Users: 1405 See the file ${PETSC_DIR}/include/petscmat.h for a complete list of 1406 data within the MatInfo context. For example, 1407 .vb 1408 MatInfo info; 1409 Mat A; 1410 double mal, nz_a, nz_u; 1411 1412 MatGetInfo(A,MAT_LOCAL,&info); 1413 mal = info.mallocs; 1414 nz_a = info.nz_allocated; 1415 .ve 1416 1417 Example for Fortran Users: 1418 Fortran users should declare info as a double precision 1419 array of dimension MAT_INFO_SIZE, and then extract the parameters 1420 of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h 1421 a complete list of parameter names. 1422 .vb 1423 double precision info(MAT_INFO_SIZE) 1424 double precision mal, nz_a 1425 Mat A 1426 integer ierr 1427 1428 call MatGetInfo(A,MAT_LOCAL,info,ierr) 1429 mal = info(MAT_INFO_MALLOCS) 1430 nz_a = info(MAT_INFO_NZ_ALLOCATED) 1431 .ve 1432 1433 Level: intermediate 1434 1435 Concepts: matrices^getting information on 1436 1437 @*/ 1438 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info) 1439 { 1440 int ierr; 1441 1442 PetscFunctionBegin; 1443 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1444 PetscValidType(mat); 1445 MatPreallocated(mat); 1446 PetscValidPointer(info); 1447 if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1448 ierr = (*mat->ops->getinfo)(mat,flag,info);CHKERRQ(ierr); 1449 PetscFunctionReturn(0); 1450 } 1451 1452 /* ----------------------------------------------------------*/ 1453 #undef __FUNCT__ 1454 #define __FUNCT__ "MatILUDTFactor" 1455 /*@C 1456 MatILUDTFactor - Performs a drop tolerance ILU factorization. 1457 1458 Collective on Mat 1459 1460 Input Parameters: 1461 + mat - the matrix 1462 . info - information about the factorization to be done 1463 . row - row permutation 1464 - col - column permutation 1465 1466 Output Parameters: 1467 . fact - the factored matrix 1468 1469 Level: developer 1470 1471 Notes: 1472 Most users should employ the simplified SLES interface for linear solvers 1473 instead of working directly with matrix algebra routines such as this. 1474 See, e.g., SLESCreate(). 1475 1476 This is currently only supported for the SeqAIJ matrix format using code 1477 from Yousef Saad's SPARSEKIT2 package (translated to C with f2c) and/or 1478 Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright 1479 and thus can be distributed with PETSc. 1480 1481 Concepts: matrices^ILUDT factorization 1482 1483 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo 1484 @*/ 1485 int MatILUDTFactor(Mat mat,MatFactorInfo *info,IS row,IS col,Mat *fact) 1486 { 1487 int ierr; 1488 1489 PetscFunctionBegin; 1490 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1491 PetscValidType(mat); 1492 MatPreallocated(mat); 1493 PetscValidPointer(fact); 1494 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1495 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1496 if (!mat->ops->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1497 1498 ierr = PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1499 ierr = (*mat->ops->iludtfactor)(mat,info,row,col,fact);CHKERRQ(ierr); 1500 ierr = PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1501 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1502 1503 PetscFunctionReturn(0); 1504 } 1505 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "MatLUFactor" 1508 /*@ 1509 MatLUFactor - Performs in-place LU factorization of matrix. 1510 1511 Collective on Mat 1512 1513 Input Parameters: 1514 + mat - the matrix 1515 . row - row permutation 1516 . col - column permutation 1517 - info - options for factorization, includes 1518 $ fill - expected fill as ratio of original fill. 1519 $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting) 1520 $ Run with the option -log_info to determine an optimal value to use 1521 1522 Notes: 1523 Most users should employ the simplified SLES interface for linear solvers 1524 instead of working directly with matrix algebra routines such as this. 1525 See, e.g., SLESCreate(). 1526 1527 This changes the state of the matrix to a factored matrix; it cannot be used 1528 for example with MatSetValues() unless one first calls MatSetUnfactored(). 1529 1530 Level: developer 1531 1532 Concepts: matrices^LU factorization 1533 1534 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), 1535 MatGetOrdering(), MatSetUnfactored(), MatFactorInfo 1536 1537 @*/ 1538 int MatLUFactor(Mat mat,IS row,IS col,MatFactorInfo *info) 1539 { 1540 int ierr; 1541 1542 PetscFunctionBegin; 1543 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1544 PetscValidType(mat); 1545 MatPreallocated(mat); 1546 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1547 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1548 if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1549 1550 ierr = PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);CHKERRQ(ierr); 1551 ierr = (*mat->ops->lufactor)(mat,row,col,info);CHKERRQ(ierr); 1552 ierr = PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);CHKERRQ(ierr); 1553 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 1554 PetscFunctionReturn(0); 1555 } 1556 1557 #undef __FUNCT__ 1558 #define __FUNCT__ "MatILUFactor" 1559 /*@ 1560 MatILUFactor - Performs in-place ILU factorization of matrix. 1561 1562 Collective on Mat 1563 1564 Input Parameters: 1565 + mat - the matrix 1566 . row - row permutation 1567 . col - column permutation 1568 - info - structure containing 1569 $ levels - number of levels of fill. 1570 $ expected fill - as ratio of original fill. 1571 $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices 1572 missing diagonal entries) 1573 1574 Notes: 1575 Probably really in-place only when level of fill is zero, otherwise allocates 1576 new space to store factored matrix and deletes previous memory. 1577 1578 Most users should employ the simplified SLES interface for linear solvers 1579 instead of working directly with matrix algebra routines such as this. 1580 See, e.g., SLESCreate(). 1581 1582 Level: developer 1583 1584 Concepts: matrices^ILU factorization 1585 1586 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo 1587 @*/ 1588 int MatILUFactor(Mat mat,IS row,IS col,MatFactorInfo *info) 1589 { 1590 int ierr; 1591 1592 PetscFunctionBegin; 1593 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1594 PetscValidType(mat); 1595 MatPreallocated(mat); 1596 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 1597 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1598 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1599 if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1600 1601 ierr = PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1602 ierr = (*mat->ops->ilufactor)(mat,row,col,info);CHKERRQ(ierr); 1603 ierr = PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1604 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 1605 PetscFunctionReturn(0); 1606 } 1607 1608 #undef __FUNCT__ 1609 #define __FUNCT__ "MatLUFactorSymbolic" 1610 /*@ 1611 MatLUFactorSymbolic - Performs symbolic LU factorization of matrix. 1612 Call this routine before calling MatLUFactorNumeric(). 1613 1614 Collective on Mat 1615 1616 Input Parameters: 1617 + mat - the matrix 1618 . row, col - row and column permutations 1619 - info - options for factorization, includes 1620 $ fill - expected fill as ratio of original fill. 1621 $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting) 1622 $ Run with the option -log_info to determine an optimal value to use 1623 1624 Output Parameter: 1625 . fact - new matrix that has been symbolically factored 1626 1627 Notes: 1628 See the users manual for additional information about 1629 choosing the fill factor for better efficiency. 1630 1631 Most users should employ the simplified SLES interface for linear solvers 1632 instead of working directly with matrix algebra routines such as this. 1633 See, e.g., SLESCreate(). 1634 1635 Level: developer 1636 1637 Concepts: matrices^LU symbolic factorization 1638 1639 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo 1640 @*/ 1641 int MatLUFactorSymbolic(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact) 1642 { 1643 int ierr; 1644 1645 PetscFunctionBegin; 1646 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1647 PetscValidType(mat); 1648 MatPreallocated(mat); 1649 PetscValidPointer(fact); 1650 PetscValidHeaderSpecific(row,IS_COOKIE); 1651 PetscValidHeaderSpecific(col,IS_COOKIE); 1652 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1653 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1654 if (!mat->ops->lufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic LU",mat->type_name); 1655 1656 ierr = PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 1657 ierr = (*mat->ops->lufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr); 1658 ierr = PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 1659 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1660 PetscFunctionReturn(0); 1661 } 1662 1663 #undef __FUNCT__ 1664 #define __FUNCT__ "MatLUFactorNumeric" 1665 /*@ 1666 MatLUFactorNumeric - Performs numeric LU factorization of a matrix. 1667 Call this routine after first calling MatLUFactorSymbolic(). 1668 1669 Collective on Mat 1670 1671 Input Parameters: 1672 + mat - the matrix 1673 - fact - the matrix generated for the factor, from MatLUFactorSymbolic() 1674 1675 Notes: 1676 See MatLUFactor() for in-place factorization. See 1677 MatCholeskyFactorNumeric() for the symmetric, positive definite case. 1678 1679 Most users should employ the simplified SLES interface for linear solvers 1680 instead of working directly with matrix algebra routines such as this. 1681 See, e.g., SLESCreate(). 1682 1683 Level: developer 1684 1685 Concepts: matrices^LU numeric factorization 1686 1687 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor() 1688 @*/ 1689 int MatLUFactorNumeric(Mat mat,Mat *fact) 1690 { 1691 int ierr; 1692 1693 PetscFunctionBegin; 1694 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1695 PetscValidType(mat); 1696 MatPreallocated(mat); 1697 PetscValidPointer(fact); 1698 PetscValidHeaderSpecific(*fact,MAT_COOKIE); 1699 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1700 if (mat->M != (*fact)->M || mat->N != (*fact)->N) { 1701 SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d", 1702 mat->M,(*fact)->M,mat->N,(*fact)->N); 1703 } 1704 if (!(*fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1705 1706 ierr = PetscLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1707 ierr = (*(*fact)->ops->lufactornumeric)(mat,fact);CHKERRQ(ierr); 1708 ierr = PetscLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1709 1710 ierr = MatView_Private(*fact);CHKERRQ(ierr); 1711 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1712 PetscFunctionReturn(0); 1713 } 1714 1715 #undef __FUNCT__ 1716 #define __FUNCT__ "MatCholeskyFactor" 1717 /*@ 1718 MatCholeskyFactor - Performs in-place Cholesky factorization of a 1719 symmetric matrix. 1720 1721 Collective on Mat 1722 1723 Input Parameters: 1724 + mat - the matrix 1725 . perm - row and column permutations 1726 - f - expected fill as ratio of original fill 1727 1728 Notes: 1729 See MatLUFactor() for the nonsymmetric case. See also 1730 MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric(). 1731 1732 Most users should employ the simplified SLES interface for linear solvers 1733 instead of working directly with matrix algebra routines such as this. 1734 See, e.g., SLESCreate(). 1735 1736 Level: developer 1737 1738 Concepts: matrices^Cholesky factorization 1739 1740 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric() 1741 MatGetOrdering() 1742 1743 @*/ 1744 int MatCholeskyFactor(Mat mat,IS perm,MatFactorInfo *info) 1745 { 1746 int ierr; 1747 1748 PetscFunctionBegin; 1749 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1750 PetscValidType(mat); 1751 MatPreallocated(mat); 1752 PetscValidHeaderSpecific(perm,IS_COOKIE); 1753 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 1754 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1755 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1756 if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1757 1758 ierr = PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);CHKERRQ(ierr); 1759 ierr = (*mat->ops->choleskyfactor)(mat,perm,info);CHKERRQ(ierr); 1760 ierr = PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);CHKERRQ(ierr); 1761 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 1762 PetscFunctionReturn(0); 1763 } 1764 1765 #undef __FUNCT__ 1766 #define __FUNCT__ "MatCholeskyFactorSymbolic" 1767 /*@ 1768 MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization 1769 of a symmetric matrix. 1770 1771 Collective on Mat 1772 1773 Input Parameters: 1774 + mat - the matrix 1775 . perm - row and column permutations 1776 - info - options for factorization, includes 1777 $ fill - expected fill as ratio of original fill. 1778 $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting) 1779 $ Run with the option -log_info to determine an optimal value to use 1780 1781 Output Parameter: 1782 . fact - the factored matrix 1783 1784 Notes: 1785 See MatLUFactorSymbolic() for the nonsymmetric case. See also 1786 MatCholeskyFactor() and MatCholeskyFactorNumeric(). 1787 1788 Most users should employ the simplified SLES interface for linear solvers 1789 instead of working directly with matrix algebra routines such as this. 1790 See, e.g., SLESCreate(). 1791 1792 Level: developer 1793 1794 Concepts: matrices^Cholesky symbolic factorization 1795 1796 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric() 1797 MatGetOrdering() 1798 1799 @*/ 1800 int MatCholeskyFactorSymbolic(Mat mat,IS perm,MatFactorInfo *info,Mat *fact) 1801 { 1802 int ierr; 1803 1804 PetscFunctionBegin; 1805 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1806 PetscValidType(mat); 1807 MatPreallocated(mat); 1808 PetscValidPointer(fact); 1809 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 1810 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1811 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1812 if (!mat->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1813 1814 ierr = PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 1815 ierr = (*mat->ops->choleskyfactorsymbolic)(mat,perm,info,fact);CHKERRQ(ierr); 1816 ierr = PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 1817 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1818 PetscFunctionReturn(0); 1819 } 1820 1821 #undef __FUNCT__ 1822 #define __FUNCT__ "MatCholeskyFactorNumeric" 1823 /*@ 1824 MatCholeskyFactorNumeric - Performs numeric Cholesky factorization 1825 of a symmetric matrix. Call this routine after first calling 1826 MatCholeskyFactorSymbolic(). 1827 1828 Collective on Mat 1829 1830 Input Parameter: 1831 . mat - the initial matrix 1832 1833 Output Parameter: 1834 . fact - the factored matrix 1835 1836 Notes: 1837 Most users should employ the simplified SLES interface for linear solvers 1838 instead of working directly with matrix algebra routines such as this. 1839 See, e.g., SLESCreate(). 1840 1841 Level: developer 1842 1843 Concepts: matrices^Cholesky numeric factorization 1844 1845 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric() 1846 @*/ 1847 int MatCholeskyFactorNumeric(Mat mat,Mat *fact) 1848 { 1849 int ierr; 1850 1851 PetscFunctionBegin; 1852 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1853 PetscValidType(mat); 1854 MatPreallocated(mat); 1855 PetscValidPointer(fact); 1856 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1857 if (!(*fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1858 if (mat->M != (*fact)->M || mat->N != (*fact)->N) { 1859 SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d", 1860 mat->M,(*fact)->M,mat->N,(*fact)->N); 1861 } 1862 1863 ierr = PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1864 ierr = (*(*fact)->ops->choleskyfactornumeric)(mat,fact);CHKERRQ(ierr); 1865 ierr = PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1866 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1867 PetscFunctionReturn(0); 1868 } 1869 1870 /* ----------------------------------------------------------------*/ 1871 #undef __FUNCT__ 1872 #define __FUNCT__ "MatSolve" 1873 /*@ 1874 MatSolve - Solves A x = b, given a factored matrix. 1875 1876 Collective on Mat and Vec 1877 1878 Input Parameters: 1879 + mat - the factored matrix 1880 - b - the right-hand-side vector 1881 1882 Output Parameter: 1883 . x - the result vector 1884 1885 Notes: 1886 The vectors b and x cannot be the same. I.e., one cannot 1887 call MatSolve(A,x,x). 1888 1889 Notes: 1890 Most users should employ the simplified SLES interface for linear solvers 1891 instead of working directly with matrix algebra routines such as this. 1892 See, e.g., SLESCreate(). 1893 1894 Level: developer 1895 1896 Concepts: matrices^triangular solves 1897 1898 .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd() 1899 @*/ 1900 int MatSolve(Mat mat,Vec b,Vec x) 1901 { 1902 int ierr; 1903 1904 PetscFunctionBegin; 1905 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1906 PetscValidType(mat); 1907 MatPreallocated(mat); 1908 PetscValidHeaderSpecific(b,VEC_COOKIE); 1909 PetscValidHeaderSpecific(x,VEC_COOKIE); 1910 PetscCheckSameComm(mat,b); 1911 PetscCheckSameComm(mat,x); 1912 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 1913 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 1914 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1915 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1916 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1917 if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0); 1918 1919 if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1920 ierr = PetscLogEventBegin(MAT_Solve,mat,b,x,0);CHKERRQ(ierr); 1921 ierr = (*mat->ops->solve)(mat,b,x);CHKERRQ(ierr); 1922 ierr = PetscLogEventEnd(MAT_Solve,mat,b,x,0);CHKERRQ(ierr); 1923 PetscFunctionReturn(0); 1924 } 1925 1926 #undef __FUNCT__ 1927 #define __FUNCT__ "MatForwardSolve" 1928 /* @ 1929 MatForwardSolve - Solves L x = b, given a factored matrix, A = LU. 1930 1931 Collective on Mat and Vec 1932 1933 Input Parameters: 1934 + mat - the factored matrix 1935 - b - the right-hand-side vector 1936 1937 Output Parameter: 1938 . x - the result vector 1939 1940 Notes: 1941 MatSolve() should be used for most applications, as it performs 1942 a forward solve followed by a backward solve. 1943 1944 The vectors b and x cannot be the same. I.e., one cannot 1945 call MatForwardSolve(A,x,x). 1946 1947 Most users should employ the simplified SLES interface for linear solvers 1948 instead of working directly with matrix algebra routines such as this. 1949 See, e.g., SLESCreate(). 1950 1951 Level: developer 1952 1953 Concepts: matrices^forward solves 1954 1955 .seealso: MatSolve(), MatBackwardSolve() 1956 @ */ 1957 int MatForwardSolve(Mat mat,Vec b,Vec x) 1958 { 1959 int ierr; 1960 1961 PetscFunctionBegin; 1962 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1963 PetscValidType(mat); 1964 MatPreallocated(mat); 1965 PetscValidHeaderSpecific(b,VEC_COOKIE); 1966 PetscValidHeaderSpecific(x,VEC_COOKIE); 1967 PetscCheckSameComm(mat,b); 1968 PetscCheckSameComm(mat,x); 1969 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 1970 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 1971 if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1972 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1973 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1974 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1975 1976 ierr = PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);CHKERRQ(ierr); 1977 ierr = (*mat->ops->forwardsolve)(mat,b,x);CHKERRQ(ierr); 1978 ierr = PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);CHKERRQ(ierr); 1979 PetscFunctionReturn(0); 1980 } 1981 1982 #undef __FUNCT__ 1983 #define __FUNCT__ "MatBackwardSolve" 1984 /* @ 1985 MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU. 1986 1987 Collective on Mat and Vec 1988 1989 Input Parameters: 1990 + mat - the factored matrix 1991 - b - the right-hand-side vector 1992 1993 Output Parameter: 1994 . x - the result vector 1995 1996 Notes: 1997 MatSolve() should be used for most applications, as it performs 1998 a forward solve followed by a backward solve. 1999 2000 The vectors b and x cannot be the same. I.e., one cannot 2001 call MatBackwardSolve(A,x,x). 2002 2003 Most users should employ the simplified SLES interface for linear solvers 2004 instead of working directly with matrix algebra routines such as this. 2005 See, e.g., SLESCreate(). 2006 2007 Level: developer 2008 2009 Concepts: matrices^backward solves 2010 2011 .seealso: MatSolve(), MatForwardSolve() 2012 @ */ 2013 int MatBackwardSolve(Mat mat,Vec b,Vec x) 2014 { 2015 int ierr; 2016 2017 PetscFunctionBegin; 2018 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2019 PetscValidType(mat); 2020 MatPreallocated(mat); 2021 PetscValidHeaderSpecific(b,VEC_COOKIE); 2022 PetscValidHeaderSpecific(x,VEC_COOKIE); 2023 PetscCheckSameComm(mat,b); 2024 PetscCheckSameComm(mat,x); 2025 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2026 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2027 if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2028 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2029 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2030 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2031 2032 ierr = PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);CHKERRQ(ierr); 2033 ierr = (*mat->ops->backwardsolve)(mat,b,x);CHKERRQ(ierr); 2034 ierr = PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);CHKERRQ(ierr); 2035 PetscFunctionReturn(0); 2036 } 2037 2038 #undef __FUNCT__ 2039 #define __FUNCT__ "MatSolveAdd" 2040 /*@ 2041 MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix. 2042 2043 Collective on Mat and Vec 2044 2045 Input Parameters: 2046 + mat - the factored matrix 2047 . b - the right-hand-side vector 2048 - y - the vector to be added to 2049 2050 Output Parameter: 2051 . x - the result vector 2052 2053 Notes: 2054 The vectors b and x cannot be the same. I.e., one cannot 2055 call MatSolveAdd(A,x,y,x). 2056 2057 Most users should employ the simplified SLES interface for linear solvers 2058 instead of working directly with matrix algebra routines such as this. 2059 See, e.g., SLESCreate(). 2060 2061 Level: developer 2062 2063 Concepts: matrices^triangular solves 2064 2065 .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd() 2066 @*/ 2067 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x) 2068 { 2069 PetscScalar one = 1.0; 2070 Vec tmp; 2071 int ierr; 2072 2073 PetscFunctionBegin; 2074 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2075 PetscValidType(mat); 2076 MatPreallocated(mat); 2077 PetscValidHeaderSpecific(y,VEC_COOKIE); 2078 PetscValidHeaderSpecific(b,VEC_COOKIE); 2079 PetscValidHeaderSpecific(x,VEC_COOKIE); 2080 PetscCheckSameComm(mat,b); 2081 PetscCheckSameComm(mat,y); 2082 PetscCheckSameComm(mat,x); 2083 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2084 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2085 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2086 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2087 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 2088 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2089 if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n); 2090 2091 ierr = PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);CHKERRQ(ierr); 2092 if (mat->ops->solveadd) { 2093 ierr = (*mat->ops->solveadd)(mat,b,y,x);CHKERRQ(ierr); 2094 } else { 2095 /* do the solve then the add manually */ 2096 if (x != y) { 2097 ierr = MatSolve(mat,b,x);CHKERRQ(ierr); 2098 ierr = VecAXPY(&one,y,x);CHKERRQ(ierr); 2099 } else { 2100 ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr); 2101 PetscLogObjectParent(mat,tmp); 2102 ierr = VecCopy(x,tmp);CHKERRQ(ierr); 2103 ierr = MatSolve(mat,b,x);CHKERRQ(ierr); 2104 ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr); 2105 ierr = VecDestroy(tmp);CHKERRQ(ierr); 2106 } 2107 } 2108 ierr = PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);CHKERRQ(ierr); 2109 PetscFunctionReturn(0); 2110 } 2111 2112 #undef __FUNCT__ 2113 #define __FUNCT__ "MatSolveTranspose" 2114 /*@ 2115 MatSolveTranspose - Solves A' x = b, given a factored matrix. 2116 2117 Collective on Mat and Vec 2118 2119 Input Parameters: 2120 + mat - the factored matrix 2121 - b - the right-hand-side vector 2122 2123 Output Parameter: 2124 . x - the result vector 2125 2126 Notes: 2127 The vectors b and x cannot be the same. I.e., one cannot 2128 call MatSolveTranspose(A,x,x). 2129 2130 Most users should employ the simplified SLES interface for linear solvers 2131 instead of working directly with matrix algebra routines such as this. 2132 See, e.g., SLESCreate(). 2133 2134 Level: developer 2135 2136 Concepts: matrices^triangular solves 2137 2138 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd() 2139 @*/ 2140 int MatSolveTranspose(Mat mat,Vec b,Vec x) 2141 { 2142 int ierr; 2143 2144 PetscFunctionBegin; 2145 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2146 PetscValidType(mat); 2147 MatPreallocated(mat); 2148 PetscValidHeaderSpecific(b,VEC_COOKIE); 2149 PetscValidHeaderSpecific(x,VEC_COOKIE); 2150 PetscCheckSameComm(mat,b); 2151 PetscCheckSameComm(mat,x); 2152 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2153 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2154 if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",mat->type_name); 2155 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 2156 if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N); 2157 2158 ierr = PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);CHKERRQ(ierr); 2159 ierr = (*mat->ops->solvetranspose)(mat,b,x);CHKERRQ(ierr); 2160 ierr = PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);CHKERRQ(ierr); 2161 PetscFunctionReturn(0); 2162 } 2163 2164 #undef __FUNCT__ 2165 #define __FUNCT__ "MatSolveTransposeAdd" 2166 /*@ 2167 MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a 2168 factored matrix. 2169 2170 Collective on Mat and Vec 2171 2172 Input Parameters: 2173 + mat - the factored matrix 2174 . b - the right-hand-side vector 2175 - y - the vector to be added to 2176 2177 Output Parameter: 2178 . x - the result vector 2179 2180 Notes: 2181 The vectors b and x cannot be the same. I.e., one cannot 2182 call MatSolveTransposeAdd(A,x,y,x). 2183 2184 Most users should employ the simplified SLES interface for linear solvers 2185 instead of working directly with matrix algebra routines such as this. 2186 See, e.g., SLESCreate(). 2187 2188 Level: developer 2189 2190 Concepts: matrices^triangular solves 2191 2192 .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose() 2193 @*/ 2194 int MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x) 2195 { 2196 PetscScalar one = 1.0; 2197 int ierr; 2198 Vec tmp; 2199 2200 PetscFunctionBegin; 2201 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2202 PetscValidType(mat); 2203 MatPreallocated(mat); 2204 PetscValidHeaderSpecific(y,VEC_COOKIE); 2205 PetscValidHeaderSpecific(b,VEC_COOKIE); 2206 PetscValidHeaderSpecific(x,VEC_COOKIE); 2207 PetscCheckSameComm(mat,b); 2208 PetscCheckSameComm(mat,y); 2209 PetscCheckSameComm(mat,x); 2210 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2211 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2212 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 2213 if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N); 2214 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N); 2215 if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n); 2216 2217 ierr = PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);CHKERRQ(ierr); 2218 if (mat->ops->solvetransposeadd) { 2219 ierr = (*mat->ops->solvetransposeadd)(mat,b,y,x);CHKERRQ(ierr); 2220 } else { 2221 /* do the solve then the add manually */ 2222 if (x != y) { 2223 ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr); 2224 ierr = VecAXPY(&one,y,x);CHKERRQ(ierr); 2225 } else { 2226 ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr); 2227 PetscLogObjectParent(mat,tmp); 2228 ierr = VecCopy(x,tmp);CHKERRQ(ierr); 2229 ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr); 2230 ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr); 2231 ierr = VecDestroy(tmp);CHKERRQ(ierr); 2232 } 2233 } 2234 ierr = PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);CHKERRQ(ierr); 2235 PetscFunctionReturn(0); 2236 } 2237 /* ----------------------------------------------------------------*/ 2238 2239 #undef __FUNCT__ 2240 #define __FUNCT__ "MatRelax" 2241 /*@ 2242 MatRelax - Computes one relaxation sweep. 2243 2244 Collective on Mat and Vec 2245 2246 Input Parameters: 2247 + mat - the matrix 2248 . b - the right hand side 2249 . omega - the relaxation factor 2250 . flag - flag indicating the type of SOR (see below) 2251 . shift - diagonal shift 2252 - its - the number of iterations 2253 - lits - the number of local iterations 2254 2255 Output Parameters: 2256 . x - the solution (can contain an initial guess) 2257 2258 SOR Flags: 2259 . SOR_FORWARD_SWEEP - forward SOR 2260 . SOR_BACKWARD_SWEEP - backward SOR 2261 . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR) 2262 . SOR_LOCAL_FORWARD_SWEEP - local forward SOR 2263 . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR 2264 . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR 2265 . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies 2266 upper/lower triangular part of matrix to 2267 vector (with omega) 2268 . SOR_ZERO_INITIAL_GUESS - zero initial guess 2269 2270 Notes: 2271 SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and 2272 SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings 2273 on each processor. 2274 2275 Application programmers will not generally use MatRelax() directly, 2276 but instead will employ the SLES/PC interface. 2277 2278 Notes for Advanced Users: 2279 The flags are implemented as bitwise inclusive or operations. 2280 For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP) 2281 to specify a zero initial guess for SSOR. 2282 2283 Most users should employ the simplified SLES interface for linear solvers 2284 instead of working directly with matrix algebra routines such as this. 2285 See, e.g., SLESCreate(). 2286 2287 Level: developer 2288 2289 Concepts: matrices^relaxation 2290 Concepts: matrices^SOR 2291 Concepts: matrices^Gauss-Seidel 2292 2293 @*/ 2294 int MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,int its,int lits,Vec x) 2295 { 2296 int ierr; 2297 2298 PetscFunctionBegin; 2299 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2300 PetscValidType(mat); 2301 MatPreallocated(mat); 2302 PetscValidHeaderSpecific(b,VEC_COOKIE); 2303 PetscValidHeaderSpecific(x,VEC_COOKIE); 2304 PetscCheckSameComm(mat,b); 2305 PetscCheckSameComm(mat,x); 2306 if (!mat->ops->relax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2307 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2308 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2309 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2310 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2311 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2312 2313 ierr = PetscLogEventBegin(MAT_Relax,mat,b,x,0);CHKERRQ(ierr); 2314 ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);CHKERRQ(ierr); 2315 ierr = PetscLogEventEnd(MAT_Relax,mat,b,x,0);CHKERRQ(ierr); 2316 PetscFunctionReturn(0); 2317 } 2318 2319 #undef __FUNCT__ 2320 #define __FUNCT__ "MatCopy_Basic" 2321 /* 2322 Default matrix copy routine. 2323 */ 2324 int MatCopy_Basic(Mat A,Mat B,MatStructure str) 2325 { 2326 int ierr,i,rstart,rend,nz,*cwork; 2327 PetscScalar *vwork; 2328 2329 PetscFunctionBegin; 2330 ierr = MatZeroEntries(B);CHKERRQ(ierr); 2331 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 2332 for (i=rstart; i<rend; i++) { 2333 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2334 ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr); 2335 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2336 } 2337 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2338 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2339 ierr = PetscObjectIncreaseState((PetscObject)B); CHKERRQ(ierr); 2340 PetscFunctionReturn(0); 2341 } 2342 2343 #undef __FUNCT__ 2344 #define __FUNCT__ "MatCopy" 2345 /*@C 2346 MatCopy - Copys a matrix to another matrix. 2347 2348 Collective on Mat 2349 2350 Input Parameters: 2351 + A - the matrix 2352 - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN 2353 2354 Output Parameter: 2355 . B - where the copy is put 2356 2357 Notes: 2358 If you use SAME_NONZERO_PATTERN then the two matrices had better have the 2359 same nonzero pattern or the routine will crash. 2360 2361 MatCopy() copies the matrix entries of a matrix to another existing 2362 matrix (after first zeroing the second matrix). A related routine is 2363 MatConvert(), which first creates a new matrix and then copies the data. 2364 2365 Level: intermediate 2366 2367 Concepts: matrices^copying 2368 2369 .seealso: MatConvert(), MatDuplicate() 2370 2371 @*/ 2372 int MatCopy(Mat A,Mat B,MatStructure str) 2373 { 2374 int ierr; 2375 2376 PetscFunctionBegin; 2377 PetscValidHeaderSpecific(A,MAT_COOKIE); 2378 PetscValidHeaderSpecific(B,MAT_COOKIE); 2379 PetscValidType(A); 2380 MatPreallocated(A); 2381 PetscValidType(B); 2382 MatPreallocated(B); 2383 PetscCheckSameComm(A,B); 2384 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2385 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2386 if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim (%d,%d) (%d,%d)",A->M,B->M, 2387 A->N,B->N); 2388 2389 ierr = PetscLogEventBegin(MAT_Copy,A,B,0,0);CHKERRQ(ierr); 2390 if (A->ops->copy) { 2391 ierr = (*A->ops->copy)(A,B,str);CHKERRQ(ierr); 2392 } else { /* generic conversion */ 2393 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2394 } 2395 if (A->mapping) { 2396 if (B->mapping) {ierr = ISLocalToGlobalMappingDestroy(B->mapping);CHKERRQ(ierr);B->mapping = 0;} 2397 ierr = MatSetLocalToGlobalMapping(B,A->mapping);CHKERRQ(ierr); 2398 } 2399 if (A->bmapping) { 2400 if (B->bmapping) {ierr = ISLocalToGlobalMappingDestroy(B->bmapping);CHKERRQ(ierr);B->bmapping = 0;} 2401 ierr = MatSetLocalToGlobalMappingBlock(B,A->mapping);CHKERRQ(ierr); 2402 } 2403 ierr = PetscLogEventEnd(MAT_Copy,A,B,0,0);CHKERRQ(ierr); 2404 ierr = PetscObjectIncreaseState((PetscObject)B); CHKERRQ(ierr); 2405 PetscFunctionReturn(0); 2406 } 2407 2408 #include "petscsys.h" 2409 PetscTruth MatConvertRegisterAllCalled = PETSC_FALSE; 2410 PetscFList MatConvertList = 0; 2411 2412 #undef __FUNCT__ 2413 #define __FUNCT__ "MatConvertRegister" 2414 /*@C 2415 MatConvertRegister - Allows one to register a routine that converts a sparse matrix 2416 from one format to another. 2417 2418 Not Collective 2419 2420 Input Parameters: 2421 + type - the type of matrix (defined in include/petscmat.h), for example, MATSEQAIJ. 2422 - Converter - the function that reads the matrix from the binary file. 2423 2424 Level: developer 2425 2426 .seealso: MatConvertRegisterAll(), MatConvert() 2427 2428 @*/ 2429 int MatConvertRegister(char *sname,char *path,char *name,int (*function)(Mat,MatType,Mat*)) 2430 { 2431 int ierr; 2432 char fullname[PETSC_MAX_PATH_LEN]; 2433 2434 PetscFunctionBegin; 2435 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 2436 ierr = PetscFListAdd(&MatConvertList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 2437 PetscFunctionReturn(0); 2438 } 2439 2440 #undef __FUNCT__ 2441 #define __FUNCT__ "MatConvert" 2442 /*@C 2443 MatConvert - Converts a matrix to another matrix, either of the same 2444 or different type. 2445 2446 Collective on Mat 2447 2448 Input Parameters: 2449 + mat - the matrix 2450 - newtype - new matrix type. Use MATSAME to create a new matrix of the 2451 same type as the original matrix. 2452 2453 Output Parameter: 2454 . M - pointer to place new matrix 2455 2456 Notes: 2457 MatConvert() first creates a new matrix and then copies the data from 2458 the first matrix. A related routine is MatCopy(), which copies the matrix 2459 entries of one matrix to another already existing matrix context. 2460 2461 Level: intermediate 2462 2463 Concepts: matrices^converting between storage formats 2464 2465 .seealso: MatCopy(), MatDuplicate() 2466 @*/ 2467 int MatConvert(Mat mat,MatType newtype,Mat *M) 2468 { 2469 int ierr; 2470 PetscTruth sametype,issame,flg; 2471 char convname[256],mtype[256]; 2472 2473 PetscFunctionBegin; 2474 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2475 PetscValidType(mat); 2476 MatPreallocated(mat); 2477 PetscValidPointer(M); 2478 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2479 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2480 2481 ierr = PetscOptionsGetString(PETSC_NULL,"-matconvert_type",mtype,256,&flg);CHKERRQ(ierr); 2482 if (flg) { 2483 newtype = mtype; 2484 } 2485 ierr = PetscLogEventBegin(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2486 2487 ierr = PetscTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr); 2488 ierr = PetscStrcmp(newtype,"same",&issame);CHKERRQ(ierr); 2489 if ((sametype || issame) && mat->ops->duplicate) { 2490 ierr = (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);CHKERRQ(ierr); 2491 } else { 2492 int (*conv)(Mat,MatType,Mat*)=PETSC_NULL; 2493 /* 2494 Order of precedence: 2495 1) See if a specialized converter is known to the current matrix. 2496 2) See if a specialized converter is known to the desired matrix class. 2497 3) See if a good general converter is registered for the desired class 2498 (as of 6/27/03 only MATMPIADJ falls into this category). 2499 4) See if a good general converter is known for the current matrix. 2500 5) Use a really basic converter. 2501 */ 2502 ierr = PetscStrcpy(convname,"MatConvert_");CHKERRQ(ierr); 2503 ierr = PetscStrcat(convname,mat->type_name);CHKERRQ(ierr); 2504 ierr = PetscStrcat(convname,"_");CHKERRQ(ierr); 2505 ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr); 2506 ierr = PetscStrcat(convname,"_C");CHKERRQ(ierr); 2507 ierr = PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2508 if (!conv) { 2509 Mat B; 2510 ierr = MatCreate(mat->comm,0,0,0,0,&B);CHKERRQ(ierr); 2511 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 2512 ierr = PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2513 ierr = MatDestroy(B);CHKERRQ(ierr); 2514 if (!conv) { 2515 if (!MatConvertRegisterAllCalled) { 2516 ierr = MatConvertRegisterAll(PETSC_NULL);CHKERRQ(ierr); 2517 } 2518 ierr = PetscFListFind(mat->comm,MatConvertList,newtype,(void(**)(void))&conv);CHKERRQ(ierr); 2519 if (!conv) { 2520 if (mat->ops->convert) { 2521 conv = mat->ops->convert; 2522 } else { 2523 conv = MatConvert_Basic; 2524 } 2525 } 2526 } 2527 } 2528 ierr = (*conv)(mat,newtype,M);CHKERRQ(ierr); 2529 } 2530 ierr = PetscLogEventEnd(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2531 ierr = PetscObjectIncreaseState((PetscObject)*M); CHKERRQ(ierr); 2532 PetscFunctionReturn(0); 2533 } 2534 2535 2536 #undef __FUNCT__ 2537 #define __FUNCT__ "MatDuplicate" 2538 /*@C 2539 MatDuplicate - Duplicates a matrix including the non-zero structure. 2540 2541 Collective on Mat 2542 2543 Input Parameters: 2544 + mat - the matrix 2545 - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero 2546 values as well or not 2547 2548 Output Parameter: 2549 . M - pointer to place new matrix 2550 2551 Level: intermediate 2552 2553 Concepts: matrices^duplicating 2554 2555 .seealso: MatCopy(), MatConvert() 2556 @*/ 2557 int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M) 2558 { 2559 int ierr; 2560 2561 PetscFunctionBegin; 2562 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2563 PetscValidType(mat); 2564 MatPreallocated(mat); 2565 PetscValidPointer(M); 2566 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2567 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2568 2569 *M = 0; 2570 ierr = PetscLogEventBegin(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2571 if (!mat->ops->duplicate) { 2572 SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type"); 2573 } 2574 ierr = (*mat->ops->duplicate)(mat,op,M);CHKERRQ(ierr); 2575 if (mat->mapping) { 2576 ierr = MatSetLocalToGlobalMapping(*M,mat->mapping);CHKERRQ(ierr); 2577 } 2578 if (mat->bmapping) { 2579 ierr = MatSetLocalToGlobalMappingBlock(*M,mat->mapping);CHKERRQ(ierr); 2580 } 2581 ierr = PetscLogEventEnd(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2582 ierr = PetscObjectIncreaseState((PetscObject)*M); CHKERRQ(ierr); 2583 PetscFunctionReturn(0); 2584 } 2585 2586 #undef __FUNCT__ 2587 #define __FUNCT__ "MatGetDiagonal" 2588 /*@ 2589 MatGetDiagonal - Gets the diagonal of a matrix. 2590 2591 Collective on Mat and Vec 2592 2593 Input Parameters: 2594 + mat - the matrix 2595 - v - the vector for storing the diagonal 2596 2597 Output Parameter: 2598 . v - the diagonal of the matrix 2599 2600 Notes: 2601 For the SeqAIJ matrix format, this routine may also be called 2602 on a LU factored matrix; in that case it routines the reciprocal of 2603 the diagonal entries in U. It returns the entries permuted by the 2604 row and column permutation used during the symbolic factorization. 2605 2606 Level: intermediate 2607 2608 Concepts: matrices^accessing diagonals 2609 2610 .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix(), MatGetRowMax() 2611 @*/ 2612 int MatGetDiagonal(Mat mat,Vec v) 2613 { 2614 int ierr; 2615 2616 PetscFunctionBegin; 2617 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2618 PetscValidType(mat); 2619 MatPreallocated(mat); 2620 PetscValidHeaderSpecific(v,VEC_COOKIE); 2621 /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */ 2622 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2623 if (!mat->ops->getdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2624 2625 ierr = (*mat->ops->getdiagonal)(mat,v);CHKERRQ(ierr); 2626 PetscFunctionReturn(0); 2627 } 2628 2629 #undef __FUNCT__ 2630 #define __FUNCT__ "MatGetRowMax" 2631 /*@ 2632 MatGetRowMax - Gets the maximum value (in absolute value) of each 2633 row of the matrix 2634 2635 Collective on Mat and Vec 2636 2637 Input Parameters: 2638 . mat - the matrix 2639 2640 Output Parameter: 2641 . v - the vector for storing the maximums 2642 2643 Level: intermediate 2644 2645 Concepts: matrices^getting row maximums 2646 2647 .seealso: MatGetDiagonal(), MatGetSubmatrices(), MatGetSubmatrix() 2648 @*/ 2649 int MatGetRowMax(Mat mat,Vec v) 2650 { 2651 int ierr; 2652 2653 PetscFunctionBegin; 2654 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2655 PetscValidType(mat); 2656 MatPreallocated(mat); 2657 PetscValidHeaderSpecific(v,VEC_COOKIE); 2658 /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */ 2659 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2660 if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2661 2662 ierr = (*mat->ops->getrowmax)(mat,v);CHKERRQ(ierr); 2663 PetscFunctionReturn(0); 2664 } 2665 2666 #undef __FUNCT__ 2667 #define __FUNCT__ "MatTranspose" 2668 /*@C 2669 MatTranspose - Computes an in-place or out-of-place transpose of a matrix. 2670 2671 Collective on Mat 2672 2673 Input Parameter: 2674 . mat - the matrix to transpose 2675 2676 Output Parameters: 2677 . B - the transpose (or pass in PETSC_NULL for an in-place transpose) 2678 2679 Level: intermediate 2680 2681 Concepts: matrices^transposing 2682 2683 .seealso: MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose() 2684 @*/ 2685 int MatTranspose(Mat mat,Mat *B) 2686 { 2687 int ierr; 2688 2689 PetscFunctionBegin; 2690 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2691 PetscValidType(mat); 2692 MatPreallocated(mat); 2693 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2694 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2695 if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2696 2697 ierr = PetscLogEventBegin(MAT_Transpose,mat,0,0,0);CHKERRQ(ierr); 2698 ierr = (*mat->ops->transpose)(mat,B);CHKERRQ(ierr); 2699 ierr = PetscLogEventEnd(MAT_Transpose,mat,0,0,0);CHKERRQ(ierr); 2700 if (B) {ierr = PetscObjectIncreaseState((PetscObject)*B); CHKERRQ(ierr);} 2701 PetscFunctionReturn(0); 2702 } 2703 2704 #undef __FUNCT__ 2705 #define __FUNCT__ "MatIsTranspose" 2706 /*@C 2707 MatIsTranspose - Test whether a matrix is another one's transpose, 2708 or its own, in which case it tests symmetry. 2709 2710 Collective on Mat 2711 2712 Input Parameter: 2713 + A - the matrix to test 2714 - B - the matrix to test against, this can equal the first parameter 2715 2716 Output Parameters: 2717 . flg - the result 2718 2719 Notes: 2720 Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm 2721 has a running time of the order of the number of nonzeros; the parallel 2722 test involves parallel copies of the block-offdiagonal parts of the matrix. 2723 2724 Level: intermediate 2725 2726 Concepts: matrices^transposing, matrix^symmetry 2727 2728 .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian() 2729 @*/ 2730 int MatIsTranspose(Mat A,Mat B,PetscTruth *flg) 2731 { 2732 int ierr,(*f)(Mat,Mat,PetscTruth*),(*g)(Mat,Mat,PetscTruth*); 2733 2734 PetscFunctionBegin; 2735 PetscValidHeaderSpecific(A,MAT_COOKIE); 2736 PetscValidHeaderSpecific(B,MAT_COOKIE); 2737 ierr = PetscObjectQueryFunction((PetscObject)A,"MatIsTranspose_C",(void (**)(void))&f);CHKERRQ(ierr); 2738 ierr = PetscObjectQueryFunction((PetscObject)B,"MatIsTranspose_C",(void (**)(void))&g);CHKERRQ(ierr); 2739 if (f && g) { 2740 if (f==g) { 2741 ierr = (*f)(A,B,flg);CHKERRQ(ierr); 2742 } else { 2743 SETERRQ(1,"Matrices do not have the same comparator for symmetry test"); 2744 } 2745 } 2746 PetscFunctionReturn(0); 2747 } 2748 2749 #undef __FUNCT__ 2750 #define __FUNCT__ "MatPermute" 2751 /*@C 2752 MatPermute - Creates a new matrix with rows and columns permuted from the 2753 original. 2754 2755 Collective on Mat 2756 2757 Input Parameters: 2758 + mat - the matrix to permute 2759 . row - row permutation, each processor supplies only the permutation for its rows 2760 - col - column permutation, each processor needs the entire column permutation, that is 2761 this is the same size as the total number of columns in the matrix 2762 2763 Output Parameters: 2764 . B - the permuted matrix 2765 2766 Level: advanced 2767 2768 Concepts: matrices^permuting 2769 2770 .seealso: MatGetOrdering() 2771 @*/ 2772 int MatPermute(Mat mat,IS row,IS col,Mat *B) 2773 { 2774 int ierr; 2775 2776 PetscFunctionBegin; 2777 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2778 PetscValidType(mat); 2779 MatPreallocated(mat); 2780 PetscValidHeaderSpecific(row,IS_COOKIE); 2781 PetscValidHeaderSpecific(col,IS_COOKIE); 2782 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2783 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2784 if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2785 ierr = (*mat->ops->permute)(mat,row,col,B);CHKERRQ(ierr); 2786 ierr = PetscObjectIncreaseState((PetscObject)*B); CHKERRQ(ierr); 2787 PetscFunctionReturn(0); 2788 } 2789 2790 #undef __FUNCT__ 2791 #define __FUNCT__ "MatPermuteSparsify" 2792 /*@C 2793 MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the 2794 original and sparsified to the prescribed tolerance. 2795 2796 Collective on Mat 2797 2798 Input Parameters: 2799 + A - The matrix to permute 2800 . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE 2801 . frac - The half-bandwidth as a fraction of the total size, or 0.0 2802 . tol - The drop tolerance 2803 . rowp - The row permutation 2804 - colp - The column permutation 2805 2806 Output Parameter: 2807 . B - The permuted, sparsified matrix 2808 2809 Level: advanced 2810 2811 Note: 2812 The default behavior (band = PETSC_DECIDE and frac = 0.0) is to 2813 restrict the half-bandwidth of the resulting matrix to 5% of the 2814 total matrix size. 2815 2816 .keywords: matrix, permute, sparsify 2817 2818 .seealso: MatGetOrdering(), MatPermute() 2819 @*/ 2820 int MatPermuteSparsify(Mat A, int band, PetscReal frac, PetscReal tol, IS rowp, IS colp, Mat *B) 2821 { 2822 IS irowp, icolp; 2823 int *rows, *cols; 2824 int M, N, locRowStart, locRowEnd; 2825 int nz, newNz; 2826 int *cwork, *cnew; 2827 PetscScalar *vwork, *vnew; 2828 int bw, size; 2829 int row, locRow, newRow, col, newCol; 2830 int ierr; 2831 2832 PetscFunctionBegin; 2833 PetscValidHeaderSpecific(A, MAT_COOKIE); 2834 PetscValidHeaderSpecific(rowp, IS_COOKIE); 2835 PetscValidHeaderSpecific(colp, IS_COOKIE); 2836 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 2837 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 2838 if (!A->ops->permutesparsify) { 2839 ierr = MatGetSize(A, &M, &N); CHKERRQ(ierr); 2840 ierr = MatGetOwnershipRange(A, &locRowStart, &locRowEnd); CHKERRQ(ierr); 2841 ierr = ISGetSize(rowp, &size); CHKERRQ(ierr); 2842 if (size != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for row permutation, should be %d", size, M); 2843 ierr = ISGetSize(colp, &size); CHKERRQ(ierr); 2844 if (size != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for column permutation, should be %d", size, N); 2845 ierr = ISInvertPermutation(rowp, 0, &irowp); CHKERRQ(ierr); 2846 ierr = ISGetIndices(irowp, &rows); CHKERRQ(ierr); 2847 ierr = ISInvertPermutation(colp, 0, &icolp); CHKERRQ(ierr); 2848 ierr = ISGetIndices(icolp, &cols); CHKERRQ(ierr); 2849 ierr = PetscMalloc(N * sizeof(int), &cnew); CHKERRQ(ierr); 2850 ierr = PetscMalloc(N * sizeof(PetscScalar), &vnew); CHKERRQ(ierr); 2851 2852 /* Setup bandwidth to include */ 2853 if (band == PETSC_DECIDE) { 2854 if (frac <= 0.0) 2855 bw = (int) (M * 0.05); 2856 else 2857 bw = (int) (M * frac); 2858 } else { 2859 if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer"); 2860 bw = band; 2861 } 2862 2863 /* Put values into new matrix */ 2864 ierr = MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B); CHKERRQ(ierr); 2865 for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) { 2866 ierr = MatGetRow(A, row, &nz, &cwork, &vwork); CHKERRQ(ierr); 2867 newRow = rows[locRow]+locRowStart; 2868 for(col = 0, newNz = 0; col < nz; col++) { 2869 newCol = cols[cwork[col]]; 2870 if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) { 2871 cnew[newNz] = newCol; 2872 vnew[newNz] = vwork[col]; 2873 newNz++; 2874 } 2875 } 2876 ierr = MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES); CHKERRQ(ierr); 2877 ierr = MatRestoreRow(A, row, &nz, &cwork, &vwork); CHKERRQ(ierr); 2878 } 2879 ierr = PetscFree(cnew); CHKERRQ(ierr); 2880 ierr = PetscFree(vnew); CHKERRQ(ierr); 2881 ierr = MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2882 ierr = MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2883 ierr = ISRestoreIndices(irowp, &rows); CHKERRQ(ierr); 2884 ierr = ISRestoreIndices(icolp, &cols); CHKERRQ(ierr); 2885 ierr = ISDestroy(irowp); CHKERRQ(ierr); 2886 ierr = ISDestroy(icolp); CHKERRQ(ierr); 2887 } else { 2888 ierr = (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B); CHKERRQ(ierr); 2889 } 2890 ierr = PetscObjectIncreaseState((PetscObject)*B); CHKERRQ(ierr); 2891 PetscFunctionReturn(0); 2892 } 2893 2894 #undef __FUNCT__ 2895 #define __FUNCT__ "MatEqual" 2896 /*@ 2897 MatEqual - Compares two matrices. 2898 2899 Collective on Mat 2900 2901 Input Parameters: 2902 + A - the first matrix 2903 - B - the second matrix 2904 2905 Output Parameter: 2906 . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise. 2907 2908 Level: intermediate 2909 2910 Concepts: matrices^equality between 2911 @*/ 2912 int MatEqual(Mat A,Mat B,PetscTruth *flg) 2913 { 2914 int ierr; 2915 2916 PetscFunctionBegin; 2917 PetscValidHeaderSpecific(A,MAT_COOKIE); 2918 PetscValidHeaderSpecific(B,MAT_COOKIE); 2919 PetscValidType(A); 2920 MatPreallocated(A); 2921 PetscValidType(B); 2922 MatPreallocated(B); 2923 PetscValidIntPointer(flg); 2924 PetscCheckSameComm(A,B); 2925 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2926 if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2927 if (A->M != B->M || A->N != B->N) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: global dim %d %d %d %d",A->M,B->M,A->N,B->N); 2928 if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",A->type_name); 2929 if (!B->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",B->type_name); 2930 if (A->ops->equal != B->ops->equal) SETERRQ2(PETSC_ERR_ARG_INCOMP,"A is type: %s\nB is type: %s",A->type_name,B->type_name); 2931 ierr = (*A->ops->equal)(A,B,flg);CHKERRQ(ierr); 2932 PetscFunctionReturn(0); 2933 } 2934 2935 #undef __FUNCT__ 2936 #define __FUNCT__ "MatDiagonalScale" 2937 /*@ 2938 MatDiagonalScale - Scales a matrix on the left and right by diagonal 2939 matrices that are stored as vectors. Either of the two scaling 2940 matrices can be PETSC_NULL. 2941 2942 Collective on Mat 2943 2944 Input Parameters: 2945 + mat - the matrix to be scaled 2946 . l - the left scaling vector (or PETSC_NULL) 2947 - r - the right scaling vector (or PETSC_NULL) 2948 2949 Notes: 2950 MatDiagonalScale() computes A = LAR, where 2951 L = a diagonal matrix, R = a diagonal matrix 2952 2953 Level: intermediate 2954 2955 Concepts: matrices^diagonal scaling 2956 Concepts: diagonal scaling of matrices 2957 2958 .seealso: MatScale() 2959 @*/ 2960 int MatDiagonalScale(Mat mat,Vec l,Vec r) 2961 { 2962 int ierr; 2963 2964 PetscFunctionBegin; 2965 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2966 PetscValidType(mat); 2967 MatPreallocated(mat); 2968 if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2969 if (l) {PetscValidHeaderSpecific(l,VEC_COOKIE);PetscCheckSameComm(mat,l);} 2970 if (r) {PetscValidHeaderSpecific(r,VEC_COOKIE);PetscCheckSameComm(mat,r);} 2971 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2972 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2973 2974 ierr = PetscLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 2975 ierr = (*mat->ops->diagonalscale)(mat,l,r);CHKERRQ(ierr); 2976 ierr = PetscLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 2977 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 2978 PetscFunctionReturn(0); 2979 } 2980 2981 #undef __FUNCT__ 2982 #define __FUNCT__ "MatScale" 2983 /*@ 2984 MatScale - Scales all elements of a matrix by a given number. 2985 2986 Collective on Mat 2987 2988 Input Parameters: 2989 + mat - the matrix to be scaled 2990 - a - the scaling value 2991 2992 Output Parameter: 2993 . mat - the scaled matrix 2994 2995 Level: intermediate 2996 2997 Concepts: matrices^scaling all entries 2998 2999 .seealso: MatDiagonalScale() 3000 @*/ 3001 int MatScale(const PetscScalar *a,Mat mat) 3002 { 3003 int ierr; 3004 3005 PetscFunctionBegin; 3006 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3007 PetscValidType(mat); 3008 MatPreallocated(mat); 3009 PetscValidScalarPointer(a); 3010 if (!mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3011 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3012 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3013 3014 ierr = PetscLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 3015 ierr = (*mat->ops->scale)(a,mat);CHKERRQ(ierr); 3016 ierr = PetscLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 3017 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3018 PetscFunctionReturn(0); 3019 } 3020 3021 #undef __FUNCT__ 3022 #define __FUNCT__ "MatNorm" 3023 /*@ 3024 MatNorm - Calculates various norms of a matrix. 3025 3026 Collective on Mat 3027 3028 Input Parameters: 3029 + mat - the matrix 3030 - type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY 3031 3032 Output Parameters: 3033 . nrm - the resulting norm 3034 3035 Level: intermediate 3036 3037 Concepts: matrices^norm 3038 Concepts: norm^of matrix 3039 @*/ 3040 int MatNorm(Mat mat,NormType type,PetscReal *nrm) 3041 { 3042 int ierr; 3043 3044 PetscFunctionBegin; 3045 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3046 PetscValidType(mat); 3047 MatPreallocated(mat); 3048 PetscValidScalarPointer(nrm); 3049 3050 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3051 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3052 if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3053 ierr = (*mat->ops->norm)(mat,type,nrm);CHKERRQ(ierr); 3054 PetscFunctionReturn(0); 3055 } 3056 3057 /* 3058 This variable is used to prevent counting of MatAssemblyBegin() that 3059 are called from within a MatAssemblyEnd(). 3060 */ 3061 static int MatAssemblyEnd_InUse = 0; 3062 #undef __FUNCT__ 3063 #define __FUNCT__ "MatAssemblyBegin" 3064 /*@ 3065 MatAssemblyBegin - Begins assembling the matrix. This routine should 3066 be called after completing all calls to MatSetValues(). 3067 3068 Collective on Mat 3069 3070 Input Parameters: 3071 + mat - the matrix 3072 - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 3073 3074 Notes: 3075 MatSetValues() generally caches the values. The matrix is ready to 3076 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 3077 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 3078 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 3079 using the matrix. 3080 3081 Level: beginner 3082 3083 Concepts: matrices^assembling 3084 3085 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled() 3086 @*/ 3087 int MatAssemblyBegin(Mat mat,MatAssemblyType type) 3088 { 3089 int ierr; 3090 3091 PetscFunctionBegin; 3092 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3093 PetscValidType(mat); 3094 MatPreallocated(mat); 3095 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?"); 3096 if (mat->assembled) { 3097 mat->was_assembled = PETSC_TRUE; 3098 mat->assembled = PETSC_FALSE; 3099 } 3100 if (!MatAssemblyEnd_InUse) { 3101 ierr = PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);CHKERRQ(ierr); 3102 if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);} 3103 ierr = PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);CHKERRQ(ierr); 3104 } else { 3105 if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);} 3106 } 3107 PetscFunctionReturn(0); 3108 } 3109 3110 #undef __FUNCT__ 3111 #define __FUNCT__ "MatAssembed" 3112 /*@ 3113 MatAssembled - Indicates if a matrix has been assembled and is ready for 3114 use; for example, in matrix-vector product. 3115 3116 Collective on Mat 3117 3118 Input Parameter: 3119 . mat - the matrix 3120 3121 Output Parameter: 3122 . assembled - PETSC_TRUE or PETSC_FALSE 3123 3124 Level: advanced 3125 3126 Concepts: matrices^assembled? 3127 3128 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin() 3129 @*/ 3130 int MatAssembled(Mat mat,PetscTruth *assembled) 3131 { 3132 PetscFunctionBegin; 3133 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3134 PetscValidType(mat); 3135 MatPreallocated(mat); 3136 *assembled = mat->assembled; 3137 PetscFunctionReturn(0); 3138 } 3139 3140 #undef __FUNCT__ 3141 #define __FUNCT__ "MatView_Private" 3142 /* 3143 Processes command line options to determine if/how a matrix 3144 is to be viewed. Called by MatAssemblyEnd() and MatLoad(). 3145 */ 3146 int MatView_Private(Mat mat) 3147 { 3148 int ierr; 3149 PetscTruth flg; 3150 static PetscTruth incall = PETSC_FALSE; 3151 3152 PetscFunctionBegin; 3153 if (incall) PetscFunctionReturn(0); 3154 incall = PETSC_TRUE; 3155 ierr = PetscOptionsBegin(mat->comm,mat->prefix,"Matrix Options","Mat");CHKERRQ(ierr); 3156 ierr = PetscOptionsName("-mat_view_info","Information on matrix size","MatView",&flg);CHKERRQ(ierr); 3157 if (flg) { 3158 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 3159 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3160 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3161 } 3162 ierr = PetscOptionsName("-mat_view_info_detailed","Nonzeros in the matrix","MatView",&flg);CHKERRQ(ierr); 3163 if (flg) { 3164 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO_DETAIL);CHKERRQ(ierr); 3165 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3166 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3167 } 3168 ierr = PetscOptionsName("-mat_view","Print matrix to stdout","MatView",&flg);CHKERRQ(ierr); 3169 if (flg) { 3170 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3171 } 3172 ierr = PetscOptionsName("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",&flg);CHKERRQ(ierr); 3173 if (flg) { 3174 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); 3175 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3176 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3177 } 3178 ierr = PetscOptionsName("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",&flg);CHKERRQ(ierr); 3179 if (flg) { 3180 ierr = MatView(mat,PETSC_VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr); 3181 ierr = PetscViewerFlush(PETSC_VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr); 3182 } 3183 ierr = PetscOptionsName("-mat_view_binary","Save matrix to file in binary format","MatView",&flg);CHKERRQ(ierr); 3184 if (flg) { 3185 ierr = MatView(mat,PETSC_VIEWER_BINARY_(mat->comm));CHKERRQ(ierr); 3186 ierr = PetscViewerFlush(PETSC_VIEWER_BINARY_(mat->comm));CHKERRQ(ierr); 3187 } 3188 ierr = PetscOptionsEnd();CHKERRQ(ierr); 3189 /* cannot have inside PetscOptionsBegin() because uses PetscOptionsBegin() */ 3190 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_draw",&flg);CHKERRQ(ierr); 3191 if (flg) { 3192 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_contour",&flg);CHKERRQ(ierr); 3193 if (flg) { 3194 PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 3195 } 3196 ierr = MatView(mat,PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3197 ierr = PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3198 if (flg) { 3199 PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3200 } 3201 } 3202 incall = PETSC_FALSE; 3203 PetscFunctionReturn(0); 3204 } 3205 3206 #undef __FUNCT__ 3207 #define __FUNCT__ "MatAssemblyEnd" 3208 /*@ 3209 MatAssemblyEnd - Completes assembling the matrix. This routine should 3210 be called after MatAssemblyBegin(). 3211 3212 Collective on Mat 3213 3214 Input Parameters: 3215 + mat - the matrix 3216 - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 3217 3218 Options Database Keys: 3219 + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly() 3220 . -mat_view_info_detailed - Prints more detailed info 3221 . -mat_view - Prints matrix in ASCII format 3222 . -mat_view_matlab - Prints matrix in Matlab format 3223 . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX(). 3224 . -display <name> - Sets display name (default is host) 3225 . -draw_pause <sec> - Sets number of seconds to pause after display 3226 . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual) 3227 . -viewer_socket_machine <machine> 3228 . -viewer_socket_port <port> 3229 . -mat_view_binary - save matrix to file in binary format 3230 - -viewer_binary_filename <name> 3231 3232 Notes: 3233 MatSetValues() generally caches the values. The matrix is ready to 3234 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 3235 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 3236 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 3237 using the matrix. 3238 3239 Level: beginner 3240 3241 .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen() 3242 @*/ 3243 int MatAssemblyEnd(Mat mat,MatAssemblyType type) 3244 { 3245 int ierr; 3246 static int inassm = 0; 3247 PetscTruth flg; 3248 3249 PetscFunctionBegin; 3250 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3251 PetscValidType(mat); 3252 MatPreallocated(mat); 3253 3254 inassm++; 3255 MatAssemblyEnd_InUse++; 3256 if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */ 3257 ierr = PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);CHKERRQ(ierr); 3258 if (mat->ops->assemblyend) { 3259 ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr); 3260 } 3261 ierr = PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);CHKERRQ(ierr); 3262 } else { 3263 if (mat->ops->assemblyend) { 3264 ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr); 3265 } 3266 } 3267 3268 /* Flush assembly is not a true assembly */ 3269 if (type != MAT_FLUSH_ASSEMBLY) { 3270 mat->assembled = PETSC_TRUE; mat->num_ass++; 3271 } 3272 mat->insertmode = NOT_SET_VALUES; 3273 MatAssemblyEnd_InUse--; 3274 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3275 if (!mat->symmetric_eternal) { 3276 mat->symmetric_set = PETSC_FALSE; 3277 mat->hermitian_set = PETSC_FALSE; 3278 mat->structurally_symmetric_set = PETSC_FALSE; 3279 } 3280 if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) { 3281 ierr = MatView_Private(mat);CHKERRQ(ierr); 3282 } 3283 inassm--; 3284 ierr = PetscOptionsHasName(mat->prefix,"-help",&flg);CHKERRQ(ierr); 3285 if (flg) { 3286 ierr = MatPrintHelp(mat);CHKERRQ(ierr); 3287 } 3288 PetscFunctionReturn(0); 3289 } 3290 3291 3292 #undef __FUNCT__ 3293 #define __FUNCT__ "MatCompress" 3294 /*@ 3295 MatCompress - Tries to store the matrix in as little space as 3296 possible. May fail if memory is already fully used, since it 3297 tries to allocate new space. 3298 3299 Collective on Mat 3300 3301 Input Parameters: 3302 . mat - the matrix 3303 3304 Level: advanced 3305 3306 @*/ 3307 int MatCompress(Mat mat) 3308 { 3309 int ierr; 3310 3311 PetscFunctionBegin; 3312 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3313 PetscValidType(mat); 3314 MatPreallocated(mat); 3315 if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);} 3316 PetscFunctionReturn(0); 3317 } 3318 3319 #undef __FUNCT__ 3320 #define __FUNCT__ "MatSetOption" 3321 /*@ 3322 MatSetOption - Sets a parameter option for a matrix. Some options 3323 may be specific to certain storage formats. Some options 3324 determine how values will be inserted (or added). Sorted, 3325 row-oriented input will generally assemble the fastest. The default 3326 is row-oriented, nonsorted input. 3327 3328 Collective on Mat 3329 3330 Input Parameters: 3331 + mat - the matrix 3332 - option - the option, one of those listed below (and possibly others), 3333 e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR 3334 3335 Options Describing Matrix Structure: 3336 + MAT_SYMMETRIC - symmetric in terms of both structure and value 3337 . MAT_HERMITIAN - transpose is the complex conjugation 3338 . MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure 3339 . MAT_NOT_SYMMETRIC - not symmetric in value 3340 . MAT_NOT_HERMITIAN - transpose is not the complex conjugation 3341 . MAT_NOT_STRUCTURALLY_SYMMETRIC - not symmetric nonzero structure 3342 . MAT_SYMMETRY_ETERNAL - if you would like the symmetry/Hermitian flag 3343 you set to be kept with all future use of the matrix 3344 including after MatAssemblyBegin/End() which could 3345 potentially change the symmetry structure, i.e. you 3346 KNOW the matrix will ALWAYS have the property you set. 3347 - MAT_NOT_SYMMETRY_ETERNAL - if MatAssemblyBegin/End() is called then the 3348 flags you set will be dropped (in case potentially 3349 the symmetry etc was lost). 3350 3351 Options For Use with MatSetValues(): 3352 Insert a logically dense subblock, which can be 3353 + MAT_ROW_ORIENTED - row-oriented (default) 3354 . MAT_COLUMN_ORIENTED - column-oriented 3355 . MAT_ROWS_SORTED - sorted by row 3356 . MAT_ROWS_UNSORTED - not sorted by row (default) 3357 . MAT_COLUMNS_SORTED - sorted by column 3358 - MAT_COLUMNS_UNSORTED - not sorted by column (default) 3359 3360 Not these options reflect the data you pass in with MatSetValues(); it has 3361 nothing to do with how the data is stored internally in the matrix 3362 data structure. 3363 3364 When (re)assembling a matrix, we can restrict the input for 3365 efficiency/debugging purposes. These options include 3366 + MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be 3367 allowed if they generate a new nonzero 3368 . MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed 3369 . MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if 3370 they generate a nonzero in a new diagonal (for block diagonal format only) 3371 . MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only) 3372 . MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries 3373 . MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry 3374 - MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly 3375 3376 Notes: 3377 Some options are relevant only for particular matrix types and 3378 are thus ignored by others. Other options are not supported by 3379 certain matrix types and will generate an error message if set. 3380 3381 If using a Fortran 77 module to compute a matrix, one may need to 3382 use the column-oriented option (or convert to the row-oriented 3383 format). 3384 3385 MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion 3386 that would generate a new entry in the nonzero structure is instead 3387 ignored. Thus, if memory has not alredy been allocated for this particular 3388 data, then the insertion is ignored. For dense matrices, in which 3389 the entire array is allocated, no entries are ever ignored. 3390 Set after the first MatAssemblyEnd() 3391 3392 MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion 3393 that would generate a new entry in the nonzero structure instead produces 3394 an error. (Currently supported for AIJ and BAIJ formats only.) 3395 This is a useful flag when using SAME_NONZERO_PATTERN in calling 3396 SLESSetOperators() to ensure that the nonzero pattern truely does 3397 remain unchanged. Set after the first MatAssemblyEnd() 3398 3399 MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion 3400 that would generate a new entry that has not been preallocated will 3401 instead produce an error. (Currently supported for AIJ and BAIJ formats 3402 only.) This is a useful flag when debugging matrix memory preallocation. 3403 3404 MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for 3405 other processors should be dropped, rather than stashed. 3406 This is useful if you know that the "owning" processor is also 3407 always generating the correct matrix entries, so that PETSc need 3408 not transfer duplicate entries generated on another processor. 3409 3410 MAT_USE_HASH_TABLE indicates that a hash table be used to improve the 3411 searches during matrix assembly. When this flag is set, the hash table 3412 is created during the first Matrix Assembly. This hash table is 3413 used the next time through, during MatSetVaules()/MatSetVaulesBlocked() 3414 to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag 3415 should be used with MAT_USE_HASH_TABLE flag. This option is currently 3416 supported by MATMPIBAIJ format only. 3417 3418 MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries 3419 are kept in the nonzero structure 3420 3421 MAT_IGNORE_ZERO_ENTRIES - for AIJ matrices this will stop zero values from creating 3422 a zero location in the matrix 3423 3424 MAT_USE_INODES - indicates using inode version of the code - works with AIJ and 3425 ROWBS matrix types 3426 3427 MAT_DO_NOT_USE_INODES - indicates not using inode version of the code - works 3428 with AIJ and ROWBS matrix types 3429 3430 Level: intermediate 3431 3432 Concepts: matrices^setting options 3433 3434 @*/ 3435 int MatSetOption(Mat mat,MatOption op) 3436 { 3437 int ierr; 3438 3439 PetscFunctionBegin; 3440 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3441 PetscValidType(mat); 3442 MatPreallocated(mat); 3443 switch (op) { 3444 case MAT_SYMMETRIC: 3445 mat->symmetric = PETSC_TRUE; 3446 mat->structurally_symmetric = PETSC_TRUE; 3447 mat->symmetric_set = PETSC_TRUE; 3448 mat->structurally_symmetric_set = PETSC_TRUE; 3449 break; 3450 case MAT_HERMITIAN: 3451 mat->hermitian = PETSC_TRUE; 3452 mat->structurally_symmetric = PETSC_TRUE; 3453 mat->hermitian_set = PETSC_TRUE; 3454 mat->structurally_symmetric_set = PETSC_TRUE; 3455 break; 3456 case MAT_STRUCTURALLY_SYMMETRIC: 3457 mat->structurally_symmetric = PETSC_TRUE; 3458 mat->structurally_symmetric_set = PETSC_TRUE; 3459 break; 3460 case MAT_NOT_SYMMETRIC: 3461 mat->symmetric = PETSC_FALSE; 3462 mat->symmetric_set = PETSC_TRUE; 3463 break; 3464 case MAT_NOT_HERMITIAN: 3465 mat->hermitian = PETSC_FALSE; 3466 mat->hermitian_set = PETSC_TRUE; 3467 break; 3468 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 3469 mat->structurally_symmetric = PETSC_FALSE; 3470 mat->structurally_symmetric_set = PETSC_TRUE; 3471 break; 3472 case MAT_SYMMETRY_ETERNAL: 3473 mat->symmetric_eternal = PETSC_TRUE; 3474 case MAT_NOT_SYMMETRY_ETERNAL: 3475 mat->symmetric_eternal = PETSC_FALSE; 3476 default: 3477 break; 3478 } 3479 if (mat->ops->setoption) { 3480 ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr); 3481 } 3482 PetscFunctionReturn(0); 3483 } 3484 3485 #undef __FUNCT__ 3486 #define __FUNCT__ "MatZeroEntries" 3487 /*@ 3488 MatZeroEntries - Zeros all entries of a matrix. For sparse matrices 3489 this routine retains the old nonzero structure. 3490 3491 Collective on Mat 3492 3493 Input Parameters: 3494 . mat - the matrix 3495 3496 Level: intermediate 3497 3498 Concepts: matrices^zeroing 3499 3500 .seealso: MatZeroRows() 3501 @*/ 3502 int MatZeroEntries(Mat mat) 3503 { 3504 int ierr; 3505 3506 PetscFunctionBegin; 3507 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3508 PetscValidType(mat); 3509 MatPreallocated(mat); 3510 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3511 if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3512 3513 ierr = PetscLogEventBegin(MAT_ZeroEntries,mat,0,0,0);CHKERRQ(ierr); 3514 ierr = (*mat->ops->zeroentries)(mat);CHKERRQ(ierr); 3515 ierr = PetscLogEventEnd(MAT_ZeroEntries,mat,0,0,0);CHKERRQ(ierr); 3516 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3517 PetscFunctionReturn(0); 3518 } 3519 3520 #undef __FUNCT__ 3521 #define __FUNCT__ "MatZeroRows" 3522 /*@C 3523 MatZeroRows - Zeros all entries (except possibly the main diagonal) 3524 of a set of rows of a matrix. 3525 3526 Collective on Mat 3527 3528 Input Parameters: 3529 + mat - the matrix 3530 . is - index set of rows to remove 3531 - diag - pointer to value put in all diagonals of eliminated rows. 3532 Note that diag is not a pointer to an array, but merely a 3533 pointer to a single value. 3534 3535 Notes: 3536 For the AIJ and BAIJ matrix formats this removes the old nonzero structure, 3537 but does not release memory. For the dense and block diagonal 3538 formats this does not alter the nonzero structure. 3539 3540 If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure 3541 of the matrix is not changed (even for AIJ and BAIJ matrices) the values are 3542 merely zeroed. 3543 3544 The user can set a value in the diagonal entry (or for the AIJ and 3545 row formats can optionally remove the main diagonal entry from the 3546 nonzero structure as well, by passing a null pointer (PETSC_NULL 3547 in C or PETSC_NULL_SCALAR in Fortran) as the final argument). 3548 3549 For the parallel case, all processes that share the matrix (i.e., 3550 those in the communicator used for matrix creation) MUST call this 3551 routine, regardless of whether any rows being zeroed are owned by 3552 them. 3553 3554 Level: intermediate 3555 3556 Concepts: matrices^zeroing rows 3557 3558 .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption() 3559 @*/ 3560 int MatZeroRows(Mat mat,IS is,const PetscScalar *diag) 3561 { 3562 int ierr; 3563 3564 PetscFunctionBegin; 3565 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3566 PetscValidType(mat); 3567 MatPreallocated(mat); 3568 PetscValidHeaderSpecific(is,IS_COOKIE); 3569 if (diag) PetscValidScalarPointer(diag); 3570 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3571 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3572 if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3573 3574 ierr = (*mat->ops->zerorows)(mat,is,diag);CHKERRQ(ierr); 3575 ierr = MatView_Private(mat);CHKERRQ(ierr); 3576 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3577 PetscFunctionReturn(0); 3578 } 3579 3580 #undef __FUNCT__ 3581 #define __FUNCT__ "MatZeroRowsLocal" 3582 /*@C 3583 MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal) 3584 of a set of rows of a matrix; using local numbering of rows. 3585 3586 Collective on Mat 3587 3588 Input Parameters: 3589 + mat - the matrix 3590 . is - index set of rows to remove 3591 - diag - pointer to value put in all diagonals of eliminated rows. 3592 Note that diag is not a pointer to an array, but merely a 3593 pointer to a single value. 3594 3595 Notes: 3596 Before calling MatZeroRowsLocal(), the user must first set the 3597 local-to-global mapping by calling MatSetLocalToGlobalMapping(). 3598 3599 For the AIJ matrix formats this removes the old nonzero structure, 3600 but does not release memory. For the dense and block diagonal 3601 formats this does not alter the nonzero structure. 3602 3603 If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure 3604 of the matrix is not changed (even for AIJ and BAIJ matrices) the values are 3605 merely zeroed. 3606 3607 The user can set a value in the diagonal entry (or for the AIJ and 3608 row formats can optionally remove the main diagonal entry from the 3609 nonzero structure as well, by passing a null pointer (PETSC_NULL 3610 in C or PETSC_NULL_SCALAR in Fortran) as the final argument). 3611 3612 Level: intermediate 3613 3614 Concepts: matrices^zeroing 3615 3616 .seealso: MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping 3617 @*/ 3618 int MatZeroRowsLocal(Mat mat,IS is,const PetscScalar *diag) 3619 { 3620 int ierr; 3621 IS newis; 3622 3623 PetscFunctionBegin; 3624 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3625 PetscValidType(mat); 3626 MatPreallocated(mat); 3627 PetscValidHeaderSpecific(is,IS_COOKIE); 3628 if (diag) PetscValidScalarPointer(diag); 3629 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3630 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3631 3632 if (mat->ops->zerorowslocal) { 3633 ierr = (*mat->ops->zerorowslocal)(mat,is,diag);CHKERRQ(ierr); 3634 } else { 3635 if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first"); 3636 ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);CHKERRQ(ierr); 3637 ierr = (*mat->ops->zerorows)(mat,newis,diag);CHKERRQ(ierr); 3638 ierr = ISDestroy(newis);CHKERRQ(ierr); 3639 } 3640 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3641 PetscFunctionReturn(0); 3642 } 3643 3644 #undef __FUNCT__ 3645 #define __FUNCT__ "MatGetSize" 3646 /*@ 3647 MatGetSize - Returns the numbers of rows and columns in a matrix. 3648 3649 Not Collective 3650 3651 Input Parameter: 3652 . mat - the matrix 3653 3654 Output Parameters: 3655 + m - the number of global rows 3656 - n - the number of global columns 3657 3658 Level: beginner 3659 3660 Concepts: matrices^size 3661 3662 .seealso: MatGetLocalSize() 3663 @*/ 3664 int MatGetSize(Mat mat,int *m,int* n) 3665 { 3666 PetscFunctionBegin; 3667 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3668 if (m) *m = mat->M; 3669 if (n) *n = mat->N; 3670 PetscFunctionReturn(0); 3671 } 3672 3673 #undef __FUNCT__ 3674 #define __FUNCT__ "MatGetLocalSize" 3675 /*@ 3676 MatGetLocalSize - Returns the number of rows and columns in a matrix 3677 stored locally. This information may be implementation dependent, so 3678 use with care. 3679 3680 Not Collective 3681 3682 Input Parameters: 3683 . mat - the matrix 3684 3685 Output Parameters: 3686 + m - the number of local rows 3687 - n - the number of local columns 3688 3689 Level: beginner 3690 3691 Concepts: matrices^local size 3692 3693 .seealso: MatGetSize() 3694 @*/ 3695 int MatGetLocalSize(Mat mat,int *m,int* n) 3696 { 3697 PetscFunctionBegin; 3698 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3699 if (m) *m = mat->m; 3700 if (n) *n = mat->n; 3701 PetscFunctionReturn(0); 3702 } 3703 3704 #undef __FUNCT__ 3705 #define __FUNCT__ "MatGetOwnershipRange" 3706 /*@ 3707 MatGetOwnershipRange - Returns the range of matrix rows owned by 3708 this processor, assuming that the matrix is laid out with the first 3709 n1 rows on the first processor, the next n2 rows on the second, etc. 3710 For certain parallel layouts this range may not be well defined. 3711 3712 Not Collective 3713 3714 Input Parameters: 3715 . mat - the matrix 3716 3717 Output Parameters: 3718 + m - the global index of the first local row 3719 - n - one more than the global index of the last local row 3720 3721 Level: beginner 3722 3723 Concepts: matrices^row ownership 3724 @*/ 3725 int MatGetOwnershipRange(Mat mat,int *m,int* n) 3726 { 3727 int ierr; 3728 3729 PetscFunctionBegin; 3730 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3731 PetscValidType(mat); 3732 MatPreallocated(mat); 3733 if (m) PetscValidIntPointer(m); 3734 if (n) PetscValidIntPointer(n); 3735 ierr = PetscMapGetLocalRange(mat->rmap,m,n);CHKERRQ(ierr); 3736 PetscFunctionReturn(0); 3737 } 3738 3739 #undef __FUNCT__ 3740 #define __FUNCT__ "MatILUFactorSymbolic" 3741 /*@ 3742 MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix. 3743 Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric() 3744 to complete the factorization. 3745 3746 Collective on Mat 3747 3748 Input Parameters: 3749 + mat - the matrix 3750 . row - row permutation 3751 . column - column permutation 3752 - info - structure containing 3753 $ levels - number of levels of fill. 3754 $ expected fill - as ratio of original fill. 3755 $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices 3756 missing diagonal entries) 3757 3758 Output Parameters: 3759 . fact - new matrix that has been symbolically factored 3760 3761 Notes: 3762 See the users manual for additional information about 3763 choosing the fill factor for better efficiency. 3764 3765 Most users should employ the simplified SLES interface for linear solvers 3766 instead of working directly with matrix algebra routines such as this. 3767 See, e.g., SLESCreate(). 3768 3769 Level: developer 3770 3771 Concepts: matrices^symbolic LU factorization 3772 Concepts: matrices^factorization 3773 Concepts: LU^symbolic factorization 3774 3775 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 3776 MatGetOrdering(), MatFactorInfo 3777 3778 @*/ 3779 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact) 3780 { 3781 int ierr; 3782 3783 PetscFunctionBegin; 3784 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3785 PetscValidType(mat); 3786 MatPreallocated(mat); 3787 PetscValidPointer(fact); 3788 PetscValidHeaderSpecific(row,IS_COOKIE); 3789 PetscValidHeaderSpecific(col,IS_COOKIE); 3790 if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %d",(int)info->levels); 3791 if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill); 3792 if (!mat->ops->ilufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ILU",mat->type_name); 3793 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3794 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3795 3796 ierr = PetscLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 3797 ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr); 3798 ierr = PetscLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 3799 PetscFunctionReturn(0); 3800 } 3801 3802 #undef __FUNCT__ 3803 #define __FUNCT__ "MatICCFactorSymbolic" 3804 /*@ 3805 MatICCFactorSymbolic - Performs symbolic incomplete 3806 Cholesky factorization for a symmetric matrix. Use 3807 MatCholeskyFactorNumeric() to complete the factorization. 3808 3809 Collective on Mat 3810 3811 Input Parameters: 3812 + mat - the matrix 3813 . perm - row and column permutation 3814 - info - structure containing 3815 $ levels - number of levels of fill. 3816 $ expected fill - as ratio of original fill. 3817 3818 Output Parameter: 3819 . fact - the factored matrix 3820 3821 Notes: 3822 Currently only no-fill factorization is supported. 3823 3824 Most users should employ the simplified SLES interface for linear solvers 3825 instead of working directly with matrix algebra routines such as this. 3826 See, e.g., SLESCreate(). 3827 3828 Level: developer 3829 3830 Concepts: matrices^symbolic incomplete Cholesky factorization 3831 Concepts: matrices^factorization 3832 Concepts: Cholsky^symbolic factorization 3833 3834 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor(), MatFactorInfo 3835 @*/ 3836 int MatICCFactorSymbolic(Mat mat,IS perm,MatFactorInfo *info,Mat *fact) 3837 { 3838 int ierr; 3839 3840 PetscFunctionBegin; 3841 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3842 PetscValidType(mat); 3843 MatPreallocated(mat); 3844 PetscValidPointer(fact); 3845 PetscValidHeaderSpecific(perm,IS_COOKIE); 3846 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3847 if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels negative %d",(int) info->levels); 3848 if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill); 3849 if (!mat->ops->iccfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ICC",mat->type_name); 3850 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3851 3852 ierr = PetscLogEventBegin(MAT_ICCFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 3853 ierr = (*mat->ops->iccfactorsymbolic)(mat,perm,info,fact);CHKERRQ(ierr); 3854 ierr = PetscLogEventEnd(MAT_ICCFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 3855 PetscFunctionReturn(0); 3856 } 3857 3858 #undef __FUNCT__ 3859 #define __FUNCT__ "MatGetArray" 3860 /*@C 3861 MatGetArray - Returns a pointer to the element values in the matrix. 3862 The result of this routine is dependent on the underlying matrix data 3863 structure, and may not even work for certain matrix types. You MUST 3864 call MatRestoreArray() when you no longer need to access the array. 3865 3866 Not Collective 3867 3868 Input Parameter: 3869 . mat - the matrix 3870 3871 Output Parameter: 3872 . v - the location of the values 3873 3874 3875 Fortran Note: 3876 This routine is used differently from Fortran, e.g., 3877 .vb 3878 Mat mat 3879 PetscScalar mat_array(1) 3880 PetscOffset i_mat 3881 int ierr 3882 call MatGetArray(mat,mat_array,i_mat,ierr) 3883 3884 C Access first local entry in matrix; note that array is 3885 C treated as one dimensional 3886 value = mat_array(i_mat + 1) 3887 3888 [... other code ...] 3889 call MatRestoreArray(mat,mat_array,i_mat,ierr) 3890 .ve 3891 3892 See the Fortran chapter of the users manual and 3893 petsc/src/mat/examples/tests for details. 3894 3895 Level: advanced 3896 3897 Concepts: matrices^access array 3898 3899 .seealso: MatRestoreArray(), MatGetArrayF90() 3900 @*/ 3901 int MatGetArray(Mat mat,PetscScalar *v[]) 3902 { 3903 int ierr; 3904 3905 PetscFunctionBegin; 3906 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3907 PetscValidType(mat); 3908 MatPreallocated(mat); 3909 PetscValidPointer(v); 3910 if (!mat->ops->getarray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3911 ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr); 3912 PetscFunctionReturn(0); 3913 } 3914 3915 #undef __FUNCT__ 3916 #define __FUNCT__ "MatRestoreArray" 3917 /*@C 3918 MatRestoreArray - Restores the matrix after MatGetArray() has been called. 3919 3920 Not Collective 3921 3922 Input Parameter: 3923 + mat - the matrix 3924 - v - the location of the values 3925 3926 Fortran Note: 3927 This routine is used differently from Fortran, e.g., 3928 .vb 3929 Mat mat 3930 PetscScalar mat_array(1) 3931 PetscOffset i_mat 3932 int ierr 3933 call MatGetArray(mat,mat_array,i_mat,ierr) 3934 3935 C Access first local entry in matrix; note that array is 3936 C treated as one dimensional 3937 value = mat_array(i_mat + 1) 3938 3939 [... other code ...] 3940 call MatRestoreArray(mat,mat_array,i_mat,ierr) 3941 .ve 3942 3943 See the Fortran chapter of the users manual and 3944 petsc/src/mat/examples/tests for details 3945 3946 Level: advanced 3947 3948 .seealso: MatGetArray(), MatRestoreArrayF90() 3949 @*/ 3950 int MatRestoreArray(Mat mat,PetscScalar *v[]) 3951 { 3952 int ierr; 3953 3954 PetscFunctionBegin; 3955 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3956 PetscValidType(mat); 3957 MatPreallocated(mat); 3958 PetscValidPointer(v); 3959 #if defined(PETSC_USE_BOPT_g) 3960 CHKMEMQ; 3961 #endif 3962 if (!mat->ops->restorearray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3963 ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr); 3964 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3965 PetscFunctionReturn(0); 3966 } 3967 3968 #undef __FUNCT__ 3969 #define __FUNCT__ "MatGetSubMatrices" 3970 /*@C 3971 MatGetSubMatrices - Extracts several submatrices from a matrix. If submat 3972 points to an array of valid matrices, they may be reused to store the new 3973 submatrices. 3974 3975 Collective on Mat 3976 3977 Input Parameters: 3978 + mat - the matrix 3979 . n - the number of submatrixes to be extracted (on this processor, may be zero) 3980 . irow, icol - index sets of rows and columns to extract 3981 - scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 3982 3983 Output Parameter: 3984 . submat - the array of submatrices 3985 3986 Notes: 3987 MatGetSubMatrices() can extract only sequential submatrices 3988 (from both sequential and parallel matrices). Use MatGetSubMatrix() 3989 to extract a parallel submatrix. 3990 3991 When extracting submatrices from a parallel matrix, each processor can 3992 form a different submatrix by setting the rows and columns of its 3993 individual index sets according to the local submatrix desired. 3994 3995 When finished using the submatrices, the user should destroy 3996 them with MatDestroyMatrices(). 3997 3998 MAT_REUSE_MATRIX can only be used when the nonzero structure of the 3999 original matrix has not changed from that last call to MatGetSubMatrices(). 4000 4001 This routine creates the matrices in submat; you should NOT create them before 4002 calling it. It also allocates the array of matrix pointers submat. 4003 4004 Fortran Note: 4005 The Fortran interface is slightly different from that given below; it 4006 requires one to pass in as submat a Mat (integer) array of size at least m. 4007 4008 Level: advanced 4009 4010 Concepts: matrices^accessing submatrices 4011 Concepts: submatrices 4012 4013 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal() 4014 @*/ 4015 int MatGetSubMatrices(Mat mat,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[]) 4016 { 4017 int ierr; 4018 4019 PetscFunctionBegin; 4020 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4021 PetscValidType(mat); 4022 MatPreallocated(mat); 4023 if (!mat->ops->getsubmatrices) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4024 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4025 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4026 4027 ierr = PetscLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);CHKERRQ(ierr); 4028 ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);CHKERRQ(ierr); 4029 ierr = PetscLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);CHKERRQ(ierr); 4030 PetscFunctionReturn(0); 4031 } 4032 4033 #undef __FUNCT__ 4034 #define __FUNCT__ "MatDestroyMatrices" 4035 /*@C 4036 MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices(). 4037 4038 Collective on Mat 4039 4040 Input Parameters: 4041 + n - the number of local matrices 4042 - mat - the matrices (note that this is a pointer to the array of matrices, just to match the calling 4043 sequence of MatGetSubMatrices()) 4044 4045 Level: advanced 4046 4047 Notes: Frees not only the matrices, but also the array that contains the matrices 4048 4049 .seealso: MatGetSubMatrices() 4050 @*/ 4051 int MatDestroyMatrices(int n,Mat *mat[]) 4052 { 4053 int ierr,i; 4054 4055 PetscFunctionBegin; 4056 if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of matrices %d",n); 4057 PetscValidPointer(mat); 4058 for (i=0; i<n; i++) { 4059 ierr = MatDestroy((*mat)[i]);CHKERRQ(ierr); 4060 } 4061 /* memory is allocated even if n = 0 */ 4062 ierr = PetscFree(*mat);CHKERRQ(ierr); 4063 PetscFunctionReturn(0); 4064 } 4065 4066 #undef __FUNCT__ 4067 #define __FUNCT__ "MatIncreaseOverlap" 4068 /*@ 4069 MatIncreaseOverlap - Given a set of submatrices indicated by index sets, 4070 replaces the index sets by larger ones that represent submatrices with 4071 additional overlap. 4072 4073 Collective on Mat 4074 4075 Input Parameters: 4076 + mat - the matrix 4077 . n - the number of index sets 4078 . is - the array of index sets (these index sets will changed during the call) 4079 - ov - the additional overlap requested 4080 4081 Level: developer 4082 4083 Concepts: overlap 4084 Concepts: ASM^computing overlap 4085 4086 .seealso: MatGetSubMatrices() 4087 @*/ 4088 int MatIncreaseOverlap(Mat mat,int n,IS is[],int ov) 4089 { 4090 int ierr; 4091 4092 PetscFunctionBegin; 4093 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4094 PetscValidType(mat); 4095 MatPreallocated(mat); 4096 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4097 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4098 4099 if (!ov) PetscFunctionReturn(0); 4100 if (!mat->ops->increaseoverlap) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4101 ierr = PetscLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);CHKERRQ(ierr); 4102 ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov);CHKERRQ(ierr); 4103 ierr = PetscLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);CHKERRQ(ierr); 4104 PetscFunctionReturn(0); 4105 } 4106 4107 #undef __FUNCT__ 4108 #define __FUNCT__ "MatPrintHelp" 4109 /*@ 4110 MatPrintHelp - Prints all the options for the matrix. 4111 4112 Collective on Mat 4113 4114 Input Parameter: 4115 . mat - the matrix 4116 4117 Options Database Keys: 4118 + -help - Prints matrix options 4119 - -h - Prints matrix options 4120 4121 Level: developer 4122 4123 .seealso: MatCreate(), MatCreateXXX() 4124 @*/ 4125 int MatPrintHelp(Mat mat) 4126 { 4127 static PetscTruth called = PETSC_FALSE; 4128 int ierr; 4129 4130 PetscFunctionBegin; 4131 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4132 PetscValidType(mat); 4133 MatPreallocated(mat); 4134 4135 if (!called) { 4136 if (mat->ops->printhelp) { 4137 ierr = (*mat->ops->printhelp)(mat);CHKERRQ(ierr); 4138 } 4139 called = PETSC_TRUE; 4140 } 4141 PetscFunctionReturn(0); 4142 } 4143 4144 #undef __FUNCT__ 4145 #define __FUNCT__ "MatGetBlockSize" 4146 /*@ 4147 MatGetBlockSize - Returns the matrix block size; useful especially for the 4148 block row and block diagonal formats. 4149 4150 Not Collective 4151 4152 Input Parameter: 4153 . mat - the matrix 4154 4155 Output Parameter: 4156 . bs - block size 4157 4158 Notes: 4159 Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG. 4160 Block row formats are MATSEQBAIJ, MATMPIBAIJ 4161 4162 Level: intermediate 4163 4164 Concepts: matrices^block size 4165 4166 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag() 4167 @*/ 4168 int MatGetBlockSize(Mat mat,int *bs) 4169 { 4170 int ierr; 4171 4172 PetscFunctionBegin; 4173 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4174 PetscValidType(mat); 4175 MatPreallocated(mat); 4176 PetscValidIntPointer(bs); 4177 if (!mat->ops->getblocksize) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4178 ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr); 4179 PetscFunctionReturn(0); 4180 } 4181 4182 #undef __FUNCT__ 4183 #define __FUNCT__ "MatGetRowIJ" 4184 /*@C 4185 MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices. 4186 4187 Collective on Mat 4188 4189 Input Parameters: 4190 + mat - the matrix 4191 . shift - 0 or 1 indicating we want the indices starting at 0 or 1 4192 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4193 symmetrized 4194 4195 Output Parameters: 4196 + n - number of rows in the (possibly compressed) matrix 4197 . ia - the row pointers 4198 . ja - the column indices 4199 - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned 4200 4201 Level: developer 4202 4203 .seealso: MatGetColumnIJ(), MatRestoreRowIJ() 4204 @*/ 4205 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int *ia[],int* ja[],PetscTruth *done) 4206 { 4207 int ierr; 4208 4209 PetscFunctionBegin; 4210 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4211 PetscValidType(mat); 4212 MatPreallocated(mat); 4213 if (ia) PetscValidIntPointer(ia); 4214 if (ja) PetscValidIntPointer(ja); 4215 PetscValidIntPointer(done); 4216 if (!mat->ops->getrowij) *done = PETSC_FALSE; 4217 else { 4218 *done = PETSC_TRUE; 4219 ierr = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4220 } 4221 PetscFunctionReturn(0); 4222 } 4223 4224 #undef __FUNCT__ 4225 #define __FUNCT__ "MatGetColumnIJ" 4226 /*@C 4227 MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices. 4228 4229 Collective on Mat 4230 4231 Input Parameters: 4232 + mat - the matrix 4233 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4234 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4235 symmetrized 4236 4237 Output Parameters: 4238 + n - number of columns in the (possibly compressed) matrix 4239 . ia - the column pointers 4240 . ja - the row indices 4241 - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned 4242 4243 Level: developer 4244 4245 .seealso: MatGetRowIJ(), MatRestoreColumnIJ() 4246 @*/ 4247 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int *ia[],int* ja[],PetscTruth *done) 4248 { 4249 int ierr; 4250 4251 PetscFunctionBegin; 4252 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4253 PetscValidType(mat); 4254 MatPreallocated(mat); 4255 if (ia) PetscValidIntPointer(ia); 4256 if (ja) PetscValidIntPointer(ja); 4257 PetscValidIntPointer(done); 4258 4259 if (!mat->ops->getcolumnij) *done = PETSC_FALSE; 4260 else { 4261 *done = PETSC_TRUE; 4262 ierr = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4263 } 4264 PetscFunctionReturn(0); 4265 } 4266 4267 #undef __FUNCT__ 4268 #define __FUNCT__ "MatRestoreRowIJ" 4269 /*@C 4270 MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with 4271 MatGetRowIJ(). 4272 4273 Collective on Mat 4274 4275 Input Parameters: 4276 + mat - the matrix 4277 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4278 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4279 symmetrized 4280 4281 Output Parameters: 4282 + n - size of (possibly compressed) matrix 4283 . ia - the row pointers 4284 . ja - the column indices 4285 - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 4286 4287 Level: developer 4288 4289 .seealso: MatGetRowIJ(), MatRestoreColumnIJ() 4290 @*/ 4291 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int *ia[],int* ja[],PetscTruth *done) 4292 { 4293 int ierr; 4294 4295 PetscFunctionBegin; 4296 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4297 PetscValidType(mat); 4298 MatPreallocated(mat); 4299 if (ia) PetscValidIntPointer(ia); 4300 if (ja) PetscValidIntPointer(ja); 4301 PetscValidIntPointer(done); 4302 4303 if (!mat->ops->restorerowij) *done = PETSC_FALSE; 4304 else { 4305 *done = PETSC_TRUE; 4306 ierr = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4307 } 4308 PetscFunctionReturn(0); 4309 } 4310 4311 #undef __FUNCT__ 4312 #define __FUNCT__ "MatRestoreColumnIJ" 4313 /*@C 4314 MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with 4315 MatGetColumnIJ(). 4316 4317 Collective on Mat 4318 4319 Input Parameters: 4320 + mat - the matrix 4321 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4322 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4323 symmetrized 4324 4325 Output Parameters: 4326 + n - size of (possibly compressed) matrix 4327 . ia - the column pointers 4328 . ja - the row indices 4329 - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 4330 4331 Level: developer 4332 4333 .seealso: MatGetColumnIJ(), MatRestoreRowIJ() 4334 @*/ 4335 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int *ia[],int* ja[],PetscTruth *done) 4336 { 4337 int ierr; 4338 4339 PetscFunctionBegin; 4340 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4341 PetscValidType(mat); 4342 MatPreallocated(mat); 4343 if (ia) PetscValidIntPointer(ia); 4344 if (ja) PetscValidIntPointer(ja); 4345 PetscValidIntPointer(done); 4346 4347 if (!mat->ops->restorecolumnij) *done = PETSC_FALSE; 4348 else { 4349 *done = PETSC_TRUE; 4350 ierr = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4351 } 4352 PetscFunctionReturn(0); 4353 } 4354 4355 #undef __FUNCT__ 4356 #define __FUNCT__ "MatColoringPatch" 4357 /*@C 4358 MatColoringPatch -Used inside matrix coloring routines that 4359 use MatGetRowIJ() and/or MatGetColumnIJ(). 4360 4361 Collective on Mat 4362 4363 Input Parameters: 4364 + mat - the matrix 4365 . n - number of colors 4366 - colorarray - array indicating color for each column 4367 4368 Output Parameters: 4369 . iscoloring - coloring generated using colorarray information 4370 4371 Level: developer 4372 4373 .seealso: MatGetRowIJ(), MatGetColumnIJ() 4374 4375 @*/ 4376 int MatColoringPatch(Mat mat,int n,int ncolors,const ISColoringValue colorarray[],ISColoring *iscoloring) 4377 { 4378 int ierr; 4379 4380 PetscFunctionBegin; 4381 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4382 PetscValidType(mat); 4383 MatPreallocated(mat); 4384 PetscValidIntPointer(colorarray); 4385 4386 if (!mat->ops->coloringpatch){ 4387 ierr = ISColoringCreate(mat->comm,n,colorarray,iscoloring);CHKERRQ(ierr); 4388 } else { 4389 ierr = (*mat->ops->coloringpatch)(mat,n,ncolors,colorarray,iscoloring);CHKERRQ(ierr); 4390 } 4391 PetscFunctionReturn(0); 4392 } 4393 4394 4395 #undef __FUNCT__ 4396 #define __FUNCT__ "MatSetUnfactored" 4397 /*@ 4398 MatSetUnfactored - Resets a factored matrix to be treated as unfactored. 4399 4400 Collective on Mat 4401 4402 Input Parameter: 4403 . mat - the factored matrix to be reset 4404 4405 Notes: 4406 This routine should be used only with factored matrices formed by in-place 4407 factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE 4408 format). This option can save memory, for example, when solving nonlinear 4409 systems with a matrix-free Newton-Krylov method and a matrix-based, in-place 4410 ILU(0) preconditioner. 4411 4412 Note that one can specify in-place ILU(0) factorization by calling 4413 .vb 4414 PCType(pc,PCILU); 4415 PCILUSeUseInPlace(pc); 4416 .ve 4417 or by using the options -pc_type ilu -pc_ilu_in_place 4418 4419 In-place factorization ILU(0) can also be used as a local 4420 solver for the blocks within the block Jacobi or additive Schwarz 4421 methods (runtime option: -sub_pc_ilu_in_place). See the discussion 4422 of these preconditioners in the users manual for details on setting 4423 local solver options. 4424 4425 Most users should employ the simplified SLES interface for linear solvers 4426 instead of working directly with matrix algebra routines such as this. 4427 See, e.g., SLESCreate(). 4428 4429 Level: developer 4430 4431 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace() 4432 4433 Concepts: matrices^unfactored 4434 4435 @*/ 4436 int MatSetUnfactored(Mat mat) 4437 { 4438 int ierr; 4439 4440 PetscFunctionBegin; 4441 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4442 PetscValidType(mat); 4443 MatPreallocated(mat); 4444 mat->factor = 0; 4445 if (!mat->ops->setunfactored) PetscFunctionReturn(0); 4446 ierr = (*mat->ops->setunfactored)(mat);CHKERRQ(ierr); 4447 PetscFunctionReturn(0); 4448 } 4449 4450 /*MC 4451 MatGetArrayF90 - Accesses a matrix array from Fortran90. 4452 4453 Synopsis: 4454 MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 4455 4456 Not collective 4457 4458 Input Parameter: 4459 . x - matrix 4460 4461 Output Parameters: 4462 + xx_v - the Fortran90 pointer to the array 4463 - ierr - error code 4464 4465 Example of Usage: 4466 .vb 4467 PetscScalar, pointer xx_v(:) 4468 .... 4469 call MatGetArrayF90(x,xx_v,ierr) 4470 a = xx_v(3) 4471 call MatRestoreArrayF90(x,xx_v,ierr) 4472 .ve 4473 4474 Notes: 4475 Not yet supported for all F90 compilers 4476 4477 Level: advanced 4478 4479 .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray() 4480 4481 Concepts: matrices^accessing array 4482 4483 M*/ 4484 4485 /*MC 4486 MatRestoreArrayF90 - Restores a matrix array that has been 4487 accessed with MatGetArrayF90(). 4488 4489 Synopsis: 4490 MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 4491 4492 Not collective 4493 4494 Input Parameters: 4495 + x - matrix 4496 - xx_v - the Fortran90 pointer to the array 4497 4498 Output Parameter: 4499 . ierr - error code 4500 4501 Example of Usage: 4502 .vb 4503 PetscScalar, pointer xx_v(:) 4504 .... 4505 call MatGetArrayF90(x,xx_v,ierr) 4506 a = xx_v(3) 4507 call MatRestoreArrayF90(x,xx_v,ierr) 4508 .ve 4509 4510 Notes: 4511 Not yet supported for all F90 compilers 4512 4513 Level: advanced 4514 4515 .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray() 4516 4517 M*/ 4518 4519 4520 #undef __FUNCT__ 4521 #define __FUNCT__ "MatGetSubMatrix" 4522 /*@ 4523 MatGetSubMatrix - Gets a single submatrix on the same number of processors 4524 as the original matrix. 4525 4526 Collective on Mat 4527 4528 Input Parameters: 4529 + mat - the original matrix 4530 . isrow - rows this processor should obtain 4531 . iscol - columns for all processors you wish to keep 4532 . csize - number of columns "local" to this processor (does nothing for sequential 4533 matrices). This should match the result from VecGetLocalSize(x,...) if you 4534 plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE 4535 - cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4536 4537 Output Parameter: 4538 . newmat - the new submatrix, of the same type as the old 4539 4540 Level: advanced 4541 4542 Notes: the iscol argument MUST be the same on each processor. You might be 4543 able to create the iscol argument with ISAllGather(). 4544 4545 The first time this is called you should use a cll of MAT_INITIAL_MATRIX, 4546 the MatGetSubMatrix() routine will create the newmat for you. Any additional calls 4547 to this routine with a mat of the same nonzero structure will reuse the matrix 4548 generated the first time. 4549 4550 Concepts: matrices^submatrices 4551 4552 .seealso: MatGetSubMatrices(), ISAllGather() 4553 @*/ 4554 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat) 4555 { 4556 int ierr, size; 4557 Mat *local; 4558 4559 PetscFunctionBegin; 4560 PetscValidType(mat); 4561 MatPreallocated(mat); 4562 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4563 ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr); 4564 4565 /* if original matrix is on just one processor then use submatrix generated */ 4566 if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) { 4567 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr); 4568 PetscFunctionReturn(0); 4569 } else if (!mat->ops->getsubmatrix && size == 1) { 4570 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr); 4571 *newmat = *local; 4572 ierr = PetscFree(local);CHKERRQ(ierr); 4573 PetscFunctionReturn(0); 4574 } 4575 4576 if (!mat->ops->getsubmatrix) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4577 ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr); 4578 ierr = PetscObjectIncreaseState((PetscObject)*newmat); CHKERRQ(ierr); 4579 PetscFunctionReturn(0); 4580 } 4581 4582 #undef __FUNCT__ 4583 #define __FUNCT__ "MatGetPetscMaps" 4584 /*@C 4585 MatGetPetscMaps - Returns the maps associated with the matrix. 4586 4587 Not Collective 4588 4589 Input Parameter: 4590 . mat - the matrix 4591 4592 Output Parameters: 4593 + rmap - the row (right) map 4594 - cmap - the column (left) map 4595 4596 Level: developer 4597 4598 Concepts: maps^getting from matrix 4599 4600 @*/ 4601 int MatGetPetscMaps(Mat mat,PetscMap *rmap,PetscMap *cmap) 4602 { 4603 int ierr; 4604 4605 PetscFunctionBegin; 4606 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4607 PetscValidType(mat); 4608 MatPreallocated(mat); 4609 ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr); 4610 PetscFunctionReturn(0); 4611 } 4612 4613 /* 4614 Version that works for all PETSc matrices 4615 */ 4616 #undef __FUNCT__ 4617 #define __FUNCT__ "MatGetPetscMaps_Petsc" 4618 int MatGetPetscMaps_Petsc(Mat mat,PetscMap *rmap,PetscMap *cmap) 4619 { 4620 PetscFunctionBegin; 4621 if (rmap) *rmap = mat->rmap; 4622 if (cmap) *cmap = mat->cmap; 4623 PetscFunctionReturn(0); 4624 } 4625 4626 #undef __FUNCT__ 4627 #define __FUNCT__ "MatSetStashInitialSize" 4628 /*@ 4629 MatSetStashInitialSize - sets the sizes of the matrix stash, that is 4630 used during the assembly process to store values that belong to 4631 other processors. 4632 4633 Not Collective 4634 4635 Input Parameters: 4636 + mat - the matrix 4637 . size - the initial size of the stash. 4638 - bsize - the initial size of the block-stash(if used). 4639 4640 Options Database Keys: 4641 + -matstash_initial_size <size> or <size0,size1,...sizep-1> 4642 - -matstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1> 4643 4644 Level: intermediate 4645 4646 Notes: 4647 The block-stash is used for values set with VecSetValuesBlocked() while 4648 the stash is used for values set with VecSetValues() 4649 4650 Run with the option -log_info and look for output of the form 4651 MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs. 4652 to determine the appropriate value, MM, to use for size and 4653 MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs. 4654 to determine the value, BMM to use for bsize 4655 4656 Concepts: stash^setting matrix size 4657 Concepts: matrices^stash 4658 4659 @*/ 4660 int MatSetStashInitialSize(Mat mat,int size, int bsize) 4661 { 4662 int ierr; 4663 4664 PetscFunctionBegin; 4665 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4666 PetscValidType(mat); 4667 MatPreallocated(mat); 4668 ierr = MatStashSetInitialSize_Private(&mat->stash,size);CHKERRQ(ierr); 4669 ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize);CHKERRQ(ierr); 4670 PetscFunctionReturn(0); 4671 } 4672 4673 #undef __FUNCT__ 4674 #define __FUNCT__ "MatInterpolateAdd" 4675 /*@ 4676 MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of 4677 the matrix 4678 4679 Collective on Mat 4680 4681 Input Parameters: 4682 + mat - the matrix 4683 . x,y - the vectors 4684 - w - where the result is stored 4685 4686 Level: intermediate 4687 4688 Notes: 4689 w may be the same vector as y. 4690 4691 This allows one to use either the restriction or interpolation (its transpose) 4692 matrix to do the interpolation 4693 4694 Concepts: interpolation 4695 4696 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict() 4697 4698 @*/ 4699 int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w) 4700 { 4701 int M,N,ierr; 4702 4703 PetscFunctionBegin; 4704 PetscValidType(A); 4705 MatPreallocated(A); 4706 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4707 if (N > M) { 4708 ierr = MatMultTransposeAdd(A,x,y,w);CHKERRQ(ierr); 4709 } else { 4710 ierr = MatMultAdd(A,x,y,w);CHKERRQ(ierr); 4711 } 4712 PetscFunctionReturn(0); 4713 } 4714 4715 #undef __FUNCT__ 4716 #define __FUNCT__ "MatInterpolate" 4717 /*@ 4718 MatInterpolate - y = A*x or A'*x depending on the shape of 4719 the matrix 4720 4721 Collective on Mat 4722 4723 Input Parameters: 4724 + mat - the matrix 4725 - x,y - the vectors 4726 4727 Level: intermediate 4728 4729 Notes: 4730 This allows one to use either the restriction or interpolation (its transpose) 4731 matrix to do the interpolation 4732 4733 Concepts: matrices^interpolation 4734 4735 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict() 4736 4737 @*/ 4738 int MatInterpolate(Mat A,Vec x,Vec y) 4739 { 4740 int M,N,ierr; 4741 4742 PetscFunctionBegin; 4743 PetscValidType(A); 4744 MatPreallocated(A); 4745 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4746 if (N > M) { 4747 ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr); 4748 } else { 4749 ierr = MatMult(A,x,y);CHKERRQ(ierr); 4750 } 4751 PetscFunctionReturn(0); 4752 } 4753 4754 #undef __FUNCT__ 4755 #define __FUNCT__ "MatRestrict" 4756 /*@ 4757 MatRestrict - y = A*x or A'*x 4758 4759 Collective on Mat 4760 4761 Input Parameters: 4762 + mat - the matrix 4763 - x,y - the vectors 4764 4765 Level: intermediate 4766 4767 Notes: 4768 This allows one to use either the restriction or interpolation (its transpose) 4769 matrix to do the restriction 4770 4771 Concepts: matrices^restriction 4772 4773 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate() 4774 4775 @*/ 4776 int MatRestrict(Mat A,Vec x,Vec y) 4777 { 4778 int M,N,ierr; 4779 4780 PetscFunctionBegin; 4781 PetscValidType(A); 4782 MatPreallocated(A); 4783 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4784 if (N > M) { 4785 ierr = MatMult(A,x,y);CHKERRQ(ierr); 4786 } else { 4787 ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr); 4788 } 4789 PetscFunctionReturn(0); 4790 } 4791 4792 #undef __FUNCT__ 4793 #define __FUNCT__ "MatNullSpaceAttach" 4794 /*@C 4795 MatNullSpaceAttach - attaches a null space to a matrix. 4796 This null space will be removed from the resulting vector whenever 4797 MatMult() is called 4798 4799 Collective on Mat 4800 4801 Input Parameters: 4802 + mat - the matrix 4803 - nullsp - the null space object 4804 4805 Level: developer 4806 4807 Notes: 4808 Overwrites any previous null space that may have been attached 4809 4810 Concepts: null space^attaching to matrix 4811 4812 .seealso: MatCreate(), MatNullSpaceCreate() 4813 @*/ 4814 int MatNullSpaceAttach(Mat mat,MatNullSpace nullsp) 4815 { 4816 int ierr; 4817 4818 PetscFunctionBegin; 4819 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4820 PetscValidType(mat); 4821 MatPreallocated(mat); 4822 PetscValidHeaderSpecific(nullsp,MAT_NULLSPACE_COOKIE); 4823 4824 if (mat->nullsp) { 4825 ierr = MatNullSpaceDestroy(mat->nullsp);CHKERRQ(ierr); 4826 } 4827 mat->nullsp = nullsp; 4828 ierr = PetscObjectReference((PetscObject)nullsp);CHKERRQ(ierr); 4829 PetscFunctionReturn(0); 4830 } 4831 4832 #undef __FUNCT__ 4833 #define __FUNCT__ "MatICCFactor" 4834 /*@ 4835 MatICCFactor - Performs in-place incomplete Cholesky factorization of matrix. 4836 4837 Collective on Mat 4838 4839 Input Parameters: 4840 + mat - the matrix 4841 . row - row/column permutation 4842 . fill - expected fill factor >= 1.0 4843 - level - level of fill, for ICC(k) 4844 4845 Notes: 4846 Probably really in-place only when level of fill is zero, otherwise allocates 4847 new space to store factored matrix and deletes previous memory. 4848 4849 Most users should employ the simplified SLES interface for linear solvers 4850 instead of working directly with matrix algebra routines such as this. 4851 See, e.g., SLESCreate(). 4852 4853 Level: developer 4854 4855 Concepts: matrices^incomplete Cholesky factorization 4856 Concepts: Cholesky factorization 4857 4858 .seealso: MatICCFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 4859 @*/ 4860 int MatICCFactor(Mat mat,IS row,MatFactorInfo* info) 4861 { 4862 int ierr; 4863 4864 PetscFunctionBegin; 4865 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4866 PetscValidType(mat); 4867 MatPreallocated(mat); 4868 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 4869 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4870 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4871 if (!mat->ops->iccfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4872 ierr = (*mat->ops->iccfactor)(mat,row,info);CHKERRQ(ierr); 4873 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 4874 PetscFunctionReturn(0); 4875 } 4876 4877 #undef __FUNCT__ 4878 #define __FUNCT__ "MatSetValuesAdic" 4879 /*@ 4880 MatSetValuesAdic - Sets values computed with ADIC automatic differentiation into a matrix. 4881 4882 Not Collective 4883 4884 Input Parameters: 4885 + mat - the matrix 4886 - v - the values compute with ADIC 4887 4888 Level: developer 4889 4890 Notes: 4891 Must call MatSetColoring() before using this routine. Also this matrix must already 4892 have its nonzero pattern determined. 4893 4894 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4895 MatSetValues(), MatSetColoring(), MatSetValuesAdifor() 4896 @*/ 4897 int MatSetValuesAdic(Mat mat,void *v) 4898 { 4899 int ierr; 4900 4901 PetscFunctionBegin; 4902 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4903 PetscValidType(mat); 4904 4905 if (!mat->assembled) { 4906 SETERRQ(1,"Matrix must be already assembled"); 4907 } 4908 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4909 if (!mat->ops->setvaluesadic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4910 ierr = (*mat->ops->setvaluesadic)(mat,v);CHKERRQ(ierr); 4911 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4912 ierr = MatView_Private(mat);CHKERRQ(ierr); 4913 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 4914 PetscFunctionReturn(0); 4915 } 4916 4917 4918 #undef __FUNCT__ 4919 #define __FUNCT__ "MatSetColoring" 4920 /*@ 4921 MatSetColoring - Sets a coloring used by calls to MatSetValuesAdic() 4922 4923 Not Collective 4924 4925 Input Parameters: 4926 + mat - the matrix 4927 - coloring - the coloring 4928 4929 Level: developer 4930 4931 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4932 MatSetValues(), MatSetValuesAdic() 4933 @*/ 4934 int MatSetColoring(Mat mat,ISColoring coloring) 4935 { 4936 int ierr; 4937 4938 PetscFunctionBegin; 4939 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4940 PetscValidType(mat); 4941 4942 if (!mat->assembled) { 4943 SETERRQ(1,"Matrix must be already assembled"); 4944 } 4945 if (!mat->ops->setcoloring) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4946 ierr = (*mat->ops->setcoloring)(mat,coloring);CHKERRQ(ierr); 4947 PetscFunctionReturn(0); 4948 } 4949 4950 #undef __FUNCT__ 4951 #define __FUNCT__ "MatSetValuesAdifor" 4952 /*@ 4953 MatSetValuesAdifor - Sets values computed with automatic differentiation into a matrix. 4954 4955 Not Collective 4956 4957 Input Parameters: 4958 + mat - the matrix 4959 . nl - leading dimension of v 4960 - v - the values compute with ADIFOR 4961 4962 Level: developer 4963 4964 Notes: 4965 Must call MatSetColoring() before using this routine. Also this matrix must already 4966 have its nonzero pattern determined. 4967 4968 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4969 MatSetValues(), MatSetColoring() 4970 @*/ 4971 int MatSetValuesAdifor(Mat mat,int nl,void *v) 4972 { 4973 int ierr; 4974 4975 PetscFunctionBegin; 4976 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4977 PetscValidType(mat); 4978 4979 if (!mat->assembled) { 4980 SETERRQ(1,"Matrix must be already assembled"); 4981 } 4982 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4983 if (!mat->ops->setvaluesadifor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4984 ierr = (*mat->ops->setvaluesadifor)(mat,nl,v);CHKERRQ(ierr); 4985 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4986 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 4987 PetscFunctionReturn(0); 4988 } 4989 4990 EXTERN int MatMPIAIJDiagonalScaleLocal(Mat,Vec); 4991 EXTERN int MatMPIBAIJDiagonalScaleLocal(Mat,Vec); 4992 4993 #undef __FUNCT__ 4994 #define __FUNCT__ "MatDiagonalScaleLocal" 4995 /*@ 4996 MatDiagonalScaleLocal - Scales columns of a matrix given the scaling values including the 4997 ghosted ones. 4998 4999 Not Collective 5000 5001 Input Parameters: 5002 + mat - the matrix 5003 - diag = the diagonal values, including ghost ones 5004 5005 Level: developer 5006 5007 Notes: Works only for MPIAIJ and MPIBAIJ matrices 5008 5009 .seealso: MatDiagonalScale() 5010 @*/ 5011 int MatDiagonalScaleLocal(Mat mat,Vec diag) 5012 { 5013 int ierr,size; 5014 5015 PetscFunctionBegin; 5016 PetscValidHeaderSpecific(mat,MAT_COOKIE); 5017 PetscValidHeaderSpecific(diag,VEC_COOKIE); 5018 PetscValidType(mat); 5019 5020 if (!mat->assembled) { 5021 SETERRQ(1,"Matrix must be already assembled"); 5022 } 5023 ierr = PetscLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 5024 ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr); 5025 if (size == 1) { 5026 int n,m; 5027 ierr = VecGetSize(diag,&n);CHKERRQ(ierr); 5028 ierr = MatGetSize(mat,0,&m);CHKERRQ(ierr); 5029 if (m == n) { 5030 ierr = MatDiagonalScale(mat,0,diag);CHKERRQ(ierr); 5031 } else { 5032 SETERRQ(1,"Only supprted for sequential matrices when no ghost points/periodic conditions"); 5033 } 5034 } else { 5035 int (*f)(Mat,Vec); 5036 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",(void (**)(void))&f);CHKERRQ(ierr); 5037 if (f) { 5038 ierr = (*f)(mat,diag);CHKERRQ(ierr); 5039 } else { 5040 SETERRQ(1,"Only supported for MPIAIJ and MPIBAIJ parallel matrices"); 5041 } 5042 } 5043 ierr = PetscLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 5044 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 5045 PetscFunctionReturn(0); 5046 } 5047 5048 #undef __FUNCT__ 5049 #define __FUNCT__ "MatGetInertia" 5050 /*@ 5051 MatGetInertia - Gets the inertia from a factored matrix 5052 5053 Collective on Mat 5054 5055 Input Parameter: 5056 . mat - the matrix 5057 5058 Output Parameters: 5059 + nneg - number of negative eigenvalues 5060 . nzero - number of zero eigenvalues 5061 - npos - number of positive eigenvalues 5062 5063 Level: advanced 5064 5065 Notes: Matrix must have been factored by MatCholeskyFactor() 5066 5067 5068 @*/ 5069 int MatGetInertia(Mat mat,int *nneg,int *nzero,int *npos) 5070 { 5071 int ierr; 5072 5073 PetscFunctionBegin; 5074 PetscValidHeaderSpecific(mat,MAT_COOKIE); 5075 PetscValidType(mat); 5076 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 5077 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Numeric factor mat is not assembled"); 5078 if (!mat->ops->getinertia) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 5079 ierr = (*mat->ops->getinertia)(mat,nneg,nzero,npos);CHKERRQ(ierr); 5080 PetscFunctionReturn(0); 5081 } 5082 5083 /* ----------------------------------------------------------------*/ 5084 #undef __FUNCT__ 5085 #define __FUNCT__ "MatSolves" 5086 /*@ 5087 MatSolves - Solves A x = b, given a factored matrix, for a collection of vectors 5088 5089 Collective on Mat and Vecs 5090 5091 Input Parameters: 5092 + mat - the factored matrix 5093 - b - the right-hand-side vectors 5094 5095 Output Parameter: 5096 . x - the result vectors 5097 5098 Notes: 5099 The vectors b and x cannot be the same. I.e., one cannot 5100 call MatSolves(A,x,x). 5101 5102 Notes: 5103 Most users should employ the simplified SLES interface for linear solvers 5104 instead of working directly with matrix algebra routines such as this. 5105 See, e.g., SLESCreate(). 5106 5107 Level: developer 5108 5109 Concepts: matrices^triangular solves 5110 5111 .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), MatSolve() 5112 @*/ 5113 int MatSolves(Mat mat,Vecs b,Vecs x) 5114 { 5115 int ierr; 5116 5117 PetscFunctionBegin; 5118 PetscValidHeaderSpecific(mat,MAT_COOKIE); 5119 PetscValidType(mat); 5120 MatPreallocated(mat); 5121 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 5122 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 5123 if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0); 5124 5125 if (!mat->ops->solves) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 5126 ierr = PetscLogEventBegin(MAT_Solves,mat,0,0,0);CHKERRQ(ierr); 5127 ierr = (*mat->ops->solves)(mat,b,x);CHKERRQ(ierr); 5128 ierr = PetscLogEventEnd(MAT_Solves,mat,0,0,0);CHKERRQ(ierr); 5129 PetscFunctionReturn(0); 5130 } 5131 5132 #undef __FUNCT__ 5133 #define __FUNCT__ "MatIsSymmetric" 5134 /*@C 5135 MatIsSymmetric - Test whether a matrix is symmetric 5136 5137 Collective on Mat 5138 5139 Input Parameter: 5140 . A - the matrix to test 5141 5142 Output Parameters: 5143 . flg - the result 5144 5145 Level: intermediate 5146 5147 Concepts: matrix^symmetry 5148 5149 .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsStructurallySymmetric(), MatSetOption() 5150 @*/ 5151 int MatIsSymmetric(Mat A,PetscTruth *flg) 5152 { 5153 int ierr; 5154 5155 PetscFunctionBegin; 5156 PetscValidHeaderSpecific(A,MAT_COOKIE); 5157 if (!A->symmetric_set) { 5158 if (!A->ops->issymmetric) SETERRQ(1,"Matrix does not support checking for symmetric"); 5159 ierr = (*A->ops->issymmetric)(A,&A->symmetric);CHKERRQ(ierr); 5160 A->symmetric_set = PETSC_TRUE; 5161 if (A->symmetric) { 5162 A->structurally_symmetric_set = PETSC_TRUE; 5163 A->structurally_symmetric = PETSC_TRUE; 5164 } 5165 } 5166 *flg = A->symmetric; 5167 PetscFunctionReturn(0); 5168 } 5169 5170 #undef __FUNCT__ 5171 #define __FUNCT__ "MatIsStructurallySymmetric" 5172 /*@C 5173 MatIsStructurallySymmetric - Test whether a matrix is structurally symmetric 5174 5175 Collective on Mat 5176 5177 Input Parameter: 5178 . A - the matrix to test 5179 5180 Output Parameters: 5181 . flg - the result 5182 5183 Level: intermediate 5184 5185 Concepts: matrix^symmetry 5186 5187 .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsSymmetric(), MatSetOption() 5188 @*/ 5189 int MatIsStructurallySymmetric(Mat A,PetscTruth *flg) 5190 { 5191 int ierr; 5192 5193 PetscFunctionBegin; 5194 PetscValidHeaderSpecific(A,MAT_COOKIE); 5195 if (!A->structurally_symmetric_set) { 5196 if (!A->ops->isstructurallysymmetric) SETERRQ(1,"Matrix does not support checking for structural symmetric"); 5197 ierr = (*A->ops->isstructurallysymmetric)(A,&A->structurally_symmetric);CHKERRQ(ierr); 5198 A->structurally_symmetric_set = PETSC_TRUE; 5199 } 5200 *flg = A->structurally_symmetric; 5201 PetscFunctionReturn(0); 5202 } 5203 5204 #undef __FUNCT__ 5205 #define __FUNCT__ "MatIsHermitian" 5206 /*@C 5207 MatIsHermitian - Test whether a matrix is Hermitian, i.e. it is the complex conjugate of its transpose. 5208 5209 Collective on Mat 5210 5211 Input Parameter: 5212 . A - the matrix to test 5213 5214 Output Parameters: 5215 . flg - the result 5216 5217 Level: intermediate 5218 5219 Concepts: matrix^symmetry 5220 5221 .seealso: MatTranspose(), MatIsTranspose(), MatIsSymmetric(), MatIsStructurallySymmetric(), MatSetOption() 5222 @*/ 5223 int MatIsHermitian(Mat A,PetscTruth *flg) 5224 { 5225 int ierr; 5226 5227 PetscFunctionBegin; 5228 PetscValidHeaderSpecific(A,MAT_COOKIE); 5229 if (!A->hermitian_set) { 5230 if (!A->ops->ishermitian) SETERRQ(1,"Matrix does not support checking for being Hermitian"); 5231 ierr = (*A->ops->ishermitian)(A,&A->hermitian);CHKERRQ(ierr); 5232 A->hermitian_set = PETSC_TRUE; 5233 if (A->hermitian) { 5234 A->structurally_symmetric_set = PETSC_TRUE; 5235 A->structurally_symmetric = PETSC_TRUE; 5236 } 5237 } 5238 *flg = A->hermitian; 5239 PetscFunctionReturn(0); 5240 } 5241