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