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