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