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 KSP 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 KSP 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 ierr = PetscObjectIncreaseState((PetscObject)y); CHKERRQ(ierr); 1126 PetscFunctionReturn(0); 1127 } 1128 1129 #undef __FUNCT__ 1130 #define __FUNCT__ "MatMultTranspose" 1131 /*@ 1132 MatMultTranspose - Computes matrix transpose times a vector. 1133 1134 Collective on Mat and Vec 1135 1136 Input Parameters: 1137 + mat - the matrix 1138 - x - the vector to be multilplied 1139 1140 Output Parameters: 1141 . y - the result 1142 1143 Notes: 1144 The vectors x and y cannot be the same. I.e., one cannot 1145 call MatMultTranspose(A,y,y). 1146 1147 Level: beginner 1148 1149 Concepts: matrix vector product^transpose 1150 1151 .seealso: MatMult(), MatMultAdd(), MatMultTransposeAdd() 1152 @*/ 1153 int MatMultTranspose(Mat mat,Vec x,Vec y) 1154 { 1155 int ierr; 1156 PetscTruth flg1, flg2; 1157 1158 PetscFunctionBegin; 1159 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1160 PetscValidType(mat); 1161 MatPreallocated(mat); 1162 PetscValidHeaderSpecific(x,VEC_COOKIE); 1163 PetscValidHeaderSpecific(y,VEC_COOKIE); 1164 1165 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1166 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1167 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1168 #ifndef PETSC_HAVE_CONSTRAINTS 1169 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 1170 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N); 1171 #endif 1172 1173 if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP, "Operation not supported"); 1174 ierr = PetscLogEventBegin(MAT_MultTranspose,mat,x,y,0);CHKERRQ(ierr); 1175 if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined"); 1176 1177 ierr = PetscTypeCompare((PetscObject)mat,MATSEQSBAIJ,&flg1); 1178 ierr = PetscTypeCompare((PetscObject)mat,MATMPISBAIJ,&flg2); 1179 if (flg1 || flg2) { /* mat is in sbaij format */ 1180 ierr = (*mat->ops->mult)(mat,x,y);CHKERRQ(ierr); 1181 } else { 1182 ierr = (*mat->ops->multtranspose)(mat,x,y);CHKERRQ(ierr); 1183 } 1184 ierr = PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);CHKERRQ(ierr); 1185 ierr = PetscObjectIncreaseState((PetscObject)y); CHKERRQ(ierr); 1186 PetscFunctionReturn(0); 1187 } 1188 1189 #undef __FUNCT__ 1190 #define __FUNCT__ "MatMultAdd" 1191 /*@ 1192 MatMultAdd - Computes v3 = v2 + A * v1. 1193 1194 Collective on Mat and Vec 1195 1196 Input Parameters: 1197 + mat - the matrix 1198 - v1, v2 - the vectors 1199 1200 Output Parameters: 1201 . v3 - the result 1202 1203 Notes: 1204 The vectors v1 and v3 cannot be the same. I.e., one cannot 1205 call MatMultAdd(A,v1,v2,v1). 1206 1207 Level: beginner 1208 1209 Concepts: matrix vector product^addition 1210 1211 .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd() 1212 @*/ 1213 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3) 1214 { 1215 int ierr; 1216 1217 PetscFunctionBegin; 1218 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1219 PetscValidType(mat); 1220 MatPreallocated(mat); 1221 PetscValidHeaderSpecific(v1,VEC_COOKIE); 1222 PetscValidHeaderSpecific(v2,VEC_COOKIE); 1223 PetscValidHeaderSpecific(v3,VEC_COOKIE); 1224 1225 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1226 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1227 if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N); 1228 if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N); 1229 if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N); 1230 if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n); 1231 if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n); 1232 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors"); 1233 1234 ierr = PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1235 ierr = (*mat->ops->multadd)(mat,v1,v2,v3);CHKERRQ(ierr); 1236 ierr = PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1237 ierr = PetscObjectIncreaseState((PetscObject)v3); CHKERRQ(ierr); 1238 PetscFunctionReturn(0); 1239 } 1240 1241 #undef __FUNCT__ 1242 #define __FUNCT__ "MatMultTransposeAdd" 1243 /*@ 1244 MatMultTransposeAdd - Computes v3 = v2 + A' * v1. 1245 1246 Collective on Mat and Vec 1247 1248 Input Parameters: 1249 + mat - the matrix 1250 - v1, v2 - the vectors 1251 1252 Output Parameters: 1253 . v3 - the result 1254 1255 Notes: 1256 The vectors v1 and v3 cannot be the same. I.e., one cannot 1257 call MatMultTransposeAdd(A,v1,v2,v1). 1258 1259 Level: beginner 1260 1261 Concepts: matrix vector product^transpose and addition 1262 1263 .seealso: MatMultTranspose(), MatMultAdd(), MatMult() 1264 @*/ 1265 int MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3) 1266 { 1267 int ierr; 1268 1269 PetscFunctionBegin; 1270 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1271 PetscValidType(mat); 1272 MatPreallocated(mat); 1273 PetscValidHeaderSpecific(v1,VEC_COOKIE); 1274 PetscValidHeaderSpecific(v2,VEC_COOKIE); 1275 PetscValidHeaderSpecific(v3,VEC_COOKIE); 1276 1277 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1278 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1279 if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1280 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors"); 1281 if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N); 1282 if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N); 1283 if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N); 1284 1285 ierr = PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1286 ierr = (*mat->ops->multtransposeadd)(mat,v1,v2,v3);CHKERRQ(ierr); 1287 ierr = PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1288 ierr = PetscObjectIncreaseState((PetscObject)v3); CHKERRQ(ierr); 1289 PetscFunctionReturn(0); 1290 } 1291 1292 #undef __FUNCT__ 1293 #define __FUNCT__ "MatMultConstrained" 1294 /*@ 1295 MatMultConstrained - The inner multiplication routine for a 1296 constrained matrix P^T A P. 1297 1298 Collective on Mat and Vec 1299 1300 Input Parameters: 1301 + mat - the matrix 1302 - x - the vector to be multilplied 1303 1304 Output Parameters: 1305 . y - the result 1306 1307 Notes: 1308 The vectors x and y cannot be the same. I.e., one cannot 1309 call MatMult(A,y,y). 1310 1311 Level: beginner 1312 1313 .keywords: matrix, multiply, matrix-vector product, constraint 1314 .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd() 1315 @*/ 1316 int MatMultConstrained(Mat mat,Vec x,Vec y) 1317 { 1318 int ierr; 1319 1320 PetscFunctionBegin; 1321 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1322 PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1323 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1324 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1325 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1326 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1327 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 1328 if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n); 1329 1330 ierr = PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1331 ierr = (*mat->ops->multconstrained)(mat,x,y); CHKERRQ(ierr); 1332 ierr = PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1333 ierr = PetscObjectIncreaseState((PetscObject)y); CHKERRQ(ierr); 1334 1335 PetscFunctionReturn(0); 1336 } 1337 1338 #undef __FUNCT__ 1339 #define __FUNCT__ "MatMultTransposeConstrained" 1340 /*@ 1341 MatMultTransposeConstrained - The inner multiplication routine for a 1342 constrained matrix P^T A^T P. 1343 1344 Collective on Mat and Vec 1345 1346 Input Parameters: 1347 + mat - the matrix 1348 - x - the vector to be multilplied 1349 1350 Output Parameters: 1351 . y - the result 1352 1353 Notes: 1354 The vectors x and y cannot be the same. I.e., one cannot 1355 call MatMult(A,y,y). 1356 1357 Level: beginner 1358 1359 .keywords: matrix, multiply, matrix-vector product, constraint 1360 .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd() 1361 @*/ 1362 int MatMultTransposeConstrained(Mat mat,Vec x,Vec y) 1363 { 1364 int ierr; 1365 1366 PetscFunctionBegin; 1367 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1368 PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1369 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1370 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1371 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1372 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1373 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 1374 1375 ierr = PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1376 ierr = (*mat->ops->multtransposeconstrained)(mat,x,y);CHKERRQ(ierr); 1377 ierr = PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1378 ierr = PetscObjectIncreaseState((PetscObject)y); CHKERRQ(ierr); 1379 1380 PetscFunctionReturn(0); 1381 } 1382 /* ------------------------------------------------------------*/ 1383 #undef __FUNCT__ 1384 #define __FUNCT__ "MatGetInfo" 1385 /*@C 1386 MatGetInfo - Returns information about matrix storage (number of 1387 nonzeros, memory, etc.). 1388 1389 Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used 1390 as the flag 1391 1392 Input Parameters: 1393 . mat - the matrix 1394 1395 Output Parameters: 1396 + flag - flag indicating the type of parameters to be returned 1397 (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors, 1398 MAT_GLOBAL_SUM - sum over all processors) 1399 - info - matrix information context 1400 1401 Notes: 1402 The MatInfo context contains a variety of matrix data, including 1403 number of nonzeros allocated and used, number of mallocs during 1404 matrix assembly, etc. Additional information for factored matrices 1405 is provided (such as the fill ratio, number of mallocs during 1406 factorization, etc.). Much of this info is printed to STDOUT 1407 when using the runtime options 1408 $ -log_info -mat_view_info 1409 1410 Example for C/C++ Users: 1411 See the file ${PETSC_DIR}/include/petscmat.h for a complete list of 1412 data within the MatInfo context. For example, 1413 .vb 1414 MatInfo info; 1415 Mat A; 1416 double mal, nz_a, nz_u; 1417 1418 MatGetInfo(A,MAT_LOCAL,&info); 1419 mal = info.mallocs; 1420 nz_a = info.nz_allocated; 1421 .ve 1422 1423 Example for Fortran Users: 1424 Fortran users should declare info as a double precision 1425 array of dimension MAT_INFO_SIZE, and then extract the parameters 1426 of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h 1427 a complete list of parameter names. 1428 .vb 1429 double precision info(MAT_INFO_SIZE) 1430 double precision mal, nz_a 1431 Mat A 1432 integer ierr 1433 1434 call MatGetInfo(A,MAT_LOCAL,info,ierr) 1435 mal = info(MAT_INFO_MALLOCS) 1436 nz_a = info(MAT_INFO_NZ_ALLOCATED) 1437 .ve 1438 1439 Level: intermediate 1440 1441 Concepts: matrices^getting information on 1442 1443 @*/ 1444 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info) 1445 { 1446 int ierr; 1447 1448 PetscFunctionBegin; 1449 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1450 PetscValidType(mat); 1451 MatPreallocated(mat); 1452 PetscValidPointer(info); 1453 if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1454 ierr = (*mat->ops->getinfo)(mat,flag,info);CHKERRQ(ierr); 1455 PetscFunctionReturn(0); 1456 } 1457 1458 /* ----------------------------------------------------------*/ 1459 #undef __FUNCT__ 1460 #define __FUNCT__ "MatILUDTFactor" 1461 /*@C 1462 MatILUDTFactor - Performs a drop tolerance ILU factorization. 1463 1464 Collective on Mat 1465 1466 Input Parameters: 1467 + mat - the matrix 1468 . info - information about the factorization to be done 1469 . row - row permutation 1470 - col - column permutation 1471 1472 Output Parameters: 1473 . fact - the factored matrix 1474 1475 Level: developer 1476 1477 Notes: 1478 Most users should employ the simplified KSP interface for linear solvers 1479 instead of working directly with matrix algebra routines such as this. 1480 See, e.g., KSPCreate(). 1481 1482 This is currently only supported for the SeqAIJ matrix format using code 1483 from Yousef Saad's SPARSEKIT2 package (translated to C with f2c) and/or 1484 Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright 1485 and thus can be distributed with PETSc. 1486 1487 Concepts: matrices^ILUDT factorization 1488 1489 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo 1490 @*/ 1491 int MatILUDTFactor(Mat mat,MatFactorInfo *info,IS row,IS col,Mat *fact) 1492 { 1493 int ierr; 1494 1495 PetscFunctionBegin; 1496 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1497 PetscValidType(mat); 1498 MatPreallocated(mat); 1499 PetscValidPointer(fact); 1500 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1501 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1502 if (!mat->ops->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1503 1504 ierr = PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1505 ierr = (*mat->ops->iludtfactor)(mat,info,row,col,fact);CHKERRQ(ierr); 1506 ierr = PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1507 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1508 1509 PetscFunctionReturn(0); 1510 } 1511 1512 #undef __FUNCT__ 1513 #define __FUNCT__ "MatLUFactor" 1514 /*@ 1515 MatLUFactor - Performs in-place LU factorization of matrix. 1516 1517 Collective on Mat 1518 1519 Input Parameters: 1520 + mat - the matrix 1521 . row - row permutation 1522 . col - column permutation 1523 - info - options for factorization, includes 1524 $ fill - expected fill as ratio of original fill. 1525 $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting) 1526 $ Run with the option -log_info to determine an optimal value to use 1527 1528 Notes: 1529 Most users should employ the simplified KSP interface for linear solvers 1530 instead of working directly with matrix algebra routines such as this. 1531 See, e.g., KSPCreate(). 1532 1533 This changes the state of the matrix to a factored matrix; it cannot be used 1534 for example with MatSetValues() unless one first calls MatSetUnfactored(). 1535 1536 Level: developer 1537 1538 Concepts: matrices^LU factorization 1539 1540 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), 1541 MatGetOrdering(), MatSetUnfactored(), MatFactorInfo 1542 1543 @*/ 1544 int MatLUFactor(Mat mat,IS row,IS col,MatFactorInfo *info) 1545 { 1546 int ierr; 1547 1548 PetscFunctionBegin; 1549 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1550 PetscValidType(mat); 1551 MatPreallocated(mat); 1552 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1553 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1554 if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1555 1556 ierr = PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);CHKERRQ(ierr); 1557 ierr = (*mat->ops->lufactor)(mat,row,col,info);CHKERRQ(ierr); 1558 ierr = PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);CHKERRQ(ierr); 1559 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 1560 PetscFunctionReturn(0); 1561 } 1562 1563 #undef __FUNCT__ 1564 #define __FUNCT__ "MatILUFactor" 1565 /*@ 1566 MatILUFactor - Performs in-place ILU factorization of matrix. 1567 1568 Collective on Mat 1569 1570 Input Parameters: 1571 + mat - the matrix 1572 . row - row permutation 1573 . col - column permutation 1574 - info - structure containing 1575 $ levels - number of levels of fill. 1576 $ expected fill - as ratio of original fill. 1577 $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices 1578 missing diagonal entries) 1579 1580 Notes: 1581 Probably really in-place only when level of fill is zero, otherwise allocates 1582 new space to store factored matrix and deletes previous memory. 1583 1584 Most users should employ the simplified KSP interface for linear solvers 1585 instead of working directly with matrix algebra routines such as this. 1586 See, e.g., KSPCreate(). 1587 1588 Level: developer 1589 1590 Concepts: matrices^ILU factorization 1591 1592 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo 1593 @*/ 1594 int MatILUFactor(Mat mat,IS row,IS col,MatFactorInfo *info) 1595 { 1596 int ierr; 1597 1598 PetscFunctionBegin; 1599 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1600 PetscValidType(mat); 1601 MatPreallocated(mat); 1602 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 1603 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1604 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1605 if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1606 1607 ierr = PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1608 ierr = (*mat->ops->ilufactor)(mat,row,col,info);CHKERRQ(ierr); 1609 ierr = PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1610 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 #undef __FUNCT__ 1615 #define __FUNCT__ "MatLUFactorSymbolic" 1616 /*@ 1617 MatLUFactorSymbolic - Performs symbolic LU factorization of matrix. 1618 Call this routine before calling MatLUFactorNumeric(). 1619 1620 Collective on Mat 1621 1622 Input Parameters: 1623 + mat - the matrix 1624 . row, col - row and column permutations 1625 - info - options for factorization, includes 1626 $ fill - expected fill as ratio of original fill. 1627 $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting) 1628 $ Run with the option -log_info to determine an optimal value to use 1629 1630 Output Parameter: 1631 . fact - new matrix that has been symbolically factored 1632 1633 Notes: 1634 See the users manual for additional information about 1635 choosing the fill factor for better efficiency. 1636 1637 Most users should employ the simplified KSP interface for linear solvers 1638 instead of working directly with matrix algebra routines such as this. 1639 See, e.g., KSPCreate(). 1640 1641 Level: developer 1642 1643 Concepts: matrices^LU symbolic factorization 1644 1645 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatFactorInfo 1646 @*/ 1647 int MatLUFactorSymbolic(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact) 1648 { 1649 int ierr; 1650 1651 PetscFunctionBegin; 1652 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1653 PetscValidType(mat); 1654 MatPreallocated(mat); 1655 PetscValidPointer(fact); 1656 PetscValidHeaderSpecific(row,IS_COOKIE); 1657 PetscValidHeaderSpecific(col,IS_COOKIE); 1658 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1659 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1660 if (!mat->ops->lufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic LU",mat->type_name); 1661 1662 ierr = PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 1663 ierr = (*mat->ops->lufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr); 1664 ierr = PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 1665 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1666 PetscFunctionReturn(0); 1667 } 1668 1669 #undef __FUNCT__ 1670 #define __FUNCT__ "MatLUFactorNumeric" 1671 /*@ 1672 MatLUFactorNumeric - Performs numeric LU factorization of a matrix. 1673 Call this routine after first calling MatLUFactorSymbolic(). 1674 1675 Collective on Mat 1676 1677 Input Parameters: 1678 + mat - the matrix 1679 - fact - the matrix generated for the factor, from MatLUFactorSymbolic() 1680 1681 Notes: 1682 See MatLUFactor() for in-place factorization. See 1683 MatCholeskyFactorNumeric() for the symmetric, positive definite case. 1684 1685 Most users should employ the simplified KSP interface for linear solvers 1686 instead of working directly with matrix algebra routines such as this. 1687 See, e.g., KSPCreate(). 1688 1689 Level: developer 1690 1691 Concepts: matrices^LU numeric factorization 1692 1693 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor() 1694 @*/ 1695 int MatLUFactorNumeric(Mat mat,Mat *fact) 1696 { 1697 int ierr; 1698 1699 PetscFunctionBegin; 1700 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1701 PetscValidType(mat); 1702 MatPreallocated(mat); 1703 PetscValidPointer(fact); 1704 PetscValidHeaderSpecific(*fact,MAT_COOKIE); 1705 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1706 if (mat->M != (*fact)->M || mat->N != (*fact)->N) { 1707 SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d", 1708 mat->M,(*fact)->M,mat->N,(*fact)->N); 1709 } 1710 if (!(*fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1711 1712 ierr = PetscLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1713 ierr = (*(*fact)->ops->lufactornumeric)(mat,fact);CHKERRQ(ierr); 1714 ierr = PetscLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1715 1716 ierr = MatView_Private(*fact);CHKERRQ(ierr); 1717 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1718 PetscFunctionReturn(0); 1719 } 1720 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "MatCholeskyFactor" 1723 /*@ 1724 MatCholeskyFactor - Performs in-place Cholesky factorization of a 1725 symmetric matrix. 1726 1727 Collective on Mat 1728 1729 Input Parameters: 1730 + mat - the matrix 1731 . perm - row and column permutations 1732 - f - expected fill as ratio of original fill 1733 1734 Notes: 1735 See MatLUFactor() for the nonsymmetric case. See also 1736 MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric(). 1737 1738 Most users should employ the simplified KSP interface for linear solvers 1739 instead of working directly with matrix algebra routines such as this. 1740 See, e.g., KSPCreate(). 1741 1742 Level: developer 1743 1744 Concepts: matrices^Cholesky factorization 1745 1746 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric() 1747 MatGetOrdering() 1748 1749 @*/ 1750 int MatCholeskyFactor(Mat mat,IS perm,MatFactorInfo *info) 1751 { 1752 int ierr; 1753 1754 PetscFunctionBegin; 1755 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1756 PetscValidType(mat); 1757 MatPreallocated(mat); 1758 PetscValidHeaderSpecific(perm,IS_COOKIE); 1759 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 1760 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1761 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1762 if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1763 1764 ierr = PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);CHKERRQ(ierr); 1765 ierr = (*mat->ops->choleskyfactor)(mat,perm,info);CHKERRQ(ierr); 1766 ierr = PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);CHKERRQ(ierr); 1767 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 1768 PetscFunctionReturn(0); 1769 } 1770 1771 #undef __FUNCT__ 1772 #define __FUNCT__ "MatCholeskyFactorSymbolic" 1773 /*@ 1774 MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization 1775 of a symmetric matrix. 1776 1777 Collective on Mat 1778 1779 Input Parameters: 1780 + mat - the matrix 1781 . perm - row and column permutations 1782 - info - options for factorization, includes 1783 $ fill - expected fill as ratio of original fill. 1784 $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting) 1785 $ Run with the option -log_info to determine an optimal value to use 1786 1787 Output Parameter: 1788 . fact - the factored matrix 1789 1790 Notes: 1791 See MatLUFactorSymbolic() for the nonsymmetric case. See also 1792 MatCholeskyFactor() and MatCholeskyFactorNumeric(). 1793 1794 Most users should employ the simplified KSP interface for linear solvers 1795 instead of working directly with matrix algebra routines such as this. 1796 See, e.g., KSPCreate(). 1797 1798 Level: developer 1799 1800 Concepts: matrices^Cholesky symbolic factorization 1801 1802 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric() 1803 MatGetOrdering() 1804 1805 @*/ 1806 int MatCholeskyFactorSymbolic(Mat mat,IS perm,MatFactorInfo *info,Mat *fact) 1807 { 1808 int ierr; 1809 1810 PetscFunctionBegin; 1811 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1812 PetscValidType(mat); 1813 MatPreallocated(mat); 1814 PetscValidPointer(fact); 1815 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 1816 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1817 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1818 if (!mat->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1819 1820 ierr = PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 1821 ierr = (*mat->ops->choleskyfactorsymbolic)(mat,perm,info,fact);CHKERRQ(ierr); 1822 ierr = PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 1823 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1824 PetscFunctionReturn(0); 1825 } 1826 1827 #undef __FUNCT__ 1828 #define __FUNCT__ "MatCholeskyFactorNumeric" 1829 /*@ 1830 MatCholeskyFactorNumeric - Performs numeric Cholesky factorization 1831 of a symmetric matrix. Call this routine after first calling 1832 MatCholeskyFactorSymbolic(). 1833 1834 Collective on Mat 1835 1836 Input Parameter: 1837 . mat - the initial matrix 1838 1839 Output Parameter: 1840 . fact - the factored matrix 1841 1842 Notes: 1843 Most users should employ the simplified KSP interface for linear solvers 1844 instead of working directly with matrix algebra routines such as this. 1845 See, e.g., KSPCreate(). 1846 1847 Level: developer 1848 1849 Concepts: matrices^Cholesky numeric factorization 1850 1851 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric() 1852 @*/ 1853 int MatCholeskyFactorNumeric(Mat mat,Mat *fact) 1854 { 1855 int ierr; 1856 1857 PetscFunctionBegin; 1858 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1859 PetscValidType(mat); 1860 MatPreallocated(mat); 1861 PetscValidPointer(fact); 1862 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1863 if (!(*fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1864 if (mat->M != (*fact)->M || mat->N != (*fact)->N) { 1865 SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d", 1866 mat->M,(*fact)->M,mat->N,(*fact)->N); 1867 } 1868 1869 ierr = PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1870 ierr = (*(*fact)->ops->choleskyfactornumeric)(mat,fact);CHKERRQ(ierr); 1871 ierr = PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1872 ierr = PetscObjectIncreaseState((PetscObject)*fact); CHKERRQ(ierr); 1873 PetscFunctionReturn(0); 1874 } 1875 1876 /* ----------------------------------------------------------------*/ 1877 #undef __FUNCT__ 1878 #define __FUNCT__ "MatSolve" 1879 /*@ 1880 MatSolve - Solves A x = b, given a factored matrix. 1881 1882 Collective on Mat and Vec 1883 1884 Input Parameters: 1885 + mat - the factored matrix 1886 - b - the right-hand-side vector 1887 1888 Output Parameter: 1889 . x - the result vector 1890 1891 Notes: 1892 The vectors b and x cannot be the same. I.e., one cannot 1893 call MatSolve(A,x,x). 1894 1895 Notes: 1896 Most users should employ the simplified KSP interface for linear solvers 1897 instead of working directly with matrix algebra routines such as this. 1898 See, e.g., KSPCreate(). 1899 1900 Level: developer 1901 1902 Concepts: matrices^triangular solves 1903 1904 .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd() 1905 @*/ 1906 int MatSolve(Mat mat,Vec b,Vec x) 1907 { 1908 int ierr; 1909 1910 PetscFunctionBegin; 1911 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1912 PetscValidType(mat); 1913 MatPreallocated(mat); 1914 PetscValidHeaderSpecific(b,VEC_COOKIE); 1915 PetscValidHeaderSpecific(x,VEC_COOKIE); 1916 PetscCheckSameComm(mat,b); 1917 PetscCheckSameComm(mat,x); 1918 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 1919 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 1920 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1921 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1922 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1923 if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0); 1924 1925 if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1926 ierr = PetscLogEventBegin(MAT_Solve,mat,b,x,0);CHKERRQ(ierr); 1927 ierr = (*mat->ops->solve)(mat,b,x);CHKERRQ(ierr); 1928 ierr = PetscLogEventEnd(MAT_Solve,mat,b,x,0);CHKERRQ(ierr); 1929 ierr = PetscObjectIncreaseState((PetscObject)x); CHKERRQ(ierr); 1930 PetscFunctionReturn(0); 1931 } 1932 1933 #undef __FUNCT__ 1934 #define __FUNCT__ "MatForwardSolve" 1935 /* @ 1936 MatForwardSolve - Solves L x = b, given a factored matrix, A = LU. 1937 1938 Collective on Mat and Vec 1939 1940 Input Parameters: 1941 + mat - the factored matrix 1942 - b - the right-hand-side vector 1943 1944 Output Parameter: 1945 . x - the result vector 1946 1947 Notes: 1948 MatSolve() should be used for most applications, as it performs 1949 a forward solve followed by a backward solve. 1950 1951 The vectors b and x cannot be the same. I.e., one cannot 1952 call MatForwardSolve(A,x,x). 1953 1954 Most users should employ the simplified KSP interface for linear solvers 1955 instead of working directly with matrix algebra routines such as this. 1956 See, e.g., KSPCreate(). 1957 1958 Level: developer 1959 1960 Concepts: matrices^forward solves 1961 1962 .seealso: MatSolve(), MatBackwardSolve() 1963 @ */ 1964 int MatForwardSolve(Mat mat,Vec b,Vec x) 1965 { 1966 int ierr; 1967 1968 PetscFunctionBegin; 1969 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1970 PetscValidType(mat); 1971 MatPreallocated(mat); 1972 PetscValidHeaderSpecific(b,VEC_COOKIE); 1973 PetscValidHeaderSpecific(x,VEC_COOKIE); 1974 PetscCheckSameComm(mat,b); 1975 PetscCheckSameComm(mat,x); 1976 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 1977 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 1978 if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1979 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1980 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1981 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1982 1983 ierr = PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);CHKERRQ(ierr); 1984 ierr = (*mat->ops->forwardsolve)(mat,b,x);CHKERRQ(ierr); 1985 ierr = PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);CHKERRQ(ierr); 1986 ierr = PetscObjectIncreaseState((PetscObject)x); CHKERRQ(ierr); 1987 PetscFunctionReturn(0); 1988 } 1989 1990 #undef __FUNCT__ 1991 #define __FUNCT__ "MatBackwardSolve" 1992 /* @ 1993 MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU. 1994 1995 Collective on Mat and Vec 1996 1997 Input Parameters: 1998 + mat - the factored matrix 1999 - b - the right-hand-side vector 2000 2001 Output Parameter: 2002 . x - the result vector 2003 2004 Notes: 2005 MatSolve() should be used for most applications, as it performs 2006 a forward solve followed by a backward solve. 2007 2008 The vectors b and x cannot be the same. I.e., one cannot 2009 call MatBackwardSolve(A,x,x). 2010 2011 Most users should employ the simplified KSP interface for linear solvers 2012 instead of working directly with matrix algebra routines such as this. 2013 See, e.g., KSPCreate(). 2014 2015 Level: developer 2016 2017 Concepts: matrices^backward solves 2018 2019 .seealso: MatSolve(), MatForwardSolve() 2020 @ */ 2021 int MatBackwardSolve(Mat mat,Vec b,Vec x) 2022 { 2023 int ierr; 2024 2025 PetscFunctionBegin; 2026 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2027 PetscValidType(mat); 2028 MatPreallocated(mat); 2029 PetscValidHeaderSpecific(b,VEC_COOKIE); 2030 PetscValidHeaderSpecific(x,VEC_COOKIE); 2031 PetscCheckSameComm(mat,b); 2032 PetscCheckSameComm(mat,x); 2033 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2034 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2035 if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2036 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2037 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2038 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2039 2040 ierr = PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);CHKERRQ(ierr); 2041 ierr = (*mat->ops->backwardsolve)(mat,b,x);CHKERRQ(ierr); 2042 ierr = PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);CHKERRQ(ierr); 2043 ierr = PetscObjectIncreaseState((PetscObject)x); CHKERRQ(ierr); 2044 PetscFunctionReturn(0); 2045 } 2046 2047 #undef __FUNCT__ 2048 #define __FUNCT__ "MatSolveAdd" 2049 /*@ 2050 MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix. 2051 2052 Collective on Mat and Vec 2053 2054 Input Parameters: 2055 + mat - the factored matrix 2056 . b - the right-hand-side vector 2057 - y - the vector to be added to 2058 2059 Output Parameter: 2060 . x - the result vector 2061 2062 Notes: 2063 The vectors b and x cannot be the same. I.e., one cannot 2064 call MatSolveAdd(A,x,y,x). 2065 2066 Most users should employ the simplified KSP interface for linear solvers 2067 instead of working directly with matrix algebra routines such as this. 2068 See, e.g., KSPCreate(). 2069 2070 Level: developer 2071 2072 Concepts: matrices^triangular solves 2073 2074 .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd() 2075 @*/ 2076 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x) 2077 { 2078 PetscScalar one = 1.0; 2079 Vec tmp; 2080 int ierr; 2081 2082 PetscFunctionBegin; 2083 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2084 PetscValidType(mat); 2085 MatPreallocated(mat); 2086 PetscValidHeaderSpecific(y,VEC_COOKIE); 2087 PetscValidHeaderSpecific(b,VEC_COOKIE); 2088 PetscValidHeaderSpecific(x,VEC_COOKIE); 2089 PetscCheckSameComm(mat,b); 2090 PetscCheckSameComm(mat,y); 2091 PetscCheckSameComm(mat,x); 2092 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2093 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2094 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2095 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2096 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 2097 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2098 if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n); 2099 2100 ierr = PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);CHKERRQ(ierr); 2101 if (mat->ops->solveadd) { 2102 ierr = (*mat->ops->solveadd)(mat,b,y,x);CHKERRQ(ierr); 2103 } else { 2104 /* do the solve then the add manually */ 2105 if (x != y) { 2106 ierr = MatSolve(mat,b,x);CHKERRQ(ierr); 2107 ierr = VecAXPY(&one,y,x);CHKERRQ(ierr); 2108 } else { 2109 ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr); 2110 PetscLogObjectParent(mat,tmp); 2111 ierr = VecCopy(x,tmp);CHKERRQ(ierr); 2112 ierr = MatSolve(mat,b,x);CHKERRQ(ierr); 2113 ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr); 2114 ierr = VecDestroy(tmp);CHKERRQ(ierr); 2115 } 2116 } 2117 ierr = PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);CHKERRQ(ierr); 2118 ierr = PetscObjectIncreaseState((PetscObject)x); CHKERRQ(ierr); 2119 PetscFunctionReturn(0); 2120 } 2121 2122 #undef __FUNCT__ 2123 #define __FUNCT__ "MatSolveTranspose" 2124 /*@ 2125 MatSolveTranspose - Solves A' x = b, given a factored matrix. 2126 2127 Collective on Mat and Vec 2128 2129 Input Parameters: 2130 + mat - the factored matrix 2131 - b - the right-hand-side vector 2132 2133 Output Parameter: 2134 . x - the result vector 2135 2136 Notes: 2137 The vectors b and x cannot be the same. I.e., one cannot 2138 call MatSolveTranspose(A,x,x). 2139 2140 Most users should employ the simplified KSP interface for linear solvers 2141 instead of working directly with matrix algebra routines such as this. 2142 See, e.g., KSPCreate(). 2143 2144 Level: developer 2145 2146 Concepts: matrices^triangular solves 2147 2148 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd() 2149 @*/ 2150 int MatSolveTranspose(Mat mat,Vec b,Vec x) 2151 { 2152 int ierr; 2153 2154 PetscFunctionBegin; 2155 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2156 PetscValidType(mat); 2157 MatPreallocated(mat); 2158 PetscValidHeaderSpecific(b,VEC_COOKIE); 2159 PetscValidHeaderSpecific(x,VEC_COOKIE); 2160 PetscCheckSameComm(mat,b); 2161 PetscCheckSameComm(mat,x); 2162 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2163 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2164 if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",mat->type_name); 2165 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 2166 if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N); 2167 2168 ierr = PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);CHKERRQ(ierr); 2169 ierr = (*mat->ops->solvetranspose)(mat,b,x);CHKERRQ(ierr); 2170 ierr = PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);CHKERRQ(ierr); 2171 ierr = PetscObjectIncreaseState((PetscObject)x); CHKERRQ(ierr); 2172 PetscFunctionReturn(0); 2173 } 2174 2175 #undef __FUNCT__ 2176 #define __FUNCT__ "MatSolveTransposeAdd" 2177 /*@ 2178 MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a 2179 factored matrix. 2180 2181 Collective on Mat and Vec 2182 2183 Input Parameters: 2184 + mat - the factored matrix 2185 . b - the right-hand-side vector 2186 - y - the vector to be added to 2187 2188 Output Parameter: 2189 . x - the result vector 2190 2191 Notes: 2192 The vectors b and x cannot be the same. I.e., one cannot 2193 call MatSolveTransposeAdd(A,x,y,x). 2194 2195 Most users should employ the simplified KSP interface for linear solvers 2196 instead of working directly with matrix algebra routines such as this. 2197 See, e.g., KSPCreate(). 2198 2199 Level: developer 2200 2201 Concepts: matrices^triangular solves 2202 2203 .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose() 2204 @*/ 2205 int MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x) 2206 { 2207 PetscScalar one = 1.0; 2208 int ierr; 2209 Vec tmp; 2210 2211 PetscFunctionBegin; 2212 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2213 PetscValidType(mat); 2214 MatPreallocated(mat); 2215 PetscValidHeaderSpecific(y,VEC_COOKIE); 2216 PetscValidHeaderSpecific(b,VEC_COOKIE); 2217 PetscValidHeaderSpecific(x,VEC_COOKIE); 2218 PetscCheckSameComm(mat,b); 2219 PetscCheckSameComm(mat,y); 2220 PetscCheckSameComm(mat,x); 2221 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2222 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2223 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 2224 if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N); 2225 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N); 2226 if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n); 2227 2228 ierr = PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);CHKERRQ(ierr); 2229 if (mat->ops->solvetransposeadd) { 2230 ierr = (*mat->ops->solvetransposeadd)(mat,b,y,x);CHKERRQ(ierr); 2231 } else { 2232 /* do the solve then the add manually */ 2233 if (x != y) { 2234 ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr); 2235 ierr = VecAXPY(&one,y,x);CHKERRQ(ierr); 2236 } else { 2237 ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr); 2238 PetscLogObjectParent(mat,tmp); 2239 ierr = VecCopy(x,tmp);CHKERRQ(ierr); 2240 ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr); 2241 ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr); 2242 ierr = VecDestroy(tmp);CHKERRQ(ierr); 2243 } 2244 } 2245 ierr = PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);CHKERRQ(ierr); 2246 ierr = PetscObjectIncreaseState((PetscObject)x); CHKERRQ(ierr); 2247 PetscFunctionReturn(0); 2248 } 2249 /* ----------------------------------------------------------------*/ 2250 2251 #undef __FUNCT__ 2252 #define __FUNCT__ "MatRelax" 2253 /*@ 2254 MatRelax - Computes one relaxation sweep. 2255 2256 Collective on Mat and Vec 2257 2258 Input Parameters: 2259 + mat - the matrix 2260 . b - the right hand side 2261 . omega - the relaxation factor 2262 . flag - flag indicating the type of SOR (see below) 2263 . shift - diagonal shift 2264 - its - the number of iterations 2265 - lits - the number of local iterations 2266 2267 Output Parameters: 2268 . x - the solution (can contain an initial guess) 2269 2270 SOR Flags: 2271 . SOR_FORWARD_SWEEP - forward SOR 2272 . SOR_BACKWARD_SWEEP - backward SOR 2273 . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR) 2274 . SOR_LOCAL_FORWARD_SWEEP - local forward SOR 2275 . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR 2276 . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR 2277 . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies 2278 upper/lower triangular part of matrix to 2279 vector (with omega) 2280 . SOR_ZERO_INITIAL_GUESS - zero initial guess 2281 2282 Notes: 2283 SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and 2284 SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings 2285 on each processor. 2286 2287 Application programmers will not generally use MatRelax() directly, 2288 but instead will employ the KSP/PC interface. 2289 2290 Notes for Advanced Users: 2291 The flags are implemented as bitwise inclusive or operations. 2292 For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP) 2293 to specify a zero initial guess for SSOR. 2294 2295 Most users should employ the simplified KSP interface for linear solvers 2296 instead of working directly with matrix algebra routines such as this. 2297 See, e.g., KSPCreate(). 2298 2299 Level: developer 2300 2301 Concepts: matrices^relaxation 2302 Concepts: matrices^SOR 2303 Concepts: matrices^Gauss-Seidel 2304 2305 @*/ 2306 int MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,int its,int lits,Vec x) 2307 { 2308 int ierr; 2309 2310 PetscFunctionBegin; 2311 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2312 PetscValidType(mat); 2313 MatPreallocated(mat); 2314 PetscValidHeaderSpecific(b,VEC_COOKIE); 2315 PetscValidHeaderSpecific(x,VEC_COOKIE); 2316 PetscCheckSameComm(mat,b); 2317 PetscCheckSameComm(mat,x); 2318 if (!mat->ops->relax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2319 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2320 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2321 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2322 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2323 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2324 2325 ierr = PetscLogEventBegin(MAT_Relax,mat,b,x,0);CHKERRQ(ierr); 2326 ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);CHKERRQ(ierr); 2327 ierr = PetscLogEventEnd(MAT_Relax,mat,b,x,0);CHKERRQ(ierr); 2328 ierr = PetscObjectIncreaseState((PetscObject)x); CHKERRQ(ierr); 2329 PetscFunctionReturn(0); 2330 } 2331 2332 #undef __FUNCT__ 2333 #define __FUNCT__ "MatCopy_Basic" 2334 /* 2335 Default matrix copy routine. 2336 */ 2337 int MatCopy_Basic(Mat A,Mat B,MatStructure str) 2338 { 2339 int ierr,i,rstart,rend,nz,*cwork; 2340 PetscScalar *vwork; 2341 2342 PetscFunctionBegin; 2343 ierr = MatZeroEntries(B);CHKERRQ(ierr); 2344 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 2345 for (i=rstart; i<rend; i++) { 2346 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2347 ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr); 2348 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2349 } 2350 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2351 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2352 ierr = PetscObjectIncreaseState((PetscObject)B); CHKERRQ(ierr); 2353 PetscFunctionReturn(0); 2354 } 2355 2356 #undef __FUNCT__ 2357 #define __FUNCT__ "MatCopy" 2358 /*@C 2359 MatCopy - Copys a matrix to another matrix. 2360 2361 Collective on Mat 2362 2363 Input Parameters: 2364 + A - the matrix 2365 - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN 2366 2367 Output Parameter: 2368 . B - where the copy is put 2369 2370 Notes: 2371 If you use SAME_NONZERO_PATTERN then the two matrices had better have the 2372 same nonzero pattern or the routine will crash. 2373 2374 MatCopy() copies the matrix entries of a matrix to another existing 2375 matrix (after first zeroing the second matrix). A related routine is 2376 MatConvert(), which first creates a new matrix and then copies the data. 2377 2378 Level: intermediate 2379 2380 Concepts: matrices^copying 2381 2382 .seealso: MatConvert(), MatDuplicate() 2383 2384 @*/ 2385 int MatCopy(Mat A,Mat B,MatStructure str) 2386 { 2387 int ierr; 2388 2389 PetscFunctionBegin; 2390 PetscValidHeaderSpecific(A,MAT_COOKIE); 2391 PetscValidHeaderSpecific(B,MAT_COOKIE); 2392 PetscValidType(A); 2393 MatPreallocated(A); 2394 PetscValidType(B); 2395 MatPreallocated(B); 2396 PetscCheckSameComm(A,B); 2397 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2398 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2399 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, 2400 A->N,B->N); 2401 2402 ierr = PetscLogEventBegin(MAT_Copy,A,B,0,0);CHKERRQ(ierr); 2403 if (A->ops->copy) { 2404 ierr = (*A->ops->copy)(A,B,str);CHKERRQ(ierr); 2405 } else { /* generic conversion */ 2406 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2407 } 2408 if (A->mapping) { 2409 if (B->mapping) {ierr = ISLocalToGlobalMappingDestroy(B->mapping);CHKERRQ(ierr);B->mapping = 0;} 2410 ierr = MatSetLocalToGlobalMapping(B,A->mapping);CHKERRQ(ierr); 2411 } 2412 if (A->bmapping) { 2413 if (B->bmapping) {ierr = ISLocalToGlobalMappingDestroy(B->bmapping);CHKERRQ(ierr);B->bmapping = 0;} 2414 ierr = MatSetLocalToGlobalMappingBlock(B,A->mapping);CHKERRQ(ierr); 2415 } 2416 ierr = PetscLogEventEnd(MAT_Copy,A,B,0,0);CHKERRQ(ierr); 2417 ierr = PetscObjectIncreaseState((PetscObject)B); CHKERRQ(ierr); 2418 PetscFunctionReturn(0); 2419 } 2420 2421 #include "petscsys.h" 2422 PetscTruth MatConvertRegisterAllCalled = PETSC_FALSE; 2423 PetscFList MatConvertList = 0; 2424 2425 #undef __FUNCT__ 2426 #define __FUNCT__ "MatConvertRegister" 2427 /*@C 2428 MatConvertRegister - Allows one to register a routine that converts a sparse matrix 2429 from one format to another. 2430 2431 Not Collective 2432 2433 Input Parameters: 2434 + type - the type of matrix (defined in include/petscmat.h), for example, MATSEQAIJ. 2435 - Converter - the function that reads the matrix from the binary file. 2436 2437 Level: developer 2438 2439 .seealso: MatConvertRegisterAll(), MatConvert() 2440 2441 @*/ 2442 int MatConvertRegister(char *sname,char *path,char *name,int (*function)(Mat,MatType,Mat*)) 2443 { 2444 int ierr; 2445 char fullname[PETSC_MAX_PATH_LEN]; 2446 2447 PetscFunctionBegin; 2448 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 2449 ierr = PetscFListAdd(&MatConvertList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 2450 PetscFunctionReturn(0); 2451 } 2452 2453 #undef __FUNCT__ 2454 #define __FUNCT__ "MatConvert" 2455 /*@C 2456 MatConvert - Converts a matrix to another matrix, either of the same 2457 or different type. 2458 2459 Collective on Mat 2460 2461 Input Parameters: 2462 + mat - the matrix 2463 - newtype - new matrix type. Use MATSAME to create a new matrix of the 2464 same type as the original matrix. 2465 2466 Output Parameter: 2467 . M - pointer to place new matrix 2468 2469 Notes: 2470 MatConvert() first creates a new matrix and then copies the data from 2471 the first matrix. A related routine is MatCopy(), which copies the matrix 2472 entries of one matrix to another already existing matrix context. 2473 2474 Level: intermediate 2475 2476 Concepts: matrices^converting between storage formats 2477 2478 .seealso: MatCopy(), MatDuplicate() 2479 @*/ 2480 int MatConvert(Mat mat,MatType newtype,Mat *M) 2481 { 2482 int ierr; 2483 PetscTruth sametype,issame,flg; 2484 char convname[256],mtype[256]; 2485 2486 PetscFunctionBegin; 2487 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2488 PetscValidType(mat); 2489 MatPreallocated(mat); 2490 PetscValidPointer(M); 2491 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2492 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2493 2494 ierr = PetscOptionsGetString(PETSC_NULL,"-matconvert_type",mtype,256,&flg);CHKERRQ(ierr); 2495 if (flg) { 2496 newtype = mtype; 2497 } 2498 ierr = PetscLogEventBegin(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2499 2500 ierr = PetscTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr); 2501 ierr = PetscStrcmp(newtype,"same",&issame);CHKERRQ(ierr); 2502 if ((sametype || issame) && mat->ops->duplicate) { 2503 ierr = (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);CHKERRQ(ierr); 2504 } else { 2505 int (*conv)(Mat,MatType,Mat*)=PETSC_NULL; 2506 /* 2507 Order of precedence: 2508 1) See if a specialized converter is known to the current matrix. 2509 2) See if a specialized converter is known to the desired matrix class. 2510 3) See if a good general converter is registered for the desired class 2511 (as of 6/27/03 only MATMPIADJ falls into this category). 2512 4) See if a good general converter is known for the current matrix. 2513 5) Use a really basic converter. 2514 */ 2515 ierr = PetscStrcpy(convname,"MatConvert_");CHKERRQ(ierr); 2516 ierr = PetscStrcat(convname,mat->type_name);CHKERRQ(ierr); 2517 ierr = PetscStrcat(convname,"_");CHKERRQ(ierr); 2518 ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr); 2519 ierr = PetscStrcat(convname,"_C");CHKERRQ(ierr); 2520 ierr = PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2521 if (!conv) { 2522 Mat B; 2523 ierr = MatCreate(mat->comm,0,0,0,0,&B);CHKERRQ(ierr); 2524 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 2525 ierr = PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2526 ierr = MatDestroy(B);CHKERRQ(ierr); 2527 if (!conv) { 2528 if (!MatConvertRegisterAllCalled) { 2529 ierr = MatConvertRegisterAll(PETSC_NULL);CHKERRQ(ierr); 2530 } 2531 ierr = PetscFListFind(mat->comm,MatConvertList,newtype,(void(**)(void))&conv);CHKERRQ(ierr); 2532 if (!conv) { 2533 if (mat->ops->convert) { 2534 conv = mat->ops->convert; 2535 } else { 2536 conv = MatConvert_Basic; 2537 } 2538 } 2539 } 2540 } 2541 ierr = (*conv)(mat,newtype,M);CHKERRQ(ierr); 2542 } 2543 ierr = PetscLogEventEnd(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2544 ierr = PetscObjectIncreaseState((PetscObject)*M); CHKERRQ(ierr); 2545 PetscFunctionReturn(0); 2546 } 2547 2548 2549 #undef __FUNCT__ 2550 #define __FUNCT__ "MatDuplicate" 2551 /*@C 2552 MatDuplicate - Duplicates a matrix including the non-zero structure. 2553 2554 Collective on Mat 2555 2556 Input Parameters: 2557 + mat - the matrix 2558 - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero 2559 values as well or not 2560 2561 Output Parameter: 2562 . M - pointer to place new matrix 2563 2564 Level: intermediate 2565 2566 Concepts: matrices^duplicating 2567 2568 .seealso: MatCopy(), MatConvert() 2569 @*/ 2570 int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M) 2571 { 2572 int ierr; 2573 2574 PetscFunctionBegin; 2575 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2576 PetscValidType(mat); 2577 MatPreallocated(mat); 2578 PetscValidPointer(M); 2579 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2580 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2581 2582 *M = 0; 2583 ierr = PetscLogEventBegin(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2584 if (!mat->ops->duplicate) { 2585 SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type"); 2586 } 2587 ierr = (*mat->ops->duplicate)(mat,op,M);CHKERRQ(ierr); 2588 if (mat->mapping) { 2589 ierr = MatSetLocalToGlobalMapping(*M,mat->mapping);CHKERRQ(ierr); 2590 } 2591 if (mat->bmapping) { 2592 ierr = MatSetLocalToGlobalMappingBlock(*M,mat->mapping);CHKERRQ(ierr); 2593 } 2594 ierr = PetscLogEventEnd(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2595 ierr = PetscObjectIncreaseState((PetscObject)*M); CHKERRQ(ierr); 2596 PetscFunctionReturn(0); 2597 } 2598 2599 #undef __FUNCT__ 2600 #define __FUNCT__ "MatGetDiagonal" 2601 /*@ 2602 MatGetDiagonal - Gets the diagonal of a matrix. 2603 2604 Collective on Mat and Vec 2605 2606 Input Parameters: 2607 + mat - the matrix 2608 - v - the vector for storing the diagonal 2609 2610 Output Parameter: 2611 . v - the diagonal of the matrix 2612 2613 Notes: 2614 For the SeqAIJ matrix format, this routine may also be called 2615 on a LU factored matrix; in that case it routines the reciprocal of 2616 the diagonal entries in U. It returns the entries permuted by the 2617 row and column permutation used during the symbolic factorization. 2618 2619 Level: intermediate 2620 2621 Concepts: matrices^accessing diagonals 2622 2623 .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix(), MatGetRowMax() 2624 @*/ 2625 int MatGetDiagonal(Mat mat,Vec v) 2626 { 2627 int ierr; 2628 2629 PetscFunctionBegin; 2630 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2631 PetscValidType(mat); 2632 MatPreallocated(mat); 2633 PetscValidHeaderSpecific(v,VEC_COOKIE); 2634 /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */ 2635 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2636 if (!mat->ops->getdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2637 2638 ierr = (*mat->ops->getdiagonal)(mat,v);CHKERRQ(ierr); 2639 ierr = PetscObjectIncreaseState((PetscObject)v); CHKERRQ(ierr); 2640 PetscFunctionReturn(0); 2641 } 2642 2643 #undef __FUNCT__ 2644 #define __FUNCT__ "MatGetRowMax" 2645 /*@ 2646 MatGetRowMax - Gets the maximum value (in absolute value) of each 2647 row of the matrix 2648 2649 Collective on Mat and Vec 2650 2651 Input Parameters: 2652 . mat - the matrix 2653 2654 Output Parameter: 2655 . v - the vector for storing the maximums 2656 2657 Level: intermediate 2658 2659 Concepts: matrices^getting row maximums 2660 2661 .seealso: MatGetDiagonal(), MatGetSubmatrices(), MatGetSubmatrix() 2662 @*/ 2663 int MatGetRowMax(Mat mat,Vec v) 2664 { 2665 int ierr; 2666 2667 PetscFunctionBegin; 2668 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2669 PetscValidType(mat); 2670 MatPreallocated(mat); 2671 PetscValidHeaderSpecific(v,VEC_COOKIE); 2672 /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */ 2673 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2674 if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2675 2676 ierr = (*mat->ops->getrowmax)(mat,v);CHKERRQ(ierr); 2677 ierr = PetscObjectIncreaseState((PetscObject)v); CHKERRQ(ierr); 2678 PetscFunctionReturn(0); 2679 } 2680 2681 #undef __FUNCT__ 2682 #define __FUNCT__ "MatTranspose" 2683 /*@C 2684 MatTranspose - Computes an in-place or out-of-place transpose of a matrix. 2685 2686 Collective on Mat 2687 2688 Input Parameter: 2689 . mat - the matrix to transpose 2690 2691 Output Parameters: 2692 . B - the transpose (or pass in PETSC_NULL for an in-place transpose) 2693 2694 Level: intermediate 2695 2696 Concepts: matrices^transposing 2697 2698 .seealso: MatMultTranspose(), MatMultTransposeAdd(), MatIsTranspose() 2699 @*/ 2700 int MatTranspose(Mat mat,Mat *B) 2701 { 2702 int ierr; 2703 2704 PetscFunctionBegin; 2705 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2706 PetscValidType(mat); 2707 MatPreallocated(mat); 2708 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2709 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2710 if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2711 2712 ierr = PetscLogEventBegin(MAT_Transpose,mat,0,0,0);CHKERRQ(ierr); 2713 ierr = (*mat->ops->transpose)(mat,B);CHKERRQ(ierr); 2714 ierr = PetscLogEventEnd(MAT_Transpose,mat,0,0,0);CHKERRQ(ierr); 2715 if (B) {ierr = PetscObjectIncreaseState((PetscObject)*B); CHKERRQ(ierr);} 2716 PetscFunctionReturn(0); 2717 } 2718 2719 #undef __FUNCT__ 2720 #define __FUNCT__ "MatIsTranspose" 2721 /*@C 2722 MatIsTranspose - Test whether a matrix is another one's transpose, 2723 or its own, in which case it tests symmetry. 2724 2725 Collective on Mat 2726 2727 Input Parameter: 2728 + A - the matrix to test 2729 - B - the matrix to test against, this can equal the first parameter 2730 2731 Output Parameters: 2732 . flg - the result 2733 2734 Notes: 2735 Only available for SeqAIJ/MPIAIJ matrices. The sequential algorithm 2736 has a running time of the order of the number of nonzeros; the parallel 2737 test involves parallel copies of the block-offdiagonal parts of the matrix. 2738 2739 Level: intermediate 2740 2741 Concepts: matrices^transposing, matrix^symmetry 2742 2743 .seealso: MatTranspose(), MatIsSymmetric(), MatIsHermitian() 2744 @*/ 2745 int MatIsTranspose(Mat A,Mat B,PetscTruth *flg) 2746 { 2747 int ierr,(*f)(Mat,Mat,PetscTruth*),(*g)(Mat,Mat,PetscTruth*); 2748 2749 PetscFunctionBegin; 2750 PetscValidHeaderSpecific(A,MAT_COOKIE); 2751 PetscValidHeaderSpecific(B,MAT_COOKIE); 2752 ierr = PetscObjectQueryFunction((PetscObject)A,"MatIsTranspose_C",(void (**)(void))&f);CHKERRQ(ierr); 2753 ierr = PetscObjectQueryFunction((PetscObject)B,"MatIsTranspose_C",(void (**)(void))&g);CHKERRQ(ierr); 2754 if (f && g) { 2755 if (f==g) { 2756 ierr = (*f)(A,B,flg);CHKERRQ(ierr); 2757 } else { 2758 SETERRQ(1,"Matrices do not have the same comparator for symmetry test"); 2759 } 2760 } 2761 PetscFunctionReturn(0); 2762 } 2763 2764 #undef __FUNCT__ 2765 #define __FUNCT__ "MatPermute" 2766 /*@C 2767 MatPermute - Creates a new matrix with rows and columns permuted from the 2768 original. 2769 2770 Collective on Mat 2771 2772 Input Parameters: 2773 + mat - the matrix to permute 2774 . row - row permutation, each processor supplies only the permutation for its rows 2775 - col - column permutation, each processor needs the entire column permutation, that is 2776 this is the same size as the total number of columns in the matrix 2777 2778 Output Parameters: 2779 . B - the permuted matrix 2780 2781 Level: advanced 2782 2783 Concepts: matrices^permuting 2784 2785 .seealso: MatGetOrdering() 2786 @*/ 2787 int MatPermute(Mat mat,IS row,IS col,Mat *B) 2788 { 2789 int ierr; 2790 2791 PetscFunctionBegin; 2792 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2793 PetscValidType(mat); 2794 MatPreallocated(mat); 2795 PetscValidHeaderSpecific(row,IS_COOKIE); 2796 PetscValidHeaderSpecific(col,IS_COOKIE); 2797 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2798 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2799 if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2800 ierr = (*mat->ops->permute)(mat,row,col,B);CHKERRQ(ierr); 2801 ierr = PetscObjectIncreaseState((PetscObject)*B); CHKERRQ(ierr); 2802 PetscFunctionReturn(0); 2803 } 2804 2805 #undef __FUNCT__ 2806 #define __FUNCT__ "MatPermuteSparsify" 2807 /*@C 2808 MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the 2809 original and sparsified to the prescribed tolerance. 2810 2811 Collective on Mat 2812 2813 Input Parameters: 2814 + A - The matrix to permute 2815 . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE 2816 . frac - The half-bandwidth as a fraction of the total size, or 0.0 2817 . tol - The drop tolerance 2818 . rowp - The row permutation 2819 - colp - The column permutation 2820 2821 Output Parameter: 2822 . B - The permuted, sparsified matrix 2823 2824 Level: advanced 2825 2826 Note: 2827 The default behavior (band = PETSC_DECIDE and frac = 0.0) is to 2828 restrict the half-bandwidth of the resulting matrix to 5% of the 2829 total matrix size. 2830 2831 .keywords: matrix, permute, sparsify 2832 2833 .seealso: MatGetOrdering(), MatPermute() 2834 @*/ 2835 int MatPermuteSparsify(Mat A, int band, PetscReal frac, PetscReal tol, IS rowp, IS colp, Mat *B) 2836 { 2837 IS irowp, icolp; 2838 int *rows, *cols; 2839 int M, N, locRowStart, locRowEnd; 2840 int nz, newNz; 2841 int *cwork, *cnew; 2842 PetscScalar *vwork, *vnew; 2843 int bw, size; 2844 int row, locRow, newRow, col, newCol; 2845 int ierr; 2846 2847 PetscFunctionBegin; 2848 PetscValidHeaderSpecific(A, MAT_COOKIE); 2849 PetscValidHeaderSpecific(rowp, IS_COOKIE); 2850 PetscValidHeaderSpecific(colp, IS_COOKIE); 2851 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 2852 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 2853 if (!A->ops->permutesparsify) { 2854 ierr = MatGetSize(A, &M, &N); CHKERRQ(ierr); 2855 ierr = MatGetOwnershipRange(A, &locRowStart, &locRowEnd); CHKERRQ(ierr); 2856 ierr = ISGetSize(rowp, &size); CHKERRQ(ierr); 2857 if (size != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for row permutation, should be %d", size, M); 2858 ierr = ISGetSize(colp, &size); CHKERRQ(ierr); 2859 if (size != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for column permutation, should be %d", size, N); 2860 ierr = ISInvertPermutation(rowp, 0, &irowp); CHKERRQ(ierr); 2861 ierr = ISGetIndices(irowp, &rows); CHKERRQ(ierr); 2862 ierr = ISInvertPermutation(colp, 0, &icolp); CHKERRQ(ierr); 2863 ierr = ISGetIndices(icolp, &cols); CHKERRQ(ierr); 2864 ierr = PetscMalloc(N * sizeof(int), &cnew); CHKERRQ(ierr); 2865 ierr = PetscMalloc(N * sizeof(PetscScalar), &vnew); CHKERRQ(ierr); 2866 2867 /* Setup bandwidth to include */ 2868 if (band == PETSC_DECIDE) { 2869 if (frac <= 0.0) 2870 bw = (int) (M * 0.05); 2871 else 2872 bw = (int) (M * frac); 2873 } else { 2874 if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer"); 2875 bw = band; 2876 } 2877 2878 /* Put values into new matrix */ 2879 ierr = MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B); CHKERRQ(ierr); 2880 for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) { 2881 ierr = MatGetRow(A, row, &nz, &cwork, &vwork); CHKERRQ(ierr); 2882 newRow = rows[locRow]+locRowStart; 2883 for(col = 0, newNz = 0; col < nz; col++) { 2884 newCol = cols[cwork[col]]; 2885 if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) { 2886 cnew[newNz] = newCol; 2887 vnew[newNz] = vwork[col]; 2888 newNz++; 2889 } 2890 } 2891 ierr = MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES); CHKERRQ(ierr); 2892 ierr = MatRestoreRow(A, row, &nz, &cwork, &vwork); CHKERRQ(ierr); 2893 } 2894 ierr = PetscFree(cnew); CHKERRQ(ierr); 2895 ierr = PetscFree(vnew); CHKERRQ(ierr); 2896 ierr = MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2897 ierr = MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2898 ierr = ISRestoreIndices(irowp, &rows); CHKERRQ(ierr); 2899 ierr = ISRestoreIndices(icolp, &cols); CHKERRQ(ierr); 2900 ierr = ISDestroy(irowp); CHKERRQ(ierr); 2901 ierr = ISDestroy(icolp); CHKERRQ(ierr); 2902 } else { 2903 ierr = (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B); CHKERRQ(ierr); 2904 } 2905 ierr = PetscObjectIncreaseState((PetscObject)*B); CHKERRQ(ierr); 2906 PetscFunctionReturn(0); 2907 } 2908 2909 #undef __FUNCT__ 2910 #define __FUNCT__ "MatEqual" 2911 /*@ 2912 MatEqual - Compares two matrices. 2913 2914 Collective on Mat 2915 2916 Input Parameters: 2917 + A - the first matrix 2918 - B - the second matrix 2919 2920 Output Parameter: 2921 . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise. 2922 2923 Level: intermediate 2924 2925 Concepts: matrices^equality between 2926 @*/ 2927 int MatEqual(Mat A,Mat B,PetscTruth *flg) 2928 { 2929 int ierr; 2930 2931 PetscFunctionBegin; 2932 PetscValidHeaderSpecific(A,MAT_COOKIE); 2933 PetscValidHeaderSpecific(B,MAT_COOKIE); 2934 PetscValidType(A); 2935 MatPreallocated(A); 2936 PetscValidType(B); 2937 MatPreallocated(B); 2938 PetscValidIntPointer(flg); 2939 PetscCheckSameComm(A,B); 2940 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2941 if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2942 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); 2943 if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",A->type_name); 2944 if (!B->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",B->type_name); 2945 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); 2946 ierr = (*A->ops->equal)(A,B,flg);CHKERRQ(ierr); 2947 PetscFunctionReturn(0); 2948 } 2949 2950 #undef __FUNCT__ 2951 #define __FUNCT__ "MatDiagonalScale" 2952 /*@ 2953 MatDiagonalScale - Scales a matrix on the left and right by diagonal 2954 matrices that are stored as vectors. Either of the two scaling 2955 matrices can be PETSC_NULL. 2956 2957 Collective on Mat 2958 2959 Input Parameters: 2960 + mat - the matrix to be scaled 2961 . l - the left scaling vector (or PETSC_NULL) 2962 - r - the right scaling vector (or PETSC_NULL) 2963 2964 Notes: 2965 MatDiagonalScale() computes A = LAR, where 2966 L = a diagonal matrix, R = a diagonal matrix 2967 2968 Level: intermediate 2969 2970 Concepts: matrices^diagonal scaling 2971 Concepts: diagonal scaling of matrices 2972 2973 .seealso: MatScale() 2974 @*/ 2975 int MatDiagonalScale(Mat mat,Vec l,Vec r) 2976 { 2977 int ierr; 2978 2979 PetscFunctionBegin; 2980 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2981 PetscValidType(mat); 2982 MatPreallocated(mat); 2983 if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2984 if (l) {PetscValidHeaderSpecific(l,VEC_COOKIE);PetscCheckSameComm(mat,l);} 2985 if (r) {PetscValidHeaderSpecific(r,VEC_COOKIE);PetscCheckSameComm(mat,r);} 2986 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2987 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2988 2989 ierr = PetscLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 2990 ierr = (*mat->ops->diagonalscale)(mat,l,r);CHKERRQ(ierr); 2991 ierr = PetscLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 2992 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 2993 PetscFunctionReturn(0); 2994 } 2995 2996 #undef __FUNCT__ 2997 #define __FUNCT__ "MatScale" 2998 /*@ 2999 MatScale - Scales all elements of a matrix by a given number. 3000 3001 Collective on Mat 3002 3003 Input Parameters: 3004 + mat - the matrix to be scaled 3005 - a - the scaling value 3006 3007 Output Parameter: 3008 . mat - the scaled matrix 3009 3010 Level: intermediate 3011 3012 Concepts: matrices^scaling all entries 3013 3014 .seealso: MatDiagonalScale() 3015 @*/ 3016 int MatScale(const PetscScalar *a,Mat mat) 3017 { 3018 int ierr; 3019 3020 PetscFunctionBegin; 3021 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3022 PetscValidType(mat); 3023 MatPreallocated(mat); 3024 PetscValidScalarPointer(a); 3025 if (!mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3026 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3027 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3028 3029 ierr = PetscLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 3030 ierr = (*mat->ops->scale)(a,mat);CHKERRQ(ierr); 3031 ierr = PetscLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 3032 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3033 PetscFunctionReturn(0); 3034 } 3035 3036 #undef __FUNCT__ 3037 #define __FUNCT__ "MatNorm" 3038 /*@ 3039 MatNorm - Calculates various norms of a matrix. 3040 3041 Collective on Mat 3042 3043 Input Parameters: 3044 + mat - the matrix 3045 - type - the type of norm, NORM_1, NORM_FROBENIUS, NORM_INFINITY 3046 3047 Output Parameters: 3048 . nrm - the resulting norm 3049 3050 Level: intermediate 3051 3052 Concepts: matrices^norm 3053 Concepts: norm^of matrix 3054 @*/ 3055 int MatNorm(Mat mat,NormType type,PetscReal *nrm) 3056 { 3057 int ierr; 3058 3059 PetscFunctionBegin; 3060 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3061 PetscValidType(mat); 3062 MatPreallocated(mat); 3063 PetscValidScalarPointer(nrm); 3064 3065 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3066 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3067 if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3068 ierr = (*mat->ops->norm)(mat,type,nrm);CHKERRQ(ierr); 3069 PetscFunctionReturn(0); 3070 } 3071 3072 /* 3073 This variable is used to prevent counting of MatAssemblyBegin() that 3074 are called from within a MatAssemblyEnd(). 3075 */ 3076 static int MatAssemblyEnd_InUse = 0; 3077 #undef __FUNCT__ 3078 #define __FUNCT__ "MatAssemblyBegin" 3079 /*@ 3080 MatAssemblyBegin - Begins assembling the matrix. This routine should 3081 be called after completing all calls to MatSetValues(). 3082 3083 Collective on Mat 3084 3085 Input Parameters: 3086 + mat - the matrix 3087 - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 3088 3089 Notes: 3090 MatSetValues() generally caches the values. The matrix is ready to 3091 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 3092 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 3093 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 3094 using the matrix. 3095 3096 Level: beginner 3097 3098 Concepts: matrices^assembling 3099 3100 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled() 3101 @*/ 3102 int MatAssemblyBegin(Mat mat,MatAssemblyType type) 3103 { 3104 int ierr; 3105 3106 PetscFunctionBegin; 3107 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3108 PetscValidType(mat); 3109 MatPreallocated(mat); 3110 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?"); 3111 if (mat->assembled) { 3112 mat->was_assembled = PETSC_TRUE; 3113 mat->assembled = PETSC_FALSE; 3114 } 3115 if (!MatAssemblyEnd_InUse) { 3116 ierr = PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);CHKERRQ(ierr); 3117 if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);} 3118 ierr = PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);CHKERRQ(ierr); 3119 } else { 3120 if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);} 3121 } 3122 PetscFunctionReturn(0); 3123 } 3124 3125 #undef __FUNCT__ 3126 #define __FUNCT__ "MatAssembed" 3127 /*@ 3128 MatAssembled - Indicates if a matrix has been assembled and is ready for 3129 use; for example, in matrix-vector product. 3130 3131 Collective on Mat 3132 3133 Input Parameter: 3134 . mat - the matrix 3135 3136 Output Parameter: 3137 . assembled - PETSC_TRUE or PETSC_FALSE 3138 3139 Level: advanced 3140 3141 Concepts: matrices^assembled? 3142 3143 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin() 3144 @*/ 3145 int MatAssembled(Mat mat,PetscTruth *assembled) 3146 { 3147 PetscFunctionBegin; 3148 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3149 PetscValidType(mat); 3150 MatPreallocated(mat); 3151 *assembled = mat->assembled; 3152 PetscFunctionReturn(0); 3153 } 3154 3155 #undef __FUNCT__ 3156 #define __FUNCT__ "MatView_Private" 3157 /* 3158 Processes command line options to determine if/how a matrix 3159 is to be viewed. Called by MatAssemblyEnd() and MatLoad(). 3160 */ 3161 int MatView_Private(Mat mat) 3162 { 3163 int ierr; 3164 PetscTruth flg; 3165 static PetscTruth incall = PETSC_FALSE; 3166 3167 PetscFunctionBegin; 3168 if (incall) PetscFunctionReturn(0); 3169 incall = PETSC_TRUE; 3170 ierr = PetscOptionsBegin(mat->comm,mat->prefix,"Matrix Options","Mat");CHKERRQ(ierr); 3171 ierr = PetscOptionsName("-mat_view_info","Information on matrix size","MatView",&flg);CHKERRQ(ierr); 3172 if (flg) { 3173 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 3174 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3175 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3176 } 3177 ierr = PetscOptionsName("-mat_view_info_detailed","Nonzeros in the matrix","MatView",&flg);CHKERRQ(ierr); 3178 if (flg) { 3179 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO_DETAIL);CHKERRQ(ierr); 3180 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3181 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3182 } 3183 ierr = PetscOptionsName("-mat_view","Print matrix to stdout","MatView",&flg);CHKERRQ(ierr); 3184 if (flg) { 3185 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3186 } 3187 ierr = PetscOptionsName("-mat_view_matlab","Print matrix to stdout in a format Matlab can read","MatView",&flg);CHKERRQ(ierr); 3188 if (flg) { 3189 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); 3190 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3191 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3192 } 3193 ierr = PetscOptionsName("-mat_view_socket","Send matrix to socket (can be read from matlab)","MatView",&flg);CHKERRQ(ierr); 3194 if (flg) { 3195 ierr = MatView(mat,PETSC_VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr); 3196 ierr = PetscViewerFlush(PETSC_VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr); 3197 } 3198 ierr = PetscOptionsName("-mat_view_binary","Save matrix to file in binary format","MatView",&flg);CHKERRQ(ierr); 3199 if (flg) { 3200 ierr = MatView(mat,PETSC_VIEWER_BINARY_(mat->comm));CHKERRQ(ierr); 3201 ierr = PetscViewerFlush(PETSC_VIEWER_BINARY_(mat->comm));CHKERRQ(ierr); 3202 } 3203 ierr = PetscOptionsEnd();CHKERRQ(ierr); 3204 /* cannot have inside PetscOptionsBegin() because uses PetscOptionsBegin() */ 3205 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_draw",&flg);CHKERRQ(ierr); 3206 if (flg) { 3207 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_contour",&flg);CHKERRQ(ierr); 3208 if (flg) { 3209 PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 3210 } 3211 ierr = MatView(mat,PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3212 ierr = PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3213 if (flg) { 3214 PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3215 } 3216 } 3217 incall = PETSC_FALSE; 3218 PetscFunctionReturn(0); 3219 } 3220 3221 #undef __FUNCT__ 3222 #define __FUNCT__ "MatAssemblyEnd" 3223 /*@ 3224 MatAssemblyEnd - Completes assembling the matrix. This routine should 3225 be called after MatAssemblyBegin(). 3226 3227 Collective on Mat 3228 3229 Input Parameters: 3230 + mat - the matrix 3231 - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 3232 3233 Options Database Keys: 3234 + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly() 3235 . -mat_view_info_detailed - Prints more detailed info 3236 . -mat_view - Prints matrix in ASCII format 3237 . -mat_view_matlab - Prints matrix in Matlab format 3238 . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX(). 3239 . -display <name> - Sets display name (default is host) 3240 . -draw_pause <sec> - Sets number of seconds to pause after display 3241 . -mat_view_socket - Sends matrix to socket, can be accessed from Matlab (see users manual) 3242 . -viewer_socket_machine <machine> 3243 . -viewer_socket_port <port> 3244 . -mat_view_binary - save matrix to file in binary format 3245 - -viewer_binary_filename <name> 3246 3247 Notes: 3248 MatSetValues() generally caches the values. The matrix is ready to 3249 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 3250 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 3251 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 3252 using the matrix. 3253 3254 Level: beginner 3255 3256 .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled(), PetscViewerSocketOpen() 3257 @*/ 3258 int MatAssemblyEnd(Mat mat,MatAssemblyType type) 3259 { 3260 int ierr; 3261 static int inassm = 0; 3262 PetscTruth flg; 3263 3264 PetscFunctionBegin; 3265 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3266 PetscValidType(mat); 3267 MatPreallocated(mat); 3268 3269 inassm++; 3270 MatAssemblyEnd_InUse++; 3271 if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */ 3272 ierr = PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);CHKERRQ(ierr); 3273 if (mat->ops->assemblyend) { 3274 ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr); 3275 } 3276 ierr = PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);CHKERRQ(ierr); 3277 } else { 3278 if (mat->ops->assemblyend) { 3279 ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr); 3280 } 3281 } 3282 3283 /* Flush assembly is not a true assembly */ 3284 if (type != MAT_FLUSH_ASSEMBLY) { 3285 mat->assembled = PETSC_TRUE; mat->num_ass++; 3286 } 3287 mat->insertmode = NOT_SET_VALUES; 3288 MatAssemblyEnd_InUse--; 3289 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3290 if (!mat->symmetric_eternal) { 3291 mat->symmetric_set = PETSC_FALSE; 3292 mat->hermitian_set = PETSC_FALSE; 3293 mat->structurally_symmetric_set = PETSC_FALSE; 3294 } 3295 if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) { 3296 ierr = MatView_Private(mat);CHKERRQ(ierr); 3297 } 3298 inassm--; 3299 ierr = PetscOptionsHasName(mat->prefix,"-help",&flg);CHKERRQ(ierr); 3300 if (flg) { 3301 ierr = MatPrintHelp(mat);CHKERRQ(ierr); 3302 } 3303 PetscFunctionReturn(0); 3304 } 3305 3306 3307 #undef __FUNCT__ 3308 #define __FUNCT__ "MatCompress" 3309 /*@ 3310 MatCompress - Tries to store the matrix in as little space as 3311 possible. May fail if memory is already fully used, since it 3312 tries to allocate new space. 3313 3314 Collective on Mat 3315 3316 Input Parameters: 3317 . mat - the matrix 3318 3319 Level: advanced 3320 3321 @*/ 3322 int MatCompress(Mat mat) 3323 { 3324 int ierr; 3325 3326 PetscFunctionBegin; 3327 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3328 PetscValidType(mat); 3329 MatPreallocated(mat); 3330 if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);} 3331 PetscFunctionReturn(0); 3332 } 3333 3334 #undef __FUNCT__ 3335 #define __FUNCT__ "MatSetOption" 3336 /*@ 3337 MatSetOption - Sets a parameter option for a matrix. Some options 3338 may be specific to certain storage formats. Some options 3339 determine how values will be inserted (or added). Sorted, 3340 row-oriented input will generally assemble the fastest. The default 3341 is row-oriented, nonsorted input. 3342 3343 Collective on Mat 3344 3345 Input Parameters: 3346 + mat - the matrix 3347 - option - the option, one of those listed below (and possibly others), 3348 e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR 3349 3350 Options Describing Matrix Structure: 3351 + MAT_SYMMETRIC - symmetric in terms of both structure and value 3352 . MAT_HERMITIAN - transpose is the complex conjugation 3353 . MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure 3354 . MAT_NOT_SYMMETRIC - not symmetric in value 3355 . MAT_NOT_HERMITIAN - transpose is not the complex conjugation 3356 . MAT_NOT_STRUCTURALLY_SYMMETRIC - not symmetric nonzero structure 3357 . MAT_SYMMETRY_ETERNAL - if you would like the symmetry/Hermitian flag 3358 you set to be kept with all future use of the matrix 3359 including after MatAssemblyBegin/End() which could 3360 potentially change the symmetry structure, i.e. you 3361 KNOW the matrix will ALWAYS have the property you set. 3362 - MAT_NOT_SYMMETRY_ETERNAL - if MatAssemblyBegin/End() is called then the 3363 flags you set will be dropped (in case potentially 3364 the symmetry etc was lost). 3365 3366 Options For Use with MatSetValues(): 3367 Insert a logically dense subblock, which can be 3368 + MAT_ROW_ORIENTED - row-oriented (default) 3369 . MAT_COLUMN_ORIENTED - column-oriented 3370 . MAT_ROWS_SORTED - sorted by row 3371 . MAT_ROWS_UNSORTED - not sorted by row (default) 3372 . MAT_COLUMNS_SORTED - sorted by column 3373 - MAT_COLUMNS_UNSORTED - not sorted by column (default) 3374 3375 Not these options reflect the data you pass in with MatSetValues(); it has 3376 nothing to do with how the data is stored internally in the matrix 3377 data structure. 3378 3379 When (re)assembling a matrix, we can restrict the input for 3380 efficiency/debugging purposes. These options include 3381 + MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be 3382 allowed if they generate a new nonzero 3383 . MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed 3384 . MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if 3385 they generate a nonzero in a new diagonal (for block diagonal format only) 3386 . MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only) 3387 . MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries 3388 . MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry 3389 - MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly 3390 3391 Notes: 3392 Some options are relevant only for particular matrix types and 3393 are thus ignored by others. Other options are not supported by 3394 certain matrix types and will generate an error message if set. 3395 3396 If using a Fortran 77 module to compute a matrix, one may need to 3397 use the column-oriented option (or convert to the row-oriented 3398 format). 3399 3400 MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion 3401 that would generate a new entry in the nonzero structure is instead 3402 ignored. Thus, if memory has not alredy been allocated for this particular 3403 data, then the insertion is ignored. For dense matrices, in which 3404 the entire array is allocated, no entries are ever ignored. 3405 Set after the first MatAssemblyEnd() 3406 3407 MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion 3408 that would generate a new entry in the nonzero structure instead produces 3409 an error. (Currently supported for AIJ and BAIJ formats only.) 3410 This is a useful flag when using SAME_NONZERO_PATTERN in calling 3411 KSPSetOperators() to ensure that the nonzero pattern truely does 3412 remain unchanged. Set after the first MatAssemblyEnd() 3413 3414 MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion 3415 that would generate a new entry that has not been preallocated will 3416 instead produce an error. (Currently supported for AIJ and BAIJ formats 3417 only.) This is a useful flag when debugging matrix memory preallocation. 3418 3419 MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for 3420 other processors should be dropped, rather than stashed. 3421 This is useful if you know that the "owning" processor is also 3422 always generating the correct matrix entries, so that PETSc need 3423 not transfer duplicate entries generated on another processor. 3424 3425 MAT_USE_HASH_TABLE indicates that a hash table be used to improve the 3426 searches during matrix assembly. When this flag is set, the hash table 3427 is created during the first Matrix Assembly. This hash table is 3428 used the next time through, during MatSetVaules()/MatSetVaulesBlocked() 3429 to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag 3430 should be used with MAT_USE_HASH_TABLE flag. This option is currently 3431 supported by MATMPIBAIJ format only. 3432 3433 MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries 3434 are kept in the nonzero structure 3435 3436 MAT_IGNORE_ZERO_ENTRIES - for AIJ matrices this will stop zero values from creating 3437 a zero location in the matrix 3438 3439 MAT_USE_INODES - indicates using inode version of the code - works with AIJ and 3440 ROWBS matrix types 3441 3442 MAT_DO_NOT_USE_INODES - indicates not using inode version of the code - works 3443 with AIJ and ROWBS matrix types 3444 3445 Level: intermediate 3446 3447 Concepts: matrices^setting options 3448 3449 @*/ 3450 int MatSetOption(Mat mat,MatOption op) 3451 { 3452 int ierr; 3453 3454 PetscFunctionBegin; 3455 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3456 PetscValidType(mat); 3457 MatPreallocated(mat); 3458 switch (op) { 3459 case MAT_SYMMETRIC: 3460 mat->symmetric = PETSC_TRUE; 3461 mat->structurally_symmetric = PETSC_TRUE; 3462 mat->symmetric_set = PETSC_TRUE; 3463 mat->structurally_symmetric_set = PETSC_TRUE; 3464 break; 3465 case MAT_HERMITIAN: 3466 mat->hermitian = PETSC_TRUE; 3467 mat->structurally_symmetric = PETSC_TRUE; 3468 mat->hermitian_set = PETSC_TRUE; 3469 mat->structurally_symmetric_set = PETSC_TRUE; 3470 break; 3471 case MAT_STRUCTURALLY_SYMMETRIC: 3472 mat->structurally_symmetric = PETSC_TRUE; 3473 mat->structurally_symmetric_set = PETSC_TRUE; 3474 break; 3475 case MAT_NOT_SYMMETRIC: 3476 mat->symmetric = PETSC_FALSE; 3477 mat->symmetric_set = PETSC_TRUE; 3478 break; 3479 case MAT_NOT_HERMITIAN: 3480 mat->hermitian = PETSC_FALSE; 3481 mat->hermitian_set = PETSC_TRUE; 3482 break; 3483 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 3484 mat->structurally_symmetric = PETSC_FALSE; 3485 mat->structurally_symmetric_set = PETSC_TRUE; 3486 break; 3487 case MAT_SYMMETRY_ETERNAL: 3488 mat->symmetric_eternal = PETSC_TRUE; 3489 case MAT_NOT_SYMMETRY_ETERNAL: 3490 mat->symmetric_eternal = PETSC_FALSE; 3491 default: 3492 break; 3493 } 3494 if (mat->ops->setoption) { 3495 ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr); 3496 } 3497 PetscFunctionReturn(0); 3498 } 3499 3500 #undef __FUNCT__ 3501 #define __FUNCT__ "MatZeroEntries" 3502 /*@ 3503 MatZeroEntries - Zeros all entries of a matrix. For sparse matrices 3504 this routine retains the old nonzero structure. 3505 3506 Collective on Mat 3507 3508 Input Parameters: 3509 . mat - the matrix 3510 3511 Level: intermediate 3512 3513 Concepts: matrices^zeroing 3514 3515 .seealso: MatZeroRows() 3516 @*/ 3517 int MatZeroEntries(Mat mat) 3518 { 3519 int ierr; 3520 3521 PetscFunctionBegin; 3522 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3523 PetscValidType(mat); 3524 MatPreallocated(mat); 3525 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3526 if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3527 3528 ierr = PetscLogEventBegin(MAT_ZeroEntries,mat,0,0,0);CHKERRQ(ierr); 3529 ierr = (*mat->ops->zeroentries)(mat);CHKERRQ(ierr); 3530 ierr = PetscLogEventEnd(MAT_ZeroEntries,mat,0,0,0);CHKERRQ(ierr); 3531 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3532 PetscFunctionReturn(0); 3533 } 3534 3535 #undef __FUNCT__ 3536 #define __FUNCT__ "MatZeroRows" 3537 /*@C 3538 MatZeroRows - Zeros all entries (except possibly the main diagonal) 3539 of a set of rows of a matrix. 3540 3541 Collective on Mat 3542 3543 Input Parameters: 3544 + mat - the matrix 3545 . is - index set of rows to remove 3546 - diag - pointer to value put in all diagonals of eliminated rows. 3547 Note that diag is not a pointer to an array, but merely a 3548 pointer to a single value. 3549 3550 Notes: 3551 For the AIJ and BAIJ matrix formats this removes the old nonzero structure, 3552 but does not release memory. For the dense and block diagonal 3553 formats this does not alter the nonzero structure. 3554 3555 If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure 3556 of the matrix is not changed (even for AIJ and BAIJ matrices) the values are 3557 merely zeroed. 3558 3559 The user can set a value in the diagonal entry (or for the AIJ and 3560 row formats can optionally remove the main diagonal entry from the 3561 nonzero structure as well, by passing a null pointer (PETSC_NULL 3562 in C or PETSC_NULL_SCALAR in Fortran) as the final argument). 3563 3564 For the parallel case, all processes that share the matrix (i.e., 3565 those in the communicator used for matrix creation) MUST call this 3566 routine, regardless of whether any rows being zeroed are owned by 3567 them. 3568 3569 Level: intermediate 3570 3571 Concepts: matrices^zeroing rows 3572 3573 .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption() 3574 @*/ 3575 int MatZeroRows(Mat mat,IS is,const PetscScalar *diag) 3576 { 3577 int ierr; 3578 3579 PetscFunctionBegin; 3580 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3581 PetscValidType(mat); 3582 MatPreallocated(mat); 3583 PetscValidHeaderSpecific(is,IS_COOKIE); 3584 if (diag) PetscValidScalarPointer(diag); 3585 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3586 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3587 if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3588 3589 ierr = (*mat->ops->zerorows)(mat,is,diag);CHKERRQ(ierr); 3590 ierr = MatView_Private(mat);CHKERRQ(ierr); 3591 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3592 PetscFunctionReturn(0); 3593 } 3594 3595 #undef __FUNCT__ 3596 #define __FUNCT__ "MatZeroRowsLocal" 3597 /*@C 3598 MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal) 3599 of a set of rows of a matrix; using local numbering of rows. 3600 3601 Collective on Mat 3602 3603 Input Parameters: 3604 + mat - the matrix 3605 . is - index set of rows to remove 3606 - diag - pointer to value put in all diagonals of eliminated rows. 3607 Note that diag is not a pointer to an array, but merely a 3608 pointer to a single value. 3609 3610 Notes: 3611 Before calling MatZeroRowsLocal(), the user must first set the 3612 local-to-global mapping by calling MatSetLocalToGlobalMapping(). 3613 3614 For the AIJ matrix formats this removes the old nonzero structure, 3615 but does not release memory. For the dense and block diagonal 3616 formats this does not alter the nonzero structure. 3617 3618 If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure 3619 of the matrix is not changed (even for AIJ and BAIJ matrices) the values are 3620 merely zeroed. 3621 3622 The user can set a value in the diagonal entry (or for the AIJ and 3623 row formats can optionally remove the main diagonal entry from the 3624 nonzero structure as well, by passing a null pointer (PETSC_NULL 3625 in C or PETSC_NULL_SCALAR in Fortran) as the final argument). 3626 3627 Level: intermediate 3628 3629 Concepts: matrices^zeroing 3630 3631 .seealso: MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping 3632 @*/ 3633 int MatZeroRowsLocal(Mat mat,IS is,const PetscScalar *diag) 3634 { 3635 int ierr; 3636 IS newis; 3637 3638 PetscFunctionBegin; 3639 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3640 PetscValidType(mat); 3641 MatPreallocated(mat); 3642 PetscValidHeaderSpecific(is,IS_COOKIE); 3643 if (diag) PetscValidScalarPointer(diag); 3644 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3645 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3646 3647 if (mat->ops->zerorowslocal) { 3648 ierr = (*mat->ops->zerorowslocal)(mat,is,diag);CHKERRQ(ierr); 3649 } else { 3650 if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first"); 3651 ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);CHKERRQ(ierr); 3652 ierr = (*mat->ops->zerorows)(mat,newis,diag);CHKERRQ(ierr); 3653 ierr = ISDestroy(newis);CHKERRQ(ierr); 3654 } 3655 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3656 PetscFunctionReturn(0); 3657 } 3658 3659 #undef __FUNCT__ 3660 #define __FUNCT__ "MatGetSize" 3661 /*@ 3662 MatGetSize - Returns the numbers of rows and columns in a matrix. 3663 3664 Not Collective 3665 3666 Input Parameter: 3667 . mat - the matrix 3668 3669 Output Parameters: 3670 + m - the number of global rows 3671 - n - the number of global columns 3672 3673 Level: beginner 3674 3675 Concepts: matrices^size 3676 3677 .seealso: MatGetLocalSize() 3678 @*/ 3679 int MatGetSize(Mat mat,int *m,int* n) 3680 { 3681 PetscFunctionBegin; 3682 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3683 if (m) *m = mat->M; 3684 if (n) *n = mat->N; 3685 PetscFunctionReturn(0); 3686 } 3687 3688 #undef __FUNCT__ 3689 #define __FUNCT__ "MatGetLocalSize" 3690 /*@ 3691 MatGetLocalSize - Returns the number of rows and columns in a matrix 3692 stored locally. This information may be implementation dependent, so 3693 use with care. 3694 3695 Not Collective 3696 3697 Input Parameters: 3698 . mat - the matrix 3699 3700 Output Parameters: 3701 + m - the number of local rows 3702 - n - the number of local columns 3703 3704 Level: beginner 3705 3706 Concepts: matrices^local size 3707 3708 .seealso: MatGetSize() 3709 @*/ 3710 int MatGetLocalSize(Mat mat,int *m,int* n) 3711 { 3712 PetscFunctionBegin; 3713 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3714 if (m) *m = mat->m; 3715 if (n) *n = mat->n; 3716 PetscFunctionReturn(0); 3717 } 3718 3719 #undef __FUNCT__ 3720 #define __FUNCT__ "MatGetOwnershipRange" 3721 /*@ 3722 MatGetOwnershipRange - Returns the range of matrix rows owned by 3723 this processor, assuming that the matrix is laid out with the first 3724 n1 rows on the first processor, the next n2 rows on the second, etc. 3725 For certain parallel layouts this range may not be well defined. 3726 3727 Not Collective 3728 3729 Input Parameters: 3730 . mat - the matrix 3731 3732 Output Parameters: 3733 + m - the global index of the first local row 3734 - n - one more than the global index of the last local row 3735 3736 Level: beginner 3737 3738 Concepts: matrices^row ownership 3739 @*/ 3740 int MatGetOwnershipRange(Mat mat,int *m,int* n) 3741 { 3742 int ierr; 3743 3744 PetscFunctionBegin; 3745 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3746 PetscValidType(mat); 3747 MatPreallocated(mat); 3748 if (m) PetscValidIntPointer(m); 3749 if (n) PetscValidIntPointer(n); 3750 ierr = PetscMapGetLocalRange(mat->rmap,m,n);CHKERRQ(ierr); 3751 PetscFunctionReturn(0); 3752 } 3753 3754 #undef __FUNCT__ 3755 #define __FUNCT__ "MatILUFactorSymbolic" 3756 /*@ 3757 MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix. 3758 Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric() 3759 to complete the factorization. 3760 3761 Collective on Mat 3762 3763 Input Parameters: 3764 + mat - the matrix 3765 . row - row permutation 3766 . column - column permutation 3767 - info - structure containing 3768 $ levels - number of levels of fill. 3769 $ expected fill - as ratio of original fill. 3770 $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices 3771 missing diagonal entries) 3772 3773 Output Parameters: 3774 . fact - new matrix that has been symbolically factored 3775 3776 Notes: 3777 See the users manual for additional information about 3778 choosing the fill factor for better efficiency. 3779 3780 Most users should employ the simplified KSP interface for linear solvers 3781 instead of working directly with matrix algebra routines such as this. 3782 See, e.g., KSPCreate(). 3783 3784 Level: developer 3785 3786 Concepts: matrices^symbolic LU factorization 3787 Concepts: matrices^factorization 3788 Concepts: LU^symbolic factorization 3789 3790 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 3791 MatGetOrdering(), MatFactorInfo 3792 3793 @*/ 3794 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatFactorInfo *info,Mat *fact) 3795 { 3796 int ierr; 3797 3798 PetscFunctionBegin; 3799 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3800 PetscValidType(mat); 3801 MatPreallocated(mat); 3802 PetscValidPointer(fact); 3803 PetscValidHeaderSpecific(row,IS_COOKIE); 3804 PetscValidHeaderSpecific(col,IS_COOKIE); 3805 if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %d",(int)info->levels); 3806 if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill); 3807 if (!mat->ops->ilufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ILU",mat->type_name); 3808 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3809 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3810 3811 ierr = PetscLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 3812 ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr); 3813 ierr = PetscLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 3814 PetscFunctionReturn(0); 3815 } 3816 3817 #undef __FUNCT__ 3818 #define __FUNCT__ "MatICCFactorSymbolic" 3819 /*@ 3820 MatICCFactorSymbolic - Performs symbolic incomplete 3821 Cholesky factorization for a symmetric matrix. Use 3822 MatCholeskyFactorNumeric() to complete the factorization. 3823 3824 Collective on Mat 3825 3826 Input Parameters: 3827 + mat - the matrix 3828 . perm - row and column permutation 3829 - info - structure containing 3830 $ levels - number of levels of fill. 3831 $ expected fill - as ratio of original fill. 3832 3833 Output Parameter: 3834 . fact - the factored matrix 3835 3836 Notes: 3837 Currently only no-fill factorization is supported. 3838 3839 Most users should employ the simplified KSP interface for linear solvers 3840 instead of working directly with matrix algebra routines such as this. 3841 See, e.g., KSPCreate(). 3842 3843 Level: developer 3844 3845 Concepts: matrices^symbolic incomplete Cholesky factorization 3846 Concepts: matrices^factorization 3847 Concepts: Cholsky^symbolic factorization 3848 3849 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor(), MatFactorInfo 3850 @*/ 3851 int MatICCFactorSymbolic(Mat mat,IS perm,MatFactorInfo *info,Mat *fact) 3852 { 3853 int ierr; 3854 3855 PetscFunctionBegin; 3856 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3857 PetscValidType(mat); 3858 MatPreallocated(mat); 3859 PetscValidPointer(fact); 3860 PetscValidHeaderSpecific(perm,IS_COOKIE); 3861 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3862 if (info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels negative %d",(int) info->levels); 3863 if (info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill); 3864 if (!mat->ops->iccfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ICC",mat->type_name); 3865 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3866 3867 ierr = PetscLogEventBegin(MAT_ICCFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 3868 ierr = (*mat->ops->iccfactorsymbolic)(mat,perm,info,fact);CHKERRQ(ierr); 3869 ierr = PetscLogEventEnd(MAT_ICCFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 3870 PetscFunctionReturn(0); 3871 } 3872 3873 #undef __FUNCT__ 3874 #define __FUNCT__ "MatGetArray" 3875 /*@C 3876 MatGetArray - Returns a pointer to the element values in the matrix. 3877 The result of this routine is dependent on the underlying matrix data 3878 structure, and may not even work for certain matrix types. You MUST 3879 call MatRestoreArray() when you no longer need to access the array. 3880 3881 Not Collective 3882 3883 Input Parameter: 3884 . mat - the matrix 3885 3886 Output Parameter: 3887 . v - the location of the values 3888 3889 3890 Fortran Note: 3891 This routine is used differently from Fortran, e.g., 3892 .vb 3893 Mat mat 3894 PetscScalar mat_array(1) 3895 PetscOffset i_mat 3896 int ierr 3897 call MatGetArray(mat,mat_array,i_mat,ierr) 3898 3899 C Access first local entry in matrix; note that array is 3900 C treated as one dimensional 3901 value = mat_array(i_mat + 1) 3902 3903 [... other code ...] 3904 call MatRestoreArray(mat,mat_array,i_mat,ierr) 3905 .ve 3906 3907 See the Fortran chapter of the users manual and 3908 petsc/src/mat/examples/tests for details. 3909 3910 Level: advanced 3911 3912 Concepts: matrices^access array 3913 3914 .seealso: MatRestoreArray(), MatGetArrayF90() 3915 @*/ 3916 int MatGetArray(Mat mat,PetscScalar *v[]) 3917 { 3918 int ierr; 3919 3920 PetscFunctionBegin; 3921 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3922 PetscValidType(mat); 3923 MatPreallocated(mat); 3924 PetscValidPointer(v); 3925 if (!mat->ops->getarray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3926 ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr); 3927 PetscFunctionReturn(0); 3928 } 3929 3930 #undef __FUNCT__ 3931 #define __FUNCT__ "MatRestoreArray" 3932 /*@C 3933 MatRestoreArray - Restores the matrix after MatGetArray() has been called. 3934 3935 Not Collective 3936 3937 Input Parameter: 3938 + mat - the matrix 3939 - v - the location of the values 3940 3941 Fortran Note: 3942 This routine is used differently from Fortran, e.g., 3943 .vb 3944 Mat mat 3945 PetscScalar mat_array(1) 3946 PetscOffset i_mat 3947 int ierr 3948 call MatGetArray(mat,mat_array,i_mat,ierr) 3949 3950 C Access first local entry in matrix; note that array is 3951 C treated as one dimensional 3952 value = mat_array(i_mat + 1) 3953 3954 [... other code ...] 3955 call MatRestoreArray(mat,mat_array,i_mat,ierr) 3956 .ve 3957 3958 See the Fortran chapter of the users manual and 3959 petsc/src/mat/examples/tests for details 3960 3961 Level: advanced 3962 3963 .seealso: MatGetArray(), MatRestoreArrayF90() 3964 @*/ 3965 int MatRestoreArray(Mat mat,PetscScalar *v[]) 3966 { 3967 int ierr; 3968 3969 PetscFunctionBegin; 3970 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3971 PetscValidType(mat); 3972 MatPreallocated(mat); 3973 PetscValidPointer(v); 3974 #if defined(PETSC_USE_BOPT_g) 3975 CHKMEMQ; 3976 #endif 3977 if (!mat->ops->restorearray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3978 ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr); 3979 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 3980 PetscFunctionReturn(0); 3981 } 3982 3983 #undef __FUNCT__ 3984 #define __FUNCT__ "MatGetSubMatrices" 3985 /*@C 3986 MatGetSubMatrices - Extracts several submatrices from a matrix. If submat 3987 points to an array of valid matrices, they may be reused to store the new 3988 submatrices. 3989 3990 Collective on Mat 3991 3992 Input Parameters: 3993 + mat - the matrix 3994 . n - the number of submatrixes to be extracted (on this processor, may be zero) 3995 . irow, icol - index sets of rows and columns to extract 3996 - scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 3997 3998 Output Parameter: 3999 . submat - the array of submatrices 4000 4001 Notes: 4002 MatGetSubMatrices() can extract only sequential submatrices 4003 (from both sequential and parallel matrices). Use MatGetSubMatrix() 4004 to extract a parallel submatrix. 4005 4006 When extracting submatrices from a parallel matrix, each processor can 4007 form a different submatrix by setting the rows and columns of its 4008 individual index sets according to the local submatrix desired. 4009 4010 When finished using the submatrices, the user should destroy 4011 them with MatDestroyMatrices(). 4012 4013 MAT_REUSE_MATRIX can only be used when the nonzero structure of the 4014 original matrix has not changed from that last call to MatGetSubMatrices(). 4015 4016 This routine creates the matrices in submat; you should NOT create them before 4017 calling it. It also allocates the array of matrix pointers submat. 4018 4019 Fortran Note: 4020 The Fortran interface is slightly different from that given below; it 4021 requires one to pass in as submat a Mat (integer) array of size at least m. 4022 4023 Level: advanced 4024 4025 Concepts: matrices^accessing submatrices 4026 Concepts: submatrices 4027 4028 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal() 4029 @*/ 4030 int MatGetSubMatrices(Mat mat,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *submat[]) 4031 { 4032 int ierr; 4033 4034 PetscFunctionBegin; 4035 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4036 PetscValidType(mat); 4037 MatPreallocated(mat); 4038 if (!mat->ops->getsubmatrices) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4039 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4040 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4041 4042 ierr = PetscLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);CHKERRQ(ierr); 4043 ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);CHKERRQ(ierr); 4044 ierr = PetscLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);CHKERRQ(ierr); 4045 PetscFunctionReturn(0); 4046 } 4047 4048 #undef __FUNCT__ 4049 #define __FUNCT__ "MatDestroyMatrices" 4050 /*@C 4051 MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices(). 4052 4053 Collective on Mat 4054 4055 Input Parameters: 4056 + n - the number of local matrices 4057 - mat - the matrices (note that this is a pointer to the array of matrices, just to match the calling 4058 sequence of MatGetSubMatrices()) 4059 4060 Level: advanced 4061 4062 Notes: Frees not only the matrices, but also the array that contains the matrices 4063 4064 .seealso: MatGetSubMatrices() 4065 @*/ 4066 int MatDestroyMatrices(int n,Mat *mat[]) 4067 { 4068 int ierr,i; 4069 4070 PetscFunctionBegin; 4071 if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of matrices %d",n); 4072 PetscValidPointer(mat); 4073 for (i=0; i<n; i++) { 4074 ierr = MatDestroy((*mat)[i]);CHKERRQ(ierr); 4075 } 4076 /* memory is allocated even if n = 0 */ 4077 ierr = PetscFree(*mat);CHKERRQ(ierr); 4078 PetscFunctionReturn(0); 4079 } 4080 4081 #undef __FUNCT__ 4082 #define __FUNCT__ "MatIncreaseOverlap" 4083 /*@ 4084 MatIncreaseOverlap - Given a set of submatrices indicated by index sets, 4085 replaces the index sets by larger ones that represent submatrices with 4086 additional overlap. 4087 4088 Collective on Mat 4089 4090 Input Parameters: 4091 + mat - the matrix 4092 . n - the number of index sets 4093 . is - the array of index sets (these index sets will changed during the call) 4094 - ov - the additional overlap requested 4095 4096 Level: developer 4097 4098 Concepts: overlap 4099 Concepts: ASM^computing overlap 4100 4101 .seealso: MatGetSubMatrices() 4102 @*/ 4103 int MatIncreaseOverlap(Mat mat,int n,IS is[],int ov) 4104 { 4105 int ierr; 4106 4107 PetscFunctionBegin; 4108 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4109 PetscValidType(mat); 4110 MatPreallocated(mat); 4111 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4112 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4113 4114 if (!ov) PetscFunctionReturn(0); 4115 if (!mat->ops->increaseoverlap) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4116 ierr = PetscLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);CHKERRQ(ierr); 4117 ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov);CHKERRQ(ierr); 4118 ierr = PetscLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);CHKERRQ(ierr); 4119 PetscFunctionReturn(0); 4120 } 4121 4122 #undef __FUNCT__ 4123 #define __FUNCT__ "MatPrintHelp" 4124 /*@ 4125 MatPrintHelp - Prints all the options for the matrix. 4126 4127 Collective on Mat 4128 4129 Input Parameter: 4130 . mat - the matrix 4131 4132 Options Database Keys: 4133 + -help - Prints matrix options 4134 - -h - Prints matrix options 4135 4136 Level: developer 4137 4138 .seealso: MatCreate(), MatCreateXXX() 4139 @*/ 4140 int MatPrintHelp(Mat mat) 4141 { 4142 static PetscTruth called = PETSC_FALSE; 4143 int ierr; 4144 4145 PetscFunctionBegin; 4146 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4147 PetscValidType(mat); 4148 MatPreallocated(mat); 4149 4150 if (!called) { 4151 if (mat->ops->printhelp) { 4152 ierr = (*mat->ops->printhelp)(mat);CHKERRQ(ierr); 4153 } 4154 called = PETSC_TRUE; 4155 } 4156 PetscFunctionReturn(0); 4157 } 4158 4159 #undef __FUNCT__ 4160 #define __FUNCT__ "MatGetBlockSize" 4161 /*@ 4162 MatGetBlockSize - Returns the matrix block size; useful especially for the 4163 block row and block diagonal formats. 4164 4165 Not Collective 4166 4167 Input Parameter: 4168 . mat - the matrix 4169 4170 Output Parameter: 4171 . bs - block size 4172 4173 Notes: 4174 Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG. 4175 Block row formats are MATSEQBAIJ, MATMPIBAIJ 4176 4177 Level: intermediate 4178 4179 Concepts: matrices^block size 4180 4181 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag() 4182 @*/ 4183 int MatGetBlockSize(Mat mat,int *bs) 4184 { 4185 int ierr; 4186 4187 PetscFunctionBegin; 4188 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4189 PetscValidType(mat); 4190 MatPreallocated(mat); 4191 PetscValidIntPointer(bs); 4192 if (!mat->ops->getblocksize) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4193 ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr); 4194 PetscFunctionReturn(0); 4195 } 4196 4197 #undef __FUNCT__ 4198 #define __FUNCT__ "MatGetRowIJ" 4199 /*@C 4200 MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices. 4201 4202 Collective on Mat 4203 4204 Input Parameters: 4205 + mat - the matrix 4206 . shift - 0 or 1 indicating we want the indices starting at 0 or 1 4207 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4208 symmetrized 4209 4210 Output Parameters: 4211 + n - number of rows in the (possibly compressed) matrix 4212 . ia - the row pointers 4213 . ja - the column indices 4214 - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned 4215 4216 Level: developer 4217 4218 .seealso: MatGetColumnIJ(), MatRestoreRowIJ() 4219 @*/ 4220 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int *ia[],int* ja[],PetscTruth *done) 4221 { 4222 int ierr; 4223 4224 PetscFunctionBegin; 4225 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4226 PetscValidType(mat); 4227 MatPreallocated(mat); 4228 if (ia) PetscValidIntPointer(ia); 4229 if (ja) PetscValidIntPointer(ja); 4230 PetscValidIntPointer(done); 4231 if (!mat->ops->getrowij) *done = PETSC_FALSE; 4232 else { 4233 *done = PETSC_TRUE; 4234 ierr = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4235 } 4236 PetscFunctionReturn(0); 4237 } 4238 4239 #undef __FUNCT__ 4240 #define __FUNCT__ "MatGetColumnIJ" 4241 /*@C 4242 MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices. 4243 4244 Collective on Mat 4245 4246 Input Parameters: 4247 + mat - the matrix 4248 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4249 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4250 symmetrized 4251 4252 Output Parameters: 4253 + n - number of columns in the (possibly compressed) matrix 4254 . ia - the column pointers 4255 . ja - the row indices 4256 - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned 4257 4258 Level: developer 4259 4260 .seealso: MatGetRowIJ(), MatRestoreColumnIJ() 4261 @*/ 4262 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int *ia[],int* ja[],PetscTruth *done) 4263 { 4264 int ierr; 4265 4266 PetscFunctionBegin; 4267 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4268 PetscValidType(mat); 4269 MatPreallocated(mat); 4270 if (ia) PetscValidIntPointer(ia); 4271 if (ja) PetscValidIntPointer(ja); 4272 PetscValidIntPointer(done); 4273 4274 if (!mat->ops->getcolumnij) *done = PETSC_FALSE; 4275 else { 4276 *done = PETSC_TRUE; 4277 ierr = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4278 } 4279 PetscFunctionReturn(0); 4280 } 4281 4282 #undef __FUNCT__ 4283 #define __FUNCT__ "MatRestoreRowIJ" 4284 /*@C 4285 MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with 4286 MatGetRowIJ(). 4287 4288 Collective on Mat 4289 4290 Input Parameters: 4291 + mat - the matrix 4292 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4293 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4294 symmetrized 4295 4296 Output Parameters: 4297 + n - size of (possibly compressed) matrix 4298 . ia - the row pointers 4299 . ja - the column indices 4300 - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 4301 4302 Level: developer 4303 4304 .seealso: MatGetRowIJ(), MatRestoreColumnIJ() 4305 @*/ 4306 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int *ia[],int* ja[],PetscTruth *done) 4307 { 4308 int ierr; 4309 4310 PetscFunctionBegin; 4311 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4312 PetscValidType(mat); 4313 MatPreallocated(mat); 4314 if (ia) PetscValidIntPointer(ia); 4315 if (ja) PetscValidIntPointer(ja); 4316 PetscValidIntPointer(done); 4317 4318 if (!mat->ops->restorerowij) *done = PETSC_FALSE; 4319 else { 4320 *done = PETSC_TRUE; 4321 ierr = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4322 } 4323 PetscFunctionReturn(0); 4324 } 4325 4326 #undef __FUNCT__ 4327 #define __FUNCT__ "MatRestoreColumnIJ" 4328 /*@C 4329 MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with 4330 MatGetColumnIJ(). 4331 4332 Collective on Mat 4333 4334 Input Parameters: 4335 + mat - the matrix 4336 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4337 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4338 symmetrized 4339 4340 Output Parameters: 4341 + n - size of (possibly compressed) matrix 4342 . ia - the column pointers 4343 . ja - the row indices 4344 - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 4345 4346 Level: developer 4347 4348 .seealso: MatGetColumnIJ(), MatRestoreRowIJ() 4349 @*/ 4350 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int *ia[],int* ja[],PetscTruth *done) 4351 { 4352 int ierr; 4353 4354 PetscFunctionBegin; 4355 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4356 PetscValidType(mat); 4357 MatPreallocated(mat); 4358 if (ia) PetscValidIntPointer(ia); 4359 if (ja) PetscValidIntPointer(ja); 4360 PetscValidIntPointer(done); 4361 4362 if (!mat->ops->restorecolumnij) *done = PETSC_FALSE; 4363 else { 4364 *done = PETSC_TRUE; 4365 ierr = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4366 } 4367 PetscFunctionReturn(0); 4368 } 4369 4370 #undef __FUNCT__ 4371 #define __FUNCT__ "MatColoringPatch" 4372 /*@C 4373 MatColoringPatch -Used inside matrix coloring routines that 4374 use MatGetRowIJ() and/or MatGetColumnIJ(). 4375 4376 Collective on Mat 4377 4378 Input Parameters: 4379 + mat - the matrix 4380 . n - number of colors 4381 - colorarray - array indicating color for each column 4382 4383 Output Parameters: 4384 . iscoloring - coloring generated using colorarray information 4385 4386 Level: developer 4387 4388 .seealso: MatGetRowIJ(), MatGetColumnIJ() 4389 4390 @*/ 4391 int MatColoringPatch(Mat mat,int n,int ncolors,const ISColoringValue colorarray[],ISColoring *iscoloring) 4392 { 4393 int ierr; 4394 4395 PetscFunctionBegin; 4396 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4397 PetscValidType(mat); 4398 MatPreallocated(mat); 4399 PetscValidIntPointer(colorarray); 4400 4401 if (!mat->ops->coloringpatch){ 4402 ierr = ISColoringCreate(mat->comm,n,colorarray,iscoloring);CHKERRQ(ierr); 4403 } else { 4404 ierr = (*mat->ops->coloringpatch)(mat,n,ncolors,colorarray,iscoloring);CHKERRQ(ierr); 4405 } 4406 PetscFunctionReturn(0); 4407 } 4408 4409 4410 #undef __FUNCT__ 4411 #define __FUNCT__ "MatSetUnfactored" 4412 /*@ 4413 MatSetUnfactored - Resets a factored matrix to be treated as unfactored. 4414 4415 Collective on Mat 4416 4417 Input Parameter: 4418 . mat - the factored matrix to be reset 4419 4420 Notes: 4421 This routine should be used only with factored matrices formed by in-place 4422 factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE 4423 format). This option can save memory, for example, when solving nonlinear 4424 systems with a matrix-free Newton-Krylov method and a matrix-based, in-place 4425 ILU(0) preconditioner. 4426 4427 Note that one can specify in-place ILU(0) factorization by calling 4428 .vb 4429 PCType(pc,PCILU); 4430 PCILUSeUseInPlace(pc); 4431 .ve 4432 or by using the options -pc_type ilu -pc_ilu_in_place 4433 4434 In-place factorization ILU(0) can also be used as a local 4435 solver for the blocks within the block Jacobi or additive Schwarz 4436 methods (runtime option: -sub_pc_ilu_in_place). See the discussion 4437 of these preconditioners in the users manual for details on setting 4438 local solver options. 4439 4440 Most users should employ the simplified KSP interface for linear solvers 4441 instead of working directly with matrix algebra routines such as this. 4442 See, e.g., KSPCreate(). 4443 4444 Level: developer 4445 4446 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace() 4447 4448 Concepts: matrices^unfactored 4449 4450 @*/ 4451 int MatSetUnfactored(Mat mat) 4452 { 4453 int ierr; 4454 4455 PetscFunctionBegin; 4456 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4457 PetscValidType(mat); 4458 MatPreallocated(mat); 4459 mat->factor = 0; 4460 if (!mat->ops->setunfactored) PetscFunctionReturn(0); 4461 ierr = (*mat->ops->setunfactored)(mat);CHKERRQ(ierr); 4462 PetscFunctionReturn(0); 4463 } 4464 4465 /*MC 4466 MatGetArrayF90 - Accesses a matrix array from Fortran90. 4467 4468 Synopsis: 4469 MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 4470 4471 Not collective 4472 4473 Input Parameter: 4474 . x - matrix 4475 4476 Output Parameters: 4477 + xx_v - the Fortran90 pointer to the array 4478 - ierr - error code 4479 4480 Example of Usage: 4481 .vb 4482 PetscScalar, pointer xx_v(:) 4483 .... 4484 call MatGetArrayF90(x,xx_v,ierr) 4485 a = xx_v(3) 4486 call MatRestoreArrayF90(x,xx_v,ierr) 4487 .ve 4488 4489 Notes: 4490 Not yet supported for all F90 compilers 4491 4492 Level: advanced 4493 4494 .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray() 4495 4496 Concepts: matrices^accessing array 4497 4498 M*/ 4499 4500 /*MC 4501 MatRestoreArrayF90 - Restores a matrix array that has been 4502 accessed with MatGetArrayF90(). 4503 4504 Synopsis: 4505 MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 4506 4507 Not collective 4508 4509 Input Parameters: 4510 + x - matrix 4511 - xx_v - the Fortran90 pointer to the array 4512 4513 Output Parameter: 4514 . ierr - error code 4515 4516 Example of Usage: 4517 .vb 4518 PetscScalar, pointer xx_v(:) 4519 .... 4520 call MatGetArrayF90(x,xx_v,ierr) 4521 a = xx_v(3) 4522 call MatRestoreArrayF90(x,xx_v,ierr) 4523 .ve 4524 4525 Notes: 4526 Not yet supported for all F90 compilers 4527 4528 Level: advanced 4529 4530 .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray() 4531 4532 M*/ 4533 4534 4535 #undef __FUNCT__ 4536 #define __FUNCT__ "MatGetSubMatrix" 4537 /*@ 4538 MatGetSubMatrix - Gets a single submatrix on the same number of processors 4539 as the original matrix. 4540 4541 Collective on Mat 4542 4543 Input Parameters: 4544 + mat - the original matrix 4545 . isrow - rows this processor should obtain 4546 . iscol - columns for all processors you wish to keep 4547 . csize - number of columns "local" to this processor (does nothing for sequential 4548 matrices). This should match the result from VecGetLocalSize(x,...) if you 4549 plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE 4550 - cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4551 4552 Output Parameter: 4553 . newmat - the new submatrix, of the same type as the old 4554 4555 Level: advanced 4556 4557 Notes: the iscol argument MUST be the same on each processor. You might be 4558 able to create the iscol argument with ISAllGather(). 4559 4560 The first time this is called you should use a cll of MAT_INITIAL_MATRIX, 4561 the MatGetSubMatrix() routine will create the newmat for you. Any additional calls 4562 to this routine with a mat of the same nonzero structure will reuse the matrix 4563 generated the first time. 4564 4565 Concepts: matrices^submatrices 4566 4567 .seealso: MatGetSubMatrices(), ISAllGather() 4568 @*/ 4569 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat) 4570 { 4571 int ierr, size; 4572 Mat *local; 4573 4574 PetscFunctionBegin; 4575 PetscValidType(mat); 4576 MatPreallocated(mat); 4577 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4578 ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr); 4579 4580 /* if original matrix is on just one processor then use submatrix generated */ 4581 if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) { 4582 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr); 4583 PetscFunctionReturn(0); 4584 } else if (!mat->ops->getsubmatrix && size == 1) { 4585 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr); 4586 *newmat = *local; 4587 ierr = PetscFree(local);CHKERRQ(ierr); 4588 PetscFunctionReturn(0); 4589 } 4590 4591 if (!mat->ops->getsubmatrix) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4592 ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr); 4593 ierr = PetscObjectIncreaseState((PetscObject)*newmat); CHKERRQ(ierr); 4594 PetscFunctionReturn(0); 4595 } 4596 4597 #undef __FUNCT__ 4598 #define __FUNCT__ "MatGetPetscMaps" 4599 /*@C 4600 MatGetPetscMaps - Returns the maps associated with the matrix. 4601 4602 Not Collective 4603 4604 Input Parameter: 4605 . mat - the matrix 4606 4607 Output Parameters: 4608 + rmap - the row (right) map 4609 - cmap - the column (left) map 4610 4611 Level: developer 4612 4613 Concepts: maps^getting from matrix 4614 4615 @*/ 4616 int MatGetPetscMaps(Mat mat,PetscMap *rmap,PetscMap *cmap) 4617 { 4618 int ierr; 4619 4620 PetscFunctionBegin; 4621 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4622 PetscValidType(mat); 4623 MatPreallocated(mat); 4624 ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr); 4625 PetscFunctionReturn(0); 4626 } 4627 4628 /* 4629 Version that works for all PETSc matrices 4630 */ 4631 #undef __FUNCT__ 4632 #define __FUNCT__ "MatGetPetscMaps_Petsc" 4633 int MatGetPetscMaps_Petsc(Mat mat,PetscMap *rmap,PetscMap *cmap) 4634 { 4635 PetscFunctionBegin; 4636 if (rmap) *rmap = mat->rmap; 4637 if (cmap) *cmap = mat->cmap; 4638 PetscFunctionReturn(0); 4639 } 4640 4641 #undef __FUNCT__ 4642 #define __FUNCT__ "MatSetStashInitialSize" 4643 /*@ 4644 MatSetStashInitialSize - sets the sizes of the matrix stash, that is 4645 used during the assembly process to store values that belong to 4646 other processors. 4647 4648 Not Collective 4649 4650 Input Parameters: 4651 + mat - the matrix 4652 . size - the initial size of the stash. 4653 - bsize - the initial size of the block-stash(if used). 4654 4655 Options Database Keys: 4656 + -matstash_initial_size <size> or <size0,size1,...sizep-1> 4657 - -matstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1> 4658 4659 Level: intermediate 4660 4661 Notes: 4662 The block-stash is used for values set with VecSetValuesBlocked() while 4663 the stash is used for values set with VecSetValues() 4664 4665 Run with the option -log_info and look for output of the form 4666 MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs. 4667 to determine the appropriate value, MM, to use for size and 4668 MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs. 4669 to determine the value, BMM to use for bsize 4670 4671 Concepts: stash^setting matrix size 4672 Concepts: matrices^stash 4673 4674 @*/ 4675 int MatSetStashInitialSize(Mat mat,int size, int bsize) 4676 { 4677 int ierr; 4678 4679 PetscFunctionBegin; 4680 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4681 PetscValidType(mat); 4682 MatPreallocated(mat); 4683 ierr = MatStashSetInitialSize_Private(&mat->stash,size);CHKERRQ(ierr); 4684 ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize);CHKERRQ(ierr); 4685 PetscFunctionReturn(0); 4686 } 4687 4688 #undef __FUNCT__ 4689 #define __FUNCT__ "MatInterpolateAdd" 4690 /*@ 4691 MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of 4692 the matrix 4693 4694 Collective on Mat 4695 4696 Input Parameters: 4697 + mat - the matrix 4698 . x,y - the vectors 4699 - w - where the result is stored 4700 4701 Level: intermediate 4702 4703 Notes: 4704 w may be the same vector as y. 4705 4706 This allows one to use either the restriction or interpolation (its transpose) 4707 matrix to do the interpolation 4708 4709 Concepts: interpolation 4710 4711 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict() 4712 4713 @*/ 4714 int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w) 4715 { 4716 int M,N,ierr; 4717 4718 PetscFunctionBegin; 4719 PetscValidType(A); 4720 MatPreallocated(A); 4721 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4722 if (N > M) { 4723 ierr = MatMultTransposeAdd(A,x,y,w);CHKERRQ(ierr); 4724 } else { 4725 ierr = MatMultAdd(A,x,y,w);CHKERRQ(ierr); 4726 } 4727 PetscFunctionReturn(0); 4728 } 4729 4730 #undef __FUNCT__ 4731 #define __FUNCT__ "MatInterpolate" 4732 /*@ 4733 MatInterpolate - y = A*x or A'*x depending on the shape of 4734 the matrix 4735 4736 Collective on Mat 4737 4738 Input Parameters: 4739 + mat - the matrix 4740 - x,y - the vectors 4741 4742 Level: intermediate 4743 4744 Notes: 4745 This allows one to use either the restriction or interpolation (its transpose) 4746 matrix to do the interpolation 4747 4748 Concepts: matrices^interpolation 4749 4750 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict() 4751 4752 @*/ 4753 int MatInterpolate(Mat A,Vec x,Vec y) 4754 { 4755 int M,N,ierr; 4756 4757 PetscFunctionBegin; 4758 PetscValidType(A); 4759 MatPreallocated(A); 4760 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4761 if (N > M) { 4762 ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr); 4763 } else { 4764 ierr = MatMult(A,x,y);CHKERRQ(ierr); 4765 } 4766 PetscFunctionReturn(0); 4767 } 4768 4769 #undef __FUNCT__ 4770 #define __FUNCT__ "MatRestrict" 4771 /*@ 4772 MatRestrict - y = A*x or A'*x 4773 4774 Collective on Mat 4775 4776 Input Parameters: 4777 + mat - the matrix 4778 - x,y - the vectors 4779 4780 Level: intermediate 4781 4782 Notes: 4783 This allows one to use either the restriction or interpolation (its transpose) 4784 matrix to do the restriction 4785 4786 Concepts: matrices^restriction 4787 4788 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate() 4789 4790 @*/ 4791 int MatRestrict(Mat A,Vec x,Vec y) 4792 { 4793 int M,N,ierr; 4794 4795 PetscFunctionBegin; 4796 PetscValidType(A); 4797 MatPreallocated(A); 4798 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4799 if (N > M) { 4800 ierr = MatMult(A,x,y);CHKERRQ(ierr); 4801 } else { 4802 ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr); 4803 } 4804 PetscFunctionReturn(0); 4805 } 4806 4807 #undef __FUNCT__ 4808 #define __FUNCT__ "MatNullSpaceAttach" 4809 /*@C 4810 MatNullSpaceAttach - attaches a null space to a matrix. 4811 This null space will be removed from the resulting vector whenever 4812 MatMult() is called 4813 4814 Collective on Mat 4815 4816 Input Parameters: 4817 + mat - the matrix 4818 - nullsp - the null space object 4819 4820 Level: developer 4821 4822 Notes: 4823 Overwrites any previous null space that may have been attached 4824 4825 Concepts: null space^attaching to matrix 4826 4827 .seealso: MatCreate(), MatNullSpaceCreate() 4828 @*/ 4829 int MatNullSpaceAttach(Mat mat,MatNullSpace nullsp) 4830 { 4831 int ierr; 4832 4833 PetscFunctionBegin; 4834 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4835 PetscValidType(mat); 4836 MatPreallocated(mat); 4837 PetscValidHeaderSpecific(nullsp,MAT_NULLSPACE_COOKIE); 4838 4839 if (mat->nullsp) { 4840 ierr = MatNullSpaceDestroy(mat->nullsp);CHKERRQ(ierr); 4841 } 4842 mat->nullsp = nullsp; 4843 ierr = PetscObjectReference((PetscObject)nullsp);CHKERRQ(ierr); 4844 PetscFunctionReturn(0); 4845 } 4846 4847 #undef __FUNCT__ 4848 #define __FUNCT__ "MatICCFactor" 4849 /*@ 4850 MatICCFactor - Performs in-place incomplete Cholesky factorization of matrix. 4851 4852 Collective on Mat 4853 4854 Input Parameters: 4855 + mat - the matrix 4856 . row - row/column permutation 4857 . fill - expected fill factor >= 1.0 4858 - level - level of fill, for ICC(k) 4859 4860 Notes: 4861 Probably really in-place only when level of fill is zero, otherwise allocates 4862 new space to store factored matrix and deletes previous memory. 4863 4864 Most users should employ the simplified KSP interface for linear solvers 4865 instead of working directly with matrix algebra routines such as this. 4866 See, e.g., KSPCreate(). 4867 4868 Level: developer 4869 4870 Concepts: matrices^incomplete Cholesky factorization 4871 Concepts: Cholesky factorization 4872 4873 .seealso: MatICCFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 4874 @*/ 4875 int MatICCFactor(Mat mat,IS row,MatFactorInfo* info) 4876 { 4877 int ierr; 4878 4879 PetscFunctionBegin; 4880 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4881 PetscValidType(mat); 4882 MatPreallocated(mat); 4883 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 4884 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4885 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4886 if (!mat->ops->iccfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4887 ierr = (*mat->ops->iccfactor)(mat,row,info);CHKERRQ(ierr); 4888 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 4889 PetscFunctionReturn(0); 4890 } 4891 4892 #undef __FUNCT__ 4893 #define __FUNCT__ "MatSetValuesAdic" 4894 /*@ 4895 MatSetValuesAdic - Sets values computed with ADIC automatic differentiation into a matrix. 4896 4897 Not Collective 4898 4899 Input Parameters: 4900 + mat - the matrix 4901 - v - the values compute with ADIC 4902 4903 Level: developer 4904 4905 Notes: 4906 Must call MatSetColoring() before using this routine. Also this matrix must already 4907 have its nonzero pattern determined. 4908 4909 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4910 MatSetValues(), MatSetColoring(), MatSetValuesAdifor() 4911 @*/ 4912 int MatSetValuesAdic(Mat mat,void *v) 4913 { 4914 int ierr; 4915 4916 PetscFunctionBegin; 4917 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4918 PetscValidType(mat); 4919 4920 if (!mat->assembled) { 4921 SETERRQ(1,"Matrix must be already assembled"); 4922 } 4923 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4924 if (!mat->ops->setvaluesadic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4925 ierr = (*mat->ops->setvaluesadic)(mat,v);CHKERRQ(ierr); 4926 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4927 ierr = MatView_Private(mat);CHKERRQ(ierr); 4928 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 4929 PetscFunctionReturn(0); 4930 } 4931 4932 4933 #undef __FUNCT__ 4934 #define __FUNCT__ "MatSetColoring" 4935 /*@ 4936 MatSetColoring - Sets a coloring used by calls to MatSetValuesAdic() 4937 4938 Not Collective 4939 4940 Input Parameters: 4941 + mat - the matrix 4942 - coloring - the coloring 4943 4944 Level: developer 4945 4946 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4947 MatSetValues(), MatSetValuesAdic() 4948 @*/ 4949 int MatSetColoring(Mat mat,ISColoring coloring) 4950 { 4951 int ierr; 4952 4953 PetscFunctionBegin; 4954 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4955 PetscValidType(mat); 4956 4957 if (!mat->assembled) { 4958 SETERRQ(1,"Matrix must be already assembled"); 4959 } 4960 if (!mat->ops->setcoloring) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4961 ierr = (*mat->ops->setcoloring)(mat,coloring);CHKERRQ(ierr); 4962 PetscFunctionReturn(0); 4963 } 4964 4965 #undef __FUNCT__ 4966 #define __FUNCT__ "MatSetValuesAdifor" 4967 /*@ 4968 MatSetValuesAdifor - Sets values computed with automatic differentiation into a matrix. 4969 4970 Not Collective 4971 4972 Input Parameters: 4973 + mat - the matrix 4974 . nl - leading dimension of v 4975 - v - the values compute with ADIFOR 4976 4977 Level: developer 4978 4979 Notes: 4980 Must call MatSetColoring() before using this routine. Also this matrix must already 4981 have its nonzero pattern determined. 4982 4983 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4984 MatSetValues(), MatSetColoring() 4985 @*/ 4986 int MatSetValuesAdifor(Mat mat,int nl,void *v) 4987 { 4988 int ierr; 4989 4990 PetscFunctionBegin; 4991 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4992 PetscValidType(mat); 4993 4994 if (!mat->assembled) { 4995 SETERRQ(1,"Matrix must be already assembled"); 4996 } 4997 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4998 if (!mat->ops->setvaluesadifor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4999 ierr = (*mat->ops->setvaluesadifor)(mat,nl,v);CHKERRQ(ierr); 5000 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 5001 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 5002 PetscFunctionReturn(0); 5003 } 5004 5005 EXTERN int MatMPIAIJDiagonalScaleLocal(Mat,Vec); 5006 EXTERN int MatMPIBAIJDiagonalScaleLocal(Mat,Vec); 5007 5008 #undef __FUNCT__ 5009 #define __FUNCT__ "MatDiagonalScaleLocal" 5010 /*@ 5011 MatDiagonalScaleLocal - Scales columns of a matrix given the scaling values including the 5012 ghosted ones. 5013 5014 Not Collective 5015 5016 Input Parameters: 5017 + mat - the matrix 5018 - diag = the diagonal values, including ghost ones 5019 5020 Level: developer 5021 5022 Notes: Works only for MPIAIJ and MPIBAIJ matrices 5023 5024 .seealso: MatDiagonalScale() 5025 @*/ 5026 int MatDiagonalScaleLocal(Mat mat,Vec diag) 5027 { 5028 int ierr,size; 5029 5030 PetscFunctionBegin; 5031 PetscValidHeaderSpecific(mat,MAT_COOKIE); 5032 PetscValidHeaderSpecific(diag,VEC_COOKIE); 5033 PetscValidType(mat); 5034 5035 if (!mat->assembled) { 5036 SETERRQ(1,"Matrix must be already assembled"); 5037 } 5038 ierr = PetscLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 5039 ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr); 5040 if (size == 1) { 5041 int n,m; 5042 ierr = VecGetSize(diag,&n);CHKERRQ(ierr); 5043 ierr = MatGetSize(mat,0,&m);CHKERRQ(ierr); 5044 if (m == n) { 5045 ierr = MatDiagonalScale(mat,0,diag);CHKERRQ(ierr); 5046 } else { 5047 SETERRQ(1,"Only supprted for sequential matrices when no ghost points/periodic conditions"); 5048 } 5049 } else { 5050 int (*f)(Mat,Vec); 5051 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",(void (**)(void))&f);CHKERRQ(ierr); 5052 if (f) { 5053 ierr = (*f)(mat,diag);CHKERRQ(ierr); 5054 } else { 5055 SETERRQ(1,"Only supported for MPIAIJ and MPIBAIJ parallel matrices"); 5056 } 5057 } 5058 ierr = PetscLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 5059 ierr = PetscObjectIncreaseState((PetscObject)mat); CHKERRQ(ierr); 5060 PetscFunctionReturn(0); 5061 } 5062 5063 #undef __FUNCT__ 5064 #define __FUNCT__ "MatGetInertia" 5065 /*@ 5066 MatGetInertia - Gets the inertia from a factored matrix 5067 5068 Collective on Mat 5069 5070 Input Parameter: 5071 . mat - the matrix 5072 5073 Output Parameters: 5074 + nneg - number of negative eigenvalues 5075 . nzero - number of zero eigenvalues 5076 - npos - number of positive eigenvalues 5077 5078 Level: advanced 5079 5080 Notes: Matrix must have been factored by MatCholeskyFactor() 5081 5082 5083 @*/ 5084 int MatGetInertia(Mat mat,int *nneg,int *nzero,int *npos) 5085 { 5086 int ierr; 5087 5088 PetscFunctionBegin; 5089 PetscValidHeaderSpecific(mat,MAT_COOKIE); 5090 PetscValidType(mat); 5091 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 5092 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Numeric factor mat is not assembled"); 5093 if (!mat->ops->getinertia) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 5094 ierr = (*mat->ops->getinertia)(mat,nneg,nzero,npos);CHKERRQ(ierr); 5095 PetscFunctionReturn(0); 5096 } 5097 5098 /* ----------------------------------------------------------------*/ 5099 #undef __FUNCT__ 5100 #define __FUNCT__ "MatSolves" 5101 /*@ 5102 MatSolves - Solves A x = b, given a factored matrix, for a collection of vectors 5103 5104 Collective on Mat and Vecs 5105 5106 Input Parameters: 5107 + mat - the factored matrix 5108 - b - the right-hand-side vectors 5109 5110 Output Parameter: 5111 . x - the result vectors 5112 5113 Notes: 5114 The vectors b and x cannot be the same. I.e., one cannot 5115 call MatSolves(A,x,x). 5116 5117 Notes: 5118 Most users should employ the simplified KSP interface for linear solvers 5119 instead of working directly with matrix algebra routines such as this. 5120 See, e.g., KSPCreate(). 5121 5122 Level: developer 5123 5124 Concepts: matrices^triangular solves 5125 5126 .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd(), MatSolve() 5127 @*/ 5128 int MatSolves(Mat mat,Vecs b,Vecs x) 5129 { 5130 int ierr; 5131 5132 PetscFunctionBegin; 5133 PetscValidHeaderSpecific(mat,MAT_COOKIE); 5134 PetscValidType(mat); 5135 MatPreallocated(mat); 5136 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 5137 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 5138 if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0); 5139 5140 if (!mat->ops->solves) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 5141 ierr = PetscLogEventBegin(MAT_Solves,mat,0,0,0);CHKERRQ(ierr); 5142 ierr = (*mat->ops->solves)(mat,b,x);CHKERRQ(ierr); 5143 ierr = PetscLogEventEnd(MAT_Solves,mat,0,0,0);CHKERRQ(ierr); 5144 PetscFunctionReturn(0); 5145 } 5146 5147 #undef __FUNCT__ 5148 #define __FUNCT__ "MatIsSymmetric" 5149 /*@C 5150 MatIsSymmetric - Test whether a matrix is symmetric 5151 5152 Collective on Mat 5153 5154 Input Parameter: 5155 . A - the matrix to test 5156 5157 Output Parameters: 5158 . flg - the result 5159 5160 Level: intermediate 5161 5162 Concepts: matrix^symmetry 5163 5164 .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsStructurallySymmetric(), MatSetOption() 5165 @*/ 5166 int MatIsSymmetric(Mat A,PetscTruth *flg) 5167 { 5168 int ierr; 5169 5170 PetscFunctionBegin; 5171 PetscValidHeaderSpecific(A,MAT_COOKIE); 5172 if (!A->symmetric_set) { 5173 if (!A->ops->issymmetric) SETERRQ(1,"Matrix does not support checking for symmetric"); 5174 ierr = (*A->ops->issymmetric)(A,&A->symmetric);CHKERRQ(ierr); 5175 A->symmetric_set = PETSC_TRUE; 5176 if (A->symmetric) { 5177 A->structurally_symmetric_set = PETSC_TRUE; 5178 A->structurally_symmetric = PETSC_TRUE; 5179 } 5180 } 5181 *flg = A->symmetric; 5182 PetscFunctionReturn(0); 5183 } 5184 5185 #undef __FUNCT__ 5186 #define __FUNCT__ "MatIsStructurallySymmetric" 5187 /*@C 5188 MatIsStructurallySymmetric - Test whether a matrix is structurally symmetric 5189 5190 Collective on Mat 5191 5192 Input Parameter: 5193 . A - the matrix to test 5194 5195 Output Parameters: 5196 . flg - the result 5197 5198 Level: intermediate 5199 5200 Concepts: matrix^symmetry 5201 5202 .seealso: MatTranspose(), MatIsTranspose(), MatIsHermitian(), MatIsSymmetric(), MatSetOption() 5203 @*/ 5204 int MatIsStructurallySymmetric(Mat A,PetscTruth *flg) 5205 { 5206 int ierr; 5207 5208 PetscFunctionBegin; 5209 PetscValidHeaderSpecific(A,MAT_COOKIE); 5210 if (!A->structurally_symmetric_set) { 5211 if (!A->ops->isstructurallysymmetric) SETERRQ(1,"Matrix does not support checking for structural symmetric"); 5212 ierr = (*A->ops->isstructurallysymmetric)(A,&A->structurally_symmetric);CHKERRQ(ierr); 5213 A->structurally_symmetric_set = PETSC_TRUE; 5214 } 5215 *flg = A->structurally_symmetric; 5216 PetscFunctionReturn(0); 5217 } 5218 5219 #undef __FUNCT__ 5220 #define __FUNCT__ "MatIsHermitian" 5221 /*@C 5222 MatIsHermitian - Test whether a matrix is Hermitian, i.e. it is the complex conjugate of its transpose. 5223 5224 Collective on Mat 5225 5226 Input Parameter: 5227 . A - the matrix to test 5228 5229 Output Parameters: 5230 . flg - the result 5231 5232 Level: intermediate 5233 5234 Concepts: matrix^symmetry 5235 5236 .seealso: MatTranspose(), MatIsTranspose(), MatIsSymmetric(), MatIsStructurallySymmetric(), MatSetOption() 5237 @*/ 5238 int MatIsHermitian(Mat A,PetscTruth *flg) 5239 { 5240 int ierr; 5241 5242 PetscFunctionBegin; 5243 PetscValidHeaderSpecific(A,MAT_COOKIE); 5244 if (!A->hermitian_set) { 5245 if (!A->ops->ishermitian) SETERRQ(1,"Matrix does not support checking for being Hermitian"); 5246 ierr = (*A->ops->ishermitian)(A,&A->hermitian);CHKERRQ(ierr); 5247 A->hermitian_set = PETSC_TRUE; 5248 if (A->hermitian) { 5249 A->structurally_symmetric_set = PETSC_TRUE; 5250 A->structurally_symmetric = PETSC_TRUE; 5251 } 5252 } 5253 *flg = A->hermitian; 5254 PetscFunctionReturn(0); 5255 } 5256