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