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