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