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 if (!mat->ops->multtranspose) SETERRQ(PETSC_ERR_SUP,"This matrix type does not have a multiply tranpose defined"); 1173 ierr = (*mat->ops->multtranspose)(mat,x,y);CHKERRQ(ierr); 1174 ierr = PetscLogEventEnd(MAT_MultTranspose,mat,x,y,0);CHKERRQ(ierr); 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "MatMultAdd" 1180 /*@ 1181 MatMultAdd - Computes v3 = v2 + A * v1. 1182 1183 Collective on Mat and Vec 1184 1185 Input Parameters: 1186 + mat - the matrix 1187 - v1, v2 - the vectors 1188 1189 Output Parameters: 1190 . v3 - the result 1191 1192 Notes: 1193 The vectors v1 and v3 cannot be the same. I.e., one cannot 1194 call MatMultAdd(A,v1,v2,v1). 1195 1196 Level: beginner 1197 1198 Concepts: matrix vector product^addition 1199 1200 .seealso: MatMultTranspose(), MatMult(), MatMultTransposeAdd() 1201 @*/ 1202 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3) 1203 { 1204 int ierr; 1205 1206 PetscFunctionBegin; 1207 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1208 PetscValidType(mat); 1209 MatPreallocated(mat); 1210 PetscValidHeaderSpecific(v1,VEC_COOKIE); 1211 PetscValidHeaderSpecific(v2,VEC_COOKIE); 1212 PetscValidHeaderSpecific(v3,VEC_COOKIE); 1213 1214 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1215 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1216 if (mat->N != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->N,v1->N); 1217 if (mat->M != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->M,v2->N); 1218 if (mat->M != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->M,v3->N); 1219 if (mat->m != v3->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: local dim %d %d",mat->m,v3->n); 1220 if (mat->m != v2->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: local dim %d %d",mat->m,v2->n); 1221 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors"); 1222 1223 ierr = PetscLogEventBegin(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1224 ierr = (*mat->ops->multadd)(mat,v1,v2,v3);CHKERRQ(ierr); 1225 ierr = PetscLogEventEnd(MAT_MultAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1226 PetscFunctionReturn(0); 1227 } 1228 1229 #undef __FUNCT__ 1230 #define __FUNCT__ "MatMultTransposeAdd" 1231 /*@ 1232 MatMultTransposeAdd - Computes v3 = v2 + A' * v1. 1233 1234 Collective on Mat and Vec 1235 1236 Input Parameters: 1237 + mat - the matrix 1238 - v1, v2 - the vectors 1239 1240 Output Parameters: 1241 . v3 - the result 1242 1243 Notes: 1244 The vectors v1 and v3 cannot be the same. I.e., one cannot 1245 call MatMultTransposeAdd(A,v1,v2,v1). 1246 1247 Level: beginner 1248 1249 Concepts: matrix vector product^transpose and addition 1250 1251 .seealso: MatMultTranspose(), MatMultAdd(), MatMult() 1252 @*/ 1253 int MatMultTransposeAdd(Mat mat,Vec v1,Vec v2,Vec v3) 1254 { 1255 int ierr; 1256 1257 PetscFunctionBegin; 1258 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1259 PetscValidType(mat); 1260 MatPreallocated(mat); 1261 PetscValidHeaderSpecific(v1,VEC_COOKIE); 1262 PetscValidHeaderSpecific(v2,VEC_COOKIE); 1263 PetscValidHeaderSpecific(v3,VEC_COOKIE); 1264 1265 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1266 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1267 if (!mat->ops->multtransposeadd) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1268 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,"v1 and v3 must be different vectors"); 1269 if (mat->M != v1->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v1: global dim %d %d",mat->M,v1->N); 1270 if (mat->N != v2->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v2: global dim %d %d",mat->N,v2->N); 1271 if (mat->N != v3->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec v3: global dim %d %d",mat->N,v3->N); 1272 1273 ierr = PetscLogEventBegin(MAT_MultTransposeAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1274 ierr = (*mat->ops->multtransposeadd)(mat,v1,v2,v3);CHKERRQ(ierr); 1275 ierr = PetscLogEventEnd(MAT_MultTransposeAdd,mat,v1,v2,v3);CHKERRQ(ierr); 1276 PetscFunctionReturn(0); 1277 } 1278 1279 #undef __FUNCT__ 1280 #define __FUNCT__ "MatMultConstrained" 1281 /*@ 1282 MatMultConstrained - The inner multiplication routine for a 1283 constrained matrix P^T A P. 1284 1285 Collective on Mat and Vec 1286 1287 Input Parameters: 1288 + mat - the matrix 1289 - x - the vector to be multilplied 1290 1291 Output Parameters: 1292 . y - the result 1293 1294 Notes: 1295 The vectors x and y cannot be the same. I.e., one cannot 1296 call MatMult(A,y,y). 1297 1298 Level: beginner 1299 1300 .keywords: matrix, multiply, matrix-vector product, constraint 1301 .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd() 1302 @*/ 1303 int MatMultConstrained(Mat mat,Vec x,Vec y) 1304 { 1305 int ierr; 1306 1307 PetscFunctionBegin; 1308 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1309 PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1310 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1311 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1312 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1313 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1314 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 1315 if (mat->m != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: local dim %d %d",mat->m,y->n); 1316 1317 ierr = PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1318 ierr = (*mat->ops->multconstrained)(mat,x,y); CHKERRQ(ierr); 1319 ierr = PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1320 1321 PetscFunctionReturn(0); 1322 } 1323 1324 #undef __FUNCT__ 1325 #define __FUNCT__ "MatMultTransposeConstrained" 1326 /*@ 1327 MatMultTransposeConstrained - The inner multiplication routine for a 1328 constrained matrix P^T A^T P. 1329 1330 Collective on Mat and Vec 1331 1332 Input Parameters: 1333 + mat - the matrix 1334 - x - the vector to be multilplied 1335 1336 Output Parameters: 1337 . y - the result 1338 1339 Notes: 1340 The vectors x and y cannot be the same. I.e., one cannot 1341 call MatMult(A,y,y). 1342 1343 Level: beginner 1344 1345 .keywords: matrix, multiply, matrix-vector product, constraint 1346 .seealso: MatMult(), MatMultTrans(), MatMultAdd(), MatMultTransAdd() 1347 @*/ 1348 int MatMultTransposeConstrained(Mat mat,Vec x,Vec y) 1349 { 1350 int ierr; 1351 1352 PetscFunctionBegin; 1353 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1354 PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1355 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1356 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1357 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"x and y must be different vectors"); 1358 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1359 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 1360 1361 ierr = PetscLogEventBegin(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1362 ierr = (*mat->ops->multtransposeconstrained)(mat,x,y);CHKERRQ(ierr); 1363 ierr = PetscLogEventEnd(MAT_MultConstrained,mat,x,y,0);CHKERRQ(ierr); 1364 1365 PetscFunctionReturn(0); 1366 } 1367 /* ------------------------------------------------------------*/ 1368 #undef __FUNCT__ 1369 #define __FUNCT__ "MatGetInfo" 1370 /*@C 1371 MatGetInfo - Returns information about matrix storage (number of 1372 nonzeros, memory, etc.). 1373 1374 Collective on Mat if MAT_GLOBAL_MAX or MAT_GLOBAL_SUM is used 1375 as the flag 1376 1377 Input Parameters: 1378 . mat - the matrix 1379 1380 Output Parameters: 1381 + flag - flag indicating the type of parameters to be returned 1382 (MAT_LOCAL - local matrix, MAT_GLOBAL_MAX - maximum over all processors, 1383 MAT_GLOBAL_SUM - sum over all processors) 1384 - info - matrix information context 1385 1386 Notes: 1387 The MatInfo context contains a variety of matrix data, including 1388 number of nonzeros allocated and used, number of mallocs during 1389 matrix assembly, etc. Additional information for factored matrices 1390 is provided (such as the fill ratio, number of mallocs during 1391 factorization, etc.). Much of this info is printed to STDOUT 1392 when using the runtime options 1393 $ -log_info -mat_view_info 1394 1395 Example for C/C++ Users: 1396 See the file ${PETSC_DIR}/include/petscmat.h for a complete list of 1397 data within the MatInfo context. For example, 1398 .vb 1399 MatInfo info; 1400 Mat A; 1401 double mal, nz_a, nz_u; 1402 1403 MatGetInfo(A,MAT_LOCAL,&info); 1404 mal = info.mallocs; 1405 nz_a = info.nz_allocated; 1406 .ve 1407 1408 Example for Fortran Users: 1409 Fortran users should declare info as a double precision 1410 array of dimension MAT_INFO_SIZE, and then extract the parameters 1411 of interest. See the file ${PETSC_DIR}/include/finclude/petscmat.h 1412 a complete list of parameter names. 1413 .vb 1414 double precision info(MAT_INFO_SIZE) 1415 double precision mal, nz_a 1416 Mat A 1417 integer ierr 1418 1419 call MatGetInfo(A,MAT_LOCAL,info,ierr) 1420 mal = info(MAT_INFO_MALLOCS) 1421 nz_a = info(MAT_INFO_NZ_ALLOCATED) 1422 .ve 1423 1424 Level: intermediate 1425 1426 Concepts: matrices^getting information on 1427 1428 @*/ 1429 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info) 1430 { 1431 int ierr; 1432 1433 PetscFunctionBegin; 1434 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1435 PetscValidType(mat); 1436 MatPreallocated(mat); 1437 PetscValidPointer(info); 1438 if (!mat->ops->getinfo) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1439 ierr = (*mat->ops->getinfo)(mat,flag,info);CHKERRQ(ierr); 1440 PetscFunctionReturn(0); 1441 } 1442 1443 /* ----------------------------------------------------------*/ 1444 #undef __FUNCT__ 1445 #define __FUNCT__ "MatILUDTFactor" 1446 /*@C 1447 MatILUDTFactor - Performs a drop tolerance ILU factorization. 1448 1449 Collective on Mat 1450 1451 Input Parameters: 1452 + mat - the matrix 1453 . info - information about the factorization to be done 1454 . row - row permutation 1455 - col - column permutation 1456 1457 Output Parameters: 1458 . fact - the factored matrix 1459 1460 Level: developer 1461 1462 Notes: 1463 Most users should employ the simplified SLES interface for linear solvers 1464 instead of working directly with matrix algebra routines such as this. 1465 See, e.g., SLESCreate(). 1466 1467 This is currently only supported for the SeqAIJ matrix format using code 1468 from Yousef Saad's SPARSEKIT2 package (translated to C with f2c) and/or 1469 Matlab. SPARSEKIT2 is copyrighted by Yousef Saad with the GNU copyright 1470 and thus can be distributed with PETSc. 1471 1472 Concepts: matrices^ILUDT factorization 1473 1474 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatILUInfo 1475 @*/ 1476 int MatILUDTFactor(Mat mat,MatILUInfo *info,IS row,IS col,Mat *fact) 1477 { 1478 int ierr; 1479 1480 PetscFunctionBegin; 1481 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1482 PetscValidType(mat); 1483 MatPreallocated(mat); 1484 PetscValidPointer(fact); 1485 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1486 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1487 if (!mat->ops->iludtfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1488 1489 ierr = PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1490 ierr = (*mat->ops->iludtfactor)(mat,info,row,col,fact);CHKERRQ(ierr); 1491 ierr = PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1492 1493 PetscFunctionReturn(0); 1494 } 1495 1496 #undef __FUNCT__ 1497 #define __FUNCT__ "MatLUFactor" 1498 /*@ 1499 MatLUFactor - Performs in-place LU factorization of matrix. 1500 1501 Collective on Mat 1502 1503 Input Parameters: 1504 + mat - the matrix 1505 . row - row permutation 1506 . col - column permutation 1507 - info - options for factorization, includes 1508 $ fill - expected fill as ratio of original fill. 1509 $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting) 1510 $ Run with the option -log_info to determine an optimal value to use 1511 1512 Notes: 1513 Most users should employ the simplified SLES interface for linear solvers 1514 instead of working directly with matrix algebra routines such as this. 1515 See, e.g., SLESCreate(). 1516 1517 This changes the state of the matrix to a factored matrix; it cannot be used 1518 for example with MatSetValues() unless one first calls MatSetUnfactored(). 1519 1520 Level: developer 1521 1522 Concepts: matrices^LU factorization 1523 1524 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), 1525 MatGetOrdering(), MatSetUnfactored(), MatLUInfo 1526 1527 @*/ 1528 int MatLUFactor(Mat mat,IS row,IS col,MatLUInfo *info) 1529 { 1530 int ierr; 1531 1532 PetscFunctionBegin; 1533 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1534 PetscValidType(mat); 1535 MatPreallocated(mat); 1536 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1537 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1538 if (!mat->ops->lufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1539 1540 ierr = PetscLogEventBegin(MAT_LUFactor,mat,row,col,0);CHKERRQ(ierr); 1541 ierr = (*mat->ops->lufactor)(mat,row,col,info);CHKERRQ(ierr); 1542 ierr = PetscLogEventEnd(MAT_LUFactor,mat,row,col,0);CHKERRQ(ierr); 1543 PetscFunctionReturn(0); 1544 } 1545 1546 #undef __FUNCT__ 1547 #define __FUNCT__ "MatILUFactor" 1548 /*@ 1549 MatILUFactor - Performs in-place ILU factorization of matrix. 1550 1551 Collective on Mat 1552 1553 Input Parameters: 1554 + mat - the matrix 1555 . row - row permutation 1556 . col - column permutation 1557 - info - structure containing 1558 $ levels - number of levels of fill. 1559 $ expected fill - as ratio of original fill. 1560 $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices 1561 missing diagonal entries) 1562 1563 Notes: 1564 Probably really in-place only when level of fill is zero, otherwise allocates 1565 new space to store factored matrix and deletes previous memory. 1566 1567 Most users should employ the simplified SLES interface for linear solvers 1568 instead of working directly with matrix algebra routines such as this. 1569 See, e.g., SLESCreate(). 1570 1571 Level: developer 1572 1573 Concepts: matrices^ILU factorization 1574 1575 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor(), MatILUInfo 1576 @*/ 1577 int MatILUFactor(Mat mat,IS row,IS col,MatILUInfo *info) 1578 { 1579 int ierr; 1580 1581 PetscFunctionBegin; 1582 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1583 PetscValidType(mat); 1584 MatPreallocated(mat); 1585 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 1586 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1587 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1588 if (!mat->ops->ilufactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1589 1590 ierr = PetscLogEventBegin(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1591 ierr = (*mat->ops->ilufactor)(mat,row,col,info);CHKERRQ(ierr); 1592 ierr = PetscLogEventEnd(MAT_ILUFactor,mat,row,col,0);CHKERRQ(ierr); 1593 PetscFunctionReturn(0); 1594 } 1595 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "MatLUFactorSymbolic" 1598 /*@ 1599 MatLUFactorSymbolic - Performs symbolic LU factorization of matrix. 1600 Call this routine before calling MatLUFactorNumeric(). 1601 1602 Collective on Mat 1603 1604 Input Parameters: 1605 + mat - the matrix 1606 . row, col - row and column permutations 1607 - info - options for factorization, includes 1608 $ fill - expected fill as ratio of original fill. 1609 $ dtcol - pivot tolerance (0 no pivot, 1 full column pivoting) 1610 $ Run with the option -log_info to determine an optimal value to use 1611 1612 Output Parameter: 1613 . fact - new matrix that has been symbolically factored 1614 1615 Notes: 1616 See the users manual for additional information about 1617 choosing the fill factor for better efficiency. 1618 1619 Most users should employ the simplified SLES interface for linear solvers 1620 instead of working directly with matrix algebra routines such as this. 1621 See, e.g., SLESCreate(). 1622 1623 Level: developer 1624 1625 Concepts: matrices^LU symbolic factorization 1626 1627 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor(), MatLUInfo 1628 @*/ 1629 int MatLUFactorSymbolic(Mat mat,IS row,IS col,MatLUInfo *info,Mat *fact) 1630 { 1631 int ierr; 1632 1633 PetscFunctionBegin; 1634 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1635 PetscValidType(mat); 1636 MatPreallocated(mat); 1637 PetscValidPointer(fact); 1638 PetscValidHeaderSpecific(row,IS_COOKIE); 1639 PetscValidHeaderSpecific(col,IS_COOKIE); 1640 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1641 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1642 if (!mat->ops->lufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic LU",mat->type_name); 1643 1644 ierr = PetscLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 1645 ierr = (*mat->ops->lufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr); 1646 ierr = PetscLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 1647 PetscFunctionReturn(0); 1648 } 1649 1650 #undef __FUNCT__ 1651 #define __FUNCT__ "MatLUFactorNumeric" 1652 /*@ 1653 MatLUFactorNumeric - Performs numeric LU factorization of a matrix. 1654 Call this routine after first calling MatLUFactorSymbolic(). 1655 1656 Collective on Mat 1657 1658 Input Parameters: 1659 + mat - the matrix 1660 - fact - the matrix generated for the factor, from MatLUFactorSymbolic() 1661 1662 Notes: 1663 See MatLUFactor() for in-place factorization. See 1664 MatCholeskyFactorNumeric() for the symmetric, positive definite case. 1665 1666 Most users should employ the simplified SLES interface for linear solvers 1667 instead of working directly with matrix algebra routines such as this. 1668 See, e.g., SLESCreate(). 1669 1670 Level: developer 1671 1672 Concepts: matrices^LU numeric factorization 1673 1674 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor() 1675 @*/ 1676 int MatLUFactorNumeric(Mat mat,Mat *fact) 1677 { 1678 int ierr; 1679 PetscTruth flg; 1680 1681 PetscFunctionBegin; 1682 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1683 PetscValidType(mat); 1684 MatPreallocated(mat); 1685 PetscValidPointer(fact); 1686 PetscValidHeaderSpecific(*fact,MAT_COOKIE); 1687 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1688 if (mat->M != (*fact)->M || mat->N != (*fact)->N) { 1689 SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dimensions are different %d should = %d %d should = %d", 1690 mat->M,(*fact)->M,mat->N,(*fact)->N); 1691 } 1692 if (!(*fact)->ops->lufactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1693 1694 ierr = PetscLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1695 ierr = (*(*fact)->ops->lufactornumeric)(mat,fact);CHKERRQ(ierr); 1696 ierr = PetscLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1697 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_view_draw",&flg);CHKERRQ(ierr); 1698 if (flg) { 1699 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_view_contour",&flg);CHKERRQ(ierr); 1700 if (flg) { 1701 ierr = PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 1702 } 1703 ierr = MatView(*fact,PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 1704 ierr = PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 1705 if (flg) { 1706 ierr = PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 1707 } 1708 } 1709 PetscFunctionReturn(0); 1710 } 1711 1712 #undef __FUNCT__ 1713 #define __FUNCT__ "MatCholeskyFactor" 1714 /*@ 1715 MatCholeskyFactor - Performs in-place Cholesky factorization of a 1716 symmetric matrix. 1717 1718 Collective on Mat 1719 1720 Input Parameters: 1721 + mat - the matrix 1722 . perm - row and column permutations 1723 - f - expected fill as ratio of original fill 1724 1725 Notes: 1726 See MatLUFactor() for the nonsymmetric case. See also 1727 MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric(). 1728 1729 Most users should employ the simplified SLES interface for linear solvers 1730 instead of working directly with matrix algebra routines such as this. 1731 See, e.g., SLESCreate(). 1732 1733 Level: developer 1734 1735 Concepts: matrices^Cholesky factorization 1736 1737 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric() 1738 MatGetOrdering() 1739 1740 @*/ 1741 int MatCholeskyFactor(Mat mat,IS perm,PetscReal f) 1742 { 1743 int ierr; 1744 1745 PetscFunctionBegin; 1746 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1747 PetscValidType(mat); 1748 MatPreallocated(mat); 1749 PetscValidHeaderSpecific(perm,IS_COOKIE); 1750 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 1751 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1752 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1753 if (!mat->ops->choleskyfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1754 1755 ierr = PetscLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0);CHKERRQ(ierr); 1756 ierr = (*mat->ops->choleskyfactor)(mat,perm,f);CHKERRQ(ierr); 1757 ierr = PetscLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0);CHKERRQ(ierr); 1758 PetscFunctionReturn(0); 1759 } 1760 1761 #undef __FUNCT__ 1762 #define __FUNCT__ "MatCholeskyFactorSymbolic" 1763 /*@ 1764 MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization 1765 of a symmetric matrix. 1766 1767 Collective on Mat 1768 1769 Input Parameters: 1770 + mat - the matrix 1771 . perm - row and column permutations 1772 - f - expected fill as ratio of original 1773 1774 Output Parameter: 1775 . fact - the factored matrix 1776 1777 Notes: 1778 See MatLUFactorSymbolic() for the nonsymmetric case. See also 1779 MatCholeskyFactor() and MatCholeskyFactorNumeric(). 1780 1781 Most users should employ the simplified SLES interface for linear solvers 1782 instead of working directly with matrix algebra routines such as this. 1783 See, e.g., SLESCreate(). 1784 1785 Level: developer 1786 1787 Concepts: matrices^Cholesky symbolic factorization 1788 1789 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric() 1790 MatGetOrdering() 1791 1792 @*/ 1793 int MatCholeskyFactorSymbolic(Mat mat,IS perm,PetscReal f,Mat *fact) 1794 { 1795 int ierr; 1796 1797 PetscFunctionBegin; 1798 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1799 PetscValidType(mat); 1800 MatPreallocated(mat); 1801 PetscValidPointer(fact); 1802 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"Matrix must be square"); 1803 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1804 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1805 if (!mat->ops->choleskyfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1806 1807 ierr = PetscLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 1808 ierr = (*mat->ops->choleskyfactorsymbolic)(mat,perm,f,fact);CHKERRQ(ierr); 1809 ierr = PetscLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 1810 PetscFunctionReturn(0); 1811 } 1812 1813 #undef __FUNCT__ 1814 #define __FUNCT__ "MatCholeskyFactorNumeric" 1815 /*@ 1816 MatCholeskyFactorNumeric - Performs numeric Cholesky factorization 1817 of a symmetric matrix. Call this routine after first calling 1818 MatCholeskyFactorSymbolic(). 1819 1820 Collective on Mat 1821 1822 Input Parameter: 1823 . mat - the initial matrix 1824 1825 Output Parameter: 1826 . fact - the factored matrix 1827 1828 Notes: 1829 Most users should employ the simplified SLES interface for linear solvers 1830 instead of working directly with matrix algebra routines such as this. 1831 See, e.g., SLESCreate(). 1832 1833 Level: developer 1834 1835 Concepts: matrices^Cholesky numeric factorization 1836 1837 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric() 1838 @*/ 1839 int MatCholeskyFactorNumeric(Mat mat,Mat *fact) 1840 { 1841 int ierr; 1842 1843 PetscFunctionBegin; 1844 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1845 PetscValidType(mat); 1846 MatPreallocated(mat); 1847 PetscValidPointer(fact); 1848 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1849 if (!(*fact)->ops->choleskyfactornumeric) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1850 if (mat->M != (*fact)->M || mat->N != (*fact)->N) { 1851 SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat mat,Mat *fact: global dim %d should = %d %d should = %d", 1852 mat->M,(*fact)->M,mat->N,(*fact)->N); 1853 } 1854 1855 ierr = PetscLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1856 ierr = (*(*fact)->ops->choleskyfactornumeric)(mat,fact);CHKERRQ(ierr); 1857 ierr = PetscLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0);CHKERRQ(ierr); 1858 PetscFunctionReturn(0); 1859 } 1860 1861 /* ----------------------------------------------------------------*/ 1862 #undef __FUNCT__ 1863 #define __FUNCT__ "MatSolve" 1864 /*@ 1865 MatSolve - Solves A x = b, given a factored matrix. 1866 1867 Collective on Mat and Vec 1868 1869 Input Parameters: 1870 + mat - the factored matrix 1871 - b - the right-hand-side vector 1872 1873 Output Parameter: 1874 . x - the result vector 1875 1876 Notes: 1877 The vectors b and x cannot be the same. I.e., one cannot 1878 call MatSolve(A,x,x). 1879 1880 Notes: 1881 Most users should employ the simplified SLES interface for linear solvers 1882 instead of working directly with matrix algebra routines such as this. 1883 See, e.g., SLESCreate(). 1884 1885 Level: developer 1886 1887 Concepts: matrices^triangular solves 1888 1889 .seealso: MatSolveAdd(), MatSolveTranspose(), MatSolveTransposeAdd() 1890 @*/ 1891 int MatSolve(Mat mat,Vec b,Vec x) 1892 { 1893 int ierr; 1894 1895 PetscFunctionBegin; 1896 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1897 PetscValidType(mat); 1898 MatPreallocated(mat); 1899 PetscValidHeaderSpecific(b,VEC_COOKIE); 1900 PetscValidHeaderSpecific(x,VEC_COOKIE); 1901 PetscCheckSameComm(mat,b); 1902 PetscCheckSameComm(mat,x); 1903 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 1904 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 1905 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1906 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1907 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1908 if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0); 1909 1910 if (!mat->ops->solve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1911 ierr = PetscLogEventBegin(MAT_Solve,mat,b,x,0);CHKERRQ(ierr); 1912 ierr = (*mat->ops->solve)(mat,b,x);CHKERRQ(ierr); 1913 ierr = PetscLogEventEnd(MAT_Solve,mat,b,x,0);CHKERRQ(ierr); 1914 PetscFunctionReturn(0); 1915 } 1916 1917 #undef __FUNCT__ 1918 #define __FUNCT__ "MatForwardSolve" 1919 /* @ 1920 MatForwardSolve - Solves L x = b, given a factored matrix, A = LU. 1921 1922 Collective on Mat and Vec 1923 1924 Input Parameters: 1925 + mat - the factored matrix 1926 - b - the right-hand-side vector 1927 1928 Output Parameter: 1929 . x - the result vector 1930 1931 Notes: 1932 MatSolve() should be used for most applications, as it performs 1933 a forward solve followed by a backward solve. 1934 1935 The vectors b and x cannot be the same. I.e., one cannot 1936 call MatForwardSolve(A,x,x). 1937 1938 Most users should employ the simplified SLES interface for linear solvers 1939 instead of working directly with matrix algebra routines such as this. 1940 See, e.g., SLESCreate(). 1941 1942 Level: developer 1943 1944 Concepts: matrices^forward solves 1945 1946 .seealso: MatSolve(), MatBackwardSolve() 1947 @ */ 1948 int MatForwardSolve(Mat mat,Vec b,Vec x) 1949 { 1950 int ierr; 1951 1952 PetscFunctionBegin; 1953 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1954 PetscValidType(mat); 1955 MatPreallocated(mat); 1956 PetscValidHeaderSpecific(b,VEC_COOKIE); 1957 PetscValidHeaderSpecific(x,VEC_COOKIE); 1958 PetscCheckSameComm(mat,b); 1959 PetscCheckSameComm(mat,x); 1960 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 1961 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 1962 if (!mat->ops->forwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 1963 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 1964 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 1965 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 1966 1967 ierr = PetscLogEventBegin(MAT_ForwardSolve,mat,b,x,0);CHKERRQ(ierr); 1968 ierr = (*mat->ops->forwardsolve)(mat,b,x);CHKERRQ(ierr); 1969 ierr = PetscLogEventEnd(MAT_ForwardSolve,mat,b,x,0);CHKERRQ(ierr); 1970 PetscFunctionReturn(0); 1971 } 1972 1973 #undef __FUNCT__ 1974 #define __FUNCT__ "MatBackwardSolve" 1975 /* @ 1976 MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU. 1977 1978 Collective on Mat and Vec 1979 1980 Input Parameters: 1981 + mat - the factored matrix 1982 - b - the right-hand-side vector 1983 1984 Output Parameter: 1985 . x - the result vector 1986 1987 Notes: 1988 MatSolve() should be used for most applications, as it performs 1989 a forward solve followed by a backward solve. 1990 1991 The vectors b and x cannot be the same. I.e., one cannot 1992 call MatBackwardSolve(A,x,x). 1993 1994 Most users should employ the simplified SLES interface for linear solvers 1995 instead of working directly with matrix algebra routines such as this. 1996 See, e.g., SLESCreate(). 1997 1998 Level: developer 1999 2000 Concepts: matrices^backward solves 2001 2002 .seealso: MatSolve(), MatForwardSolve() 2003 @ */ 2004 int MatBackwardSolve(Mat mat,Vec b,Vec x) 2005 { 2006 int ierr; 2007 2008 PetscFunctionBegin; 2009 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2010 PetscValidType(mat); 2011 MatPreallocated(mat); 2012 PetscValidHeaderSpecific(b,VEC_COOKIE); 2013 PetscValidHeaderSpecific(x,VEC_COOKIE); 2014 PetscCheckSameComm(mat,b); 2015 PetscCheckSameComm(mat,x); 2016 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2017 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2018 if (!mat->ops->backwardsolve) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2019 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2020 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2021 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2022 2023 ierr = PetscLogEventBegin(MAT_BackwardSolve,mat,b,x,0);CHKERRQ(ierr); 2024 ierr = (*mat->ops->backwardsolve)(mat,b,x);CHKERRQ(ierr); 2025 ierr = PetscLogEventEnd(MAT_BackwardSolve,mat,b,x,0);CHKERRQ(ierr); 2026 PetscFunctionReturn(0); 2027 } 2028 2029 #undef __FUNCT__ 2030 #define __FUNCT__ "MatSolveAdd" 2031 /*@ 2032 MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix. 2033 2034 Collective on Mat and Vec 2035 2036 Input Parameters: 2037 + mat - the factored matrix 2038 . b - the right-hand-side vector 2039 - y - the vector to be added to 2040 2041 Output Parameter: 2042 . x - the result vector 2043 2044 Notes: 2045 The vectors b and x cannot be the same. I.e., one cannot 2046 call MatSolveAdd(A,x,y,x). 2047 2048 Most users should employ the simplified SLES interface for linear solvers 2049 instead of working directly with matrix algebra routines such as this. 2050 See, e.g., SLESCreate(). 2051 2052 Level: developer 2053 2054 Concepts: matrices^triangular solves 2055 2056 .seealso: MatSolve(), MatSolveTranspose(), MatSolveTransposeAdd() 2057 @*/ 2058 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x) 2059 { 2060 PetscScalar one = 1.0; 2061 Vec tmp; 2062 int ierr; 2063 2064 PetscFunctionBegin; 2065 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2066 PetscValidType(mat); 2067 MatPreallocated(mat); 2068 PetscValidHeaderSpecific(y,VEC_COOKIE); 2069 PetscValidHeaderSpecific(b,VEC_COOKIE); 2070 PetscValidHeaderSpecific(x,VEC_COOKIE); 2071 PetscCheckSameComm(mat,b); 2072 PetscCheckSameComm(mat,y); 2073 PetscCheckSameComm(mat,x); 2074 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2075 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2076 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2077 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2078 if (mat->M != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->M,y->N); 2079 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2080 if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n); 2081 2082 ierr = PetscLogEventBegin(MAT_SolveAdd,mat,b,x,y);CHKERRQ(ierr); 2083 if (mat->ops->solveadd) { 2084 ierr = (*mat->ops->solveadd)(mat,b,y,x);CHKERRQ(ierr); 2085 } else { 2086 /* do the solve then the add manually */ 2087 if (x != y) { 2088 ierr = MatSolve(mat,b,x);CHKERRQ(ierr); 2089 ierr = VecAXPY(&one,y,x);CHKERRQ(ierr); 2090 } else { 2091 ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr); 2092 PetscLogObjectParent(mat,tmp); 2093 ierr = VecCopy(x,tmp);CHKERRQ(ierr); 2094 ierr = MatSolve(mat,b,x);CHKERRQ(ierr); 2095 ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr); 2096 ierr = VecDestroy(tmp);CHKERRQ(ierr); 2097 } 2098 } 2099 ierr = PetscLogEventEnd(MAT_SolveAdd,mat,b,x,y);CHKERRQ(ierr); 2100 PetscFunctionReturn(0); 2101 } 2102 2103 #undef __FUNCT__ 2104 #define __FUNCT__ "MatSolveTranspose" 2105 /*@ 2106 MatSolveTranspose - Solves A' x = b, given a factored matrix. 2107 2108 Collective on Mat and Vec 2109 2110 Input Parameters: 2111 + mat - the factored matrix 2112 - b - the right-hand-side vector 2113 2114 Output Parameter: 2115 . x - the result vector 2116 2117 Notes: 2118 The vectors b and x cannot be the same. I.e., one cannot 2119 call MatSolveTranspose(A,x,x). 2120 2121 Most users should employ the simplified SLES interface for linear solvers 2122 instead of working directly with matrix algebra routines such as this. 2123 See, e.g., SLESCreate(). 2124 2125 Level: developer 2126 2127 Concepts: matrices^triangular solves 2128 2129 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransposeAdd() 2130 @*/ 2131 int MatSolveTranspose(Mat mat,Vec b,Vec x) 2132 { 2133 int ierr; 2134 2135 PetscFunctionBegin; 2136 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2137 PetscValidType(mat); 2138 MatPreallocated(mat); 2139 PetscValidHeaderSpecific(b,VEC_COOKIE); 2140 PetscValidHeaderSpecific(x,VEC_COOKIE); 2141 PetscCheckSameComm(mat,b); 2142 PetscCheckSameComm(mat,x); 2143 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2144 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2145 if (!mat->ops->solvetranspose) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s",mat->type_name); 2146 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 2147 if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N); 2148 2149 ierr = PetscLogEventBegin(MAT_SolveTranspose,mat,b,x,0);CHKERRQ(ierr); 2150 ierr = (*mat->ops->solvetranspose)(mat,b,x);CHKERRQ(ierr); 2151 ierr = PetscLogEventEnd(MAT_SolveTranspose,mat,b,x,0);CHKERRQ(ierr); 2152 PetscFunctionReturn(0); 2153 } 2154 2155 #undef __FUNCT__ 2156 #define __FUNCT__ "MatSolveTransposeAdd" 2157 /*@ 2158 MatSolveTransposeAdd - Computes x = y + inv(Transpose(A)) b, given a 2159 factored matrix. 2160 2161 Collective on Mat and Vec 2162 2163 Input Parameters: 2164 + mat - the factored matrix 2165 . b - the right-hand-side vector 2166 - y - the vector to be added to 2167 2168 Output Parameter: 2169 . x - the result vector 2170 2171 Notes: 2172 The vectors b and x cannot be the same. I.e., one cannot 2173 call MatSolveTransposeAdd(A,x,y,x). 2174 2175 Most users should employ the simplified SLES interface for linear solvers 2176 instead of working directly with matrix algebra routines such as this. 2177 See, e.g., SLESCreate(). 2178 2179 Level: developer 2180 2181 Concepts: matrices^triangular solves 2182 2183 .seealso: MatSolve(), MatSolveAdd(), MatSolveTranspose() 2184 @*/ 2185 int MatSolveTransposeAdd(Mat mat,Vec b,Vec y,Vec x) 2186 { 2187 PetscScalar one = 1.0; 2188 int ierr; 2189 Vec tmp; 2190 2191 PetscFunctionBegin; 2192 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2193 PetscValidType(mat); 2194 MatPreallocated(mat); 2195 PetscValidHeaderSpecific(y,VEC_COOKIE); 2196 PetscValidHeaderSpecific(b,VEC_COOKIE); 2197 PetscValidHeaderSpecific(x,VEC_COOKIE); 2198 PetscCheckSameComm(mat,b); 2199 PetscCheckSameComm(mat,y); 2200 PetscCheckSameComm(mat,x); 2201 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,"x and b must be different vectors"); 2202 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Unfactored matrix"); 2203 if (mat->M != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->M,x->N); 2204 if (mat->N != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->N,b->N); 2205 if (mat->N != y->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec y: global dim %d %d",mat->N,y->N); 2206 if (x->n != y->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Vec x,Vec y: local dim %d %d",x->n,y->n); 2207 2208 ierr = PetscLogEventBegin(MAT_SolveTransposeAdd,mat,b,x,y);CHKERRQ(ierr); 2209 if (mat->ops->solvetransposeadd) { 2210 ierr = (*mat->ops->solvetransposeadd)(mat,b,y,x);CHKERRQ(ierr); 2211 } else { 2212 /* do the solve then the add manually */ 2213 if (x != y) { 2214 ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr); 2215 ierr = VecAXPY(&one,y,x);CHKERRQ(ierr); 2216 } else { 2217 ierr = VecDuplicate(x,&tmp);CHKERRQ(ierr); 2218 PetscLogObjectParent(mat,tmp); 2219 ierr = VecCopy(x,tmp);CHKERRQ(ierr); 2220 ierr = MatSolveTranspose(mat,b,x);CHKERRQ(ierr); 2221 ierr = VecAXPY(&one,tmp,x);CHKERRQ(ierr); 2222 ierr = VecDestroy(tmp);CHKERRQ(ierr); 2223 } 2224 } 2225 ierr = PetscLogEventEnd(MAT_SolveTransposeAdd,mat,b,x,y);CHKERRQ(ierr); 2226 PetscFunctionReturn(0); 2227 } 2228 /* ----------------------------------------------------------------*/ 2229 2230 #undef __FUNCT__ 2231 #define __FUNCT__ "MatRelax" 2232 /*@ 2233 MatRelax - Computes one relaxation sweep. 2234 2235 Collective on Mat and Vec 2236 2237 Input Parameters: 2238 + mat - the matrix 2239 . b - the right hand side 2240 . omega - the relaxation factor 2241 . flag - flag indicating the type of SOR (see below) 2242 . shift - diagonal shift 2243 - its - the number of iterations 2244 - lits - the number of local iterations 2245 2246 Output Parameters: 2247 . x - the solution (can contain an initial guess) 2248 2249 SOR Flags: 2250 . SOR_FORWARD_SWEEP - forward SOR 2251 . SOR_BACKWARD_SWEEP - backward SOR 2252 . SOR_SYMMETRIC_SWEEP - SSOR (symmetric SOR) 2253 . SOR_LOCAL_FORWARD_SWEEP - local forward SOR 2254 . SOR_LOCAL_BACKWARD_SWEEP - local forward SOR 2255 . SOR_LOCAL_SYMMETRIC_SWEEP - local SSOR 2256 . SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies 2257 upper/lower triangular part of matrix to 2258 vector (with omega) 2259 . SOR_ZERO_INITIAL_GUESS - zero initial guess 2260 2261 Notes: 2262 SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and 2263 SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings 2264 on each processor. 2265 2266 Application programmers will not generally use MatRelax() directly, 2267 but instead will employ the SLES/PC interface. 2268 2269 Notes for Advanced Users: 2270 The flags are implemented as bitwise inclusive or operations. 2271 For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP) 2272 to specify a zero initial guess for SSOR. 2273 2274 Most users should employ the simplified SLES interface for linear solvers 2275 instead of working directly with matrix algebra routines such as this. 2276 See, e.g., SLESCreate(). 2277 2278 Level: developer 2279 2280 Concepts: matrices^relaxation 2281 Concepts: matrices^SOR 2282 Concepts: matrices^Gauss-Seidel 2283 2284 @*/ 2285 int MatRelax(Mat mat,Vec b,PetscReal omega,MatSORType flag,PetscReal shift,int its,int lits,Vec x) 2286 { 2287 int ierr; 2288 2289 PetscFunctionBegin; 2290 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2291 PetscValidType(mat); 2292 MatPreallocated(mat); 2293 PetscValidHeaderSpecific(b,VEC_COOKIE); 2294 PetscValidHeaderSpecific(x,VEC_COOKIE); 2295 PetscCheckSameComm(mat,b); 2296 PetscCheckSameComm(mat,x); 2297 if (!mat->ops->relax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2298 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2299 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2300 if (mat->N != x->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec x: global dim %d %d",mat->N,x->N); 2301 if (mat->M != b->N) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: global dim %d %d",mat->M,b->N); 2302 if (mat->m != b->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"Mat mat,Vec b: local dim %d %d",mat->m,b->n); 2303 2304 ierr = PetscLogEventBegin(MAT_Relax,mat,b,x,0);CHKERRQ(ierr); 2305 ierr =(*mat->ops->relax)(mat,b,omega,flag,shift,its,lits,x);CHKERRQ(ierr); 2306 ierr = PetscLogEventEnd(MAT_Relax,mat,b,x,0);CHKERRQ(ierr); 2307 PetscFunctionReturn(0); 2308 } 2309 2310 #undef __FUNCT__ 2311 #define __FUNCT__ "MatCopy_Basic" 2312 /* 2313 Default matrix copy routine. 2314 */ 2315 int MatCopy_Basic(Mat A,Mat B,MatStructure str) 2316 { 2317 int ierr,i,rstart,rend,nz,*cwork; 2318 PetscScalar *vwork; 2319 2320 PetscFunctionBegin; 2321 ierr = MatZeroEntries(B);CHKERRQ(ierr); 2322 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 2323 for (i=rstart; i<rend; i++) { 2324 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2325 ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr); 2326 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2327 } 2328 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2329 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2330 PetscFunctionReturn(0); 2331 } 2332 2333 #undef __FUNCT__ 2334 #define __FUNCT__ "MatCopy" 2335 /*@C 2336 MatCopy - Copys a matrix to another matrix. 2337 2338 Collective on Mat 2339 2340 Input Parameters: 2341 + A - the matrix 2342 - str - SAME_NONZERO_PATTERN or DIFFERENT_NONZERO_PATTERN 2343 2344 Output Parameter: 2345 . B - where the copy is put 2346 2347 Notes: 2348 If you use SAME_NONZERO_PATTERN then the two matrices had better have the 2349 same nonzero pattern or the routine will crash. 2350 2351 MatCopy() copies the matrix entries of a matrix to another existing 2352 matrix (after first zeroing the second matrix). A related routine is 2353 MatConvert(), which first creates a new matrix and then copies the data. 2354 2355 Level: intermediate 2356 2357 Concepts: matrices^copying 2358 2359 .seealso: MatConvert() 2360 @*/ 2361 int MatCopy(Mat A,Mat B,MatStructure str) 2362 { 2363 int ierr; 2364 2365 PetscFunctionBegin; 2366 PetscValidHeaderSpecific(A,MAT_COOKIE); 2367 PetscValidHeaderSpecific(B,MAT_COOKIE); 2368 PetscValidType(A); 2369 MatPreallocated(A); 2370 PetscValidType(B); 2371 MatPreallocated(B); 2372 PetscCheckSameComm(A,B); 2373 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2374 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2375 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, 2376 A->N,B->N); 2377 2378 ierr = PetscLogEventBegin(MAT_Copy,A,B,0,0);CHKERRQ(ierr); 2379 if (A->ops->copy) { 2380 ierr = (*A->ops->copy)(A,B,str);CHKERRQ(ierr); 2381 } else { /* generic conversion */ 2382 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2383 } 2384 ierr = PetscLogEventEnd(MAT_Copy,A,B,0,0);CHKERRQ(ierr); 2385 PetscFunctionReturn(0); 2386 } 2387 2388 #include "petscsys.h" 2389 PetscTruth MatConvertRegisterAllCalled = PETSC_FALSE; 2390 PetscFList MatConvertList = 0; 2391 2392 #undef __FUNCT__ 2393 #define __FUNCT__ "MatConvertRegister" 2394 /*@C 2395 MatConvertRegister - Allows one to register a routine that reads matrices 2396 from a binary file for a particular matrix type. 2397 2398 Not Collective 2399 2400 Input Parameters: 2401 + type - the type of matrix (defined in include/petscmat.h), for example, MATSEQAIJ. 2402 - Converter - the function that reads the matrix from the binary file. 2403 2404 Level: developer 2405 2406 .seealso: MatConvertRegisterAll(), MatConvert() 2407 2408 @*/ 2409 int MatConvertRegister(char *sname,char *path,char *name,int (*function)(Mat,MatType,Mat*)) 2410 { 2411 int ierr; 2412 char fullname[256]; 2413 2414 PetscFunctionBegin; 2415 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 2416 ierr = PetscFListAdd(&MatConvertList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 2417 PetscFunctionReturn(0); 2418 } 2419 2420 #undef __FUNCT__ 2421 #define __FUNCT__ "MatConvert" 2422 /*@C 2423 MatConvert - Converts a matrix to another matrix, either of the same 2424 or different type. 2425 2426 Collective on Mat 2427 2428 Input Parameters: 2429 + mat - the matrix 2430 - newtype - new matrix type. Use MATSAME to create a new matrix of the 2431 same type as the original matrix. 2432 2433 Output Parameter: 2434 . M - pointer to place new matrix 2435 2436 Notes: 2437 MatConvert() first creates a new matrix and then copies the data from 2438 the first matrix. A related routine is MatCopy(), which copies the matrix 2439 entries of one matrix to another already existing matrix context. 2440 2441 Level: intermediate 2442 2443 Concepts: matrices^converting between storage formats 2444 2445 .seealso: MatCopy(), MatDuplicate() 2446 @*/ 2447 int MatConvert(Mat mat,MatType newtype,Mat *M) 2448 { 2449 int ierr; 2450 PetscTruth sametype,issame,flg; 2451 char convname[256],mtype[256]; 2452 2453 PetscFunctionBegin; 2454 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2455 PetscValidType(mat); 2456 MatPreallocated(mat); 2457 PetscValidPointer(M); 2458 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2459 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2460 2461 ierr = PetscOptionsGetString(PETSC_NULL,"-matconvert_type",mtype,256,&flg);CHKERRQ(ierr); 2462 if (flg) { 2463 newtype = mtype; 2464 } 2465 ierr = PetscLogEventBegin(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2466 2467 ierr = PetscTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(ierr); 2468 ierr = PetscStrcmp(newtype,"same",&issame);CHKERRQ(ierr); 2469 if ((sametype || issame) && mat->ops->duplicate) { 2470 ierr = (*mat->ops->duplicate)(mat,MAT_COPY_VALUES,M);CHKERRQ(ierr); 2471 } else { 2472 int (*conv)(Mat,MatType,Mat*); 2473 if (!MatConvertRegisterAllCalled) { 2474 ierr = MatConvertRegisterAll(PETSC_NULL);CHKERRQ(ierr); 2475 } 2476 ierr = PetscFListFind(mat->comm,MatConvertList,newtype,(void(**)(void))&conv);CHKERRQ(ierr); 2477 if (conv) { 2478 ierr = (*conv)(mat,newtype,M);CHKERRQ(ierr); 2479 } else { 2480 ierr = PetscStrcpy(convname,"MatConvert_");CHKERRQ(ierr); 2481 ierr = PetscStrcat(convname,mat->type_name);CHKERRQ(ierr); 2482 ierr = PetscStrcat(convname,"_");CHKERRQ(ierr); 2483 ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr); 2484 ierr = PetscStrcat(convname,"_C");CHKERRQ(ierr); 2485 ierr = PetscObjectQueryFunction((PetscObject)mat,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2486 if (conv) { 2487 ierr = (*conv)(mat,newtype,M);CHKERRQ(ierr); 2488 } else { 2489 if (mat->ops->convert) { 2490 ierr = (*mat->ops->convert)(mat,newtype,M);CHKERRQ(ierr); 2491 } else { 2492 ierr = MatConvert_Basic(mat,newtype,M);CHKERRQ(ierr); 2493 } 2494 } 2495 } 2496 } 2497 ierr = PetscLogEventEnd(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2498 PetscFunctionReturn(0); 2499 } 2500 2501 2502 #undef __FUNCT__ 2503 #define __FUNCT__ "MatDuplicate" 2504 /*@C 2505 MatDuplicate - Duplicates a matrix including the non-zero structure. 2506 2507 Collective on Mat 2508 2509 Input Parameters: 2510 + mat - the matrix 2511 - op - either MAT_DO_NOT_COPY_VALUES or MAT_COPY_VALUES, cause it to copy nonzero 2512 values as well or not 2513 2514 Output Parameter: 2515 . M - pointer to place new matrix 2516 2517 Level: intermediate 2518 2519 Concepts: matrices^duplicating 2520 2521 .seealso: MatCopy(), MatConvert() 2522 @*/ 2523 int MatDuplicate(Mat mat,MatDuplicateOption op,Mat *M) 2524 { 2525 int ierr; 2526 2527 PetscFunctionBegin; 2528 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2529 PetscValidType(mat); 2530 MatPreallocated(mat); 2531 PetscValidPointer(M); 2532 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2533 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2534 2535 *M = 0; 2536 ierr = PetscLogEventBegin(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2537 if (!mat->ops->duplicate) { 2538 SETERRQ(PETSC_ERR_SUP,"Not written for this matrix type"); 2539 } 2540 ierr = (*mat->ops->duplicate)(mat,op,M);CHKERRQ(ierr); 2541 ierr = PetscLogEventEnd(MAT_Convert,mat,0,0,0);CHKERRQ(ierr); 2542 PetscFunctionReturn(0); 2543 } 2544 2545 #undef __FUNCT__ 2546 #define __FUNCT__ "MatGetDiagonal" 2547 /*@ 2548 MatGetDiagonal - Gets the diagonal of a matrix. 2549 2550 Collective on Mat and Vec 2551 2552 Input Parameters: 2553 + mat - the matrix 2554 - v - the vector for storing the diagonal 2555 2556 Output Parameter: 2557 . v - the diagonal of the matrix 2558 2559 Notes: 2560 For the SeqAIJ matrix format, this routine may also be called 2561 on a LU factored matrix; in that case it routines the reciprocal of 2562 the diagonal entries in U. It returns the entries permuted by the 2563 row and column permutation used during the symbolic factorization. 2564 2565 Level: intermediate 2566 2567 Concepts: matrices^accessing diagonals 2568 2569 .seealso: MatGetRow(), MatGetSubmatrices(), MatGetSubmatrix(), MatGetRowMax() 2570 @*/ 2571 int MatGetDiagonal(Mat mat,Vec v) 2572 { 2573 int ierr; 2574 2575 PetscFunctionBegin; 2576 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2577 PetscValidType(mat); 2578 MatPreallocated(mat); 2579 PetscValidHeaderSpecific(v,VEC_COOKIE); 2580 /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */ 2581 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2582 if (!mat->ops->getdiagonal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2583 2584 ierr = (*mat->ops->getdiagonal)(mat,v);CHKERRQ(ierr); 2585 PetscFunctionReturn(0); 2586 } 2587 2588 #undef __FUNCT__ 2589 #define __FUNCT__ "MatGetRowMax" 2590 /*@ 2591 MatGetRowMax - Gets the maximum value (in absolute value) of each 2592 row of the matrix 2593 2594 Collective on Mat and Vec 2595 2596 Input Parameters: 2597 . mat - the matrix 2598 2599 Output Parameter: 2600 . v - the vector for storing the maximums 2601 2602 Level: intermediate 2603 2604 Concepts: matrices^getting row maximums 2605 2606 .seealso: MatGetDiagonal(), MatGetSubmatrices(), MatGetSubmatrix() 2607 @*/ 2608 int MatGetRowMax(Mat mat,Vec v) 2609 { 2610 int ierr; 2611 2612 PetscFunctionBegin; 2613 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2614 PetscValidType(mat); 2615 MatPreallocated(mat); 2616 PetscValidHeaderSpecific(v,VEC_COOKIE); 2617 /* PetscCheckSameComm(mat,v); Could be MPI vector but Seq matrix cause of two submatrix storage */ 2618 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2619 if (!mat->ops->getrowmax) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2620 2621 ierr = (*mat->ops->getrowmax)(mat,v);CHKERRQ(ierr); 2622 PetscFunctionReturn(0); 2623 } 2624 2625 #undef __FUNCT__ 2626 #define __FUNCT__ "MatTranspose" 2627 /*@C 2628 MatTranspose - Computes an in-place or out-of-place transpose of a matrix. 2629 2630 Collective on Mat 2631 2632 Input Parameter: 2633 . mat - the matrix to transpose 2634 2635 Output Parameters: 2636 . B - the transpose (or pass in PETSC_NULL for an in-place transpose) 2637 2638 Level: intermediate 2639 2640 Concepts: matrices^transposing 2641 2642 .seealso: MatMultTranspose(), MatMultTransposeAdd() 2643 @*/ 2644 int MatTranspose(Mat mat,Mat *B) 2645 { 2646 int ierr; 2647 2648 PetscFunctionBegin; 2649 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2650 PetscValidType(mat); 2651 MatPreallocated(mat); 2652 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2653 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2654 if (!mat->ops->transpose) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2655 ierr = (*mat->ops->transpose)(mat,B);CHKERRQ(ierr); 2656 PetscFunctionReturn(0); 2657 } 2658 2659 #undef __FUNCT__ 2660 #define __FUNCT__ "MatPermute" 2661 /*@C 2662 MatPermute - Creates a new matrix with rows and columns permuted from the 2663 original. 2664 2665 Collective on Mat 2666 2667 Input Parameters: 2668 + mat - the matrix to permute 2669 . row - row permutation, each processor supplies only the permutation for its rows 2670 - col - column permutation, each processor needs the entire column permutation, that is 2671 this is the same size as the total number of columns in the matrix 2672 2673 Output Parameters: 2674 . B - the permuted matrix 2675 2676 Level: advanced 2677 2678 Concepts: matrices^permuting 2679 2680 .seealso: MatGetOrdering() 2681 @*/ 2682 int MatPermute(Mat mat,IS row,IS col,Mat *B) 2683 { 2684 int ierr; 2685 2686 PetscFunctionBegin; 2687 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2688 PetscValidType(mat); 2689 MatPreallocated(mat); 2690 PetscValidHeaderSpecific(row,IS_COOKIE); 2691 PetscValidHeaderSpecific(col,IS_COOKIE); 2692 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2693 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2694 if (!mat->ops->permute) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2695 ierr = (*mat->ops->permute)(mat,row,col,B);CHKERRQ(ierr); 2696 PetscFunctionReturn(0); 2697 } 2698 2699 #undef __FUNCT__ 2700 #define __FUNCT__ "MatPermuteSparsify" 2701 /*@C 2702 MatPermuteSparsify - Creates a new matrix with rows and columns permuted from the 2703 original and sparsified to the prescribed tolerance. 2704 2705 Collective on Mat 2706 2707 Input Parameters: 2708 + A - The matrix to permute 2709 . band - The half-bandwidth of the sparsified matrix, or PETSC_DECIDE 2710 . frac - The half-bandwidth as a fraction of the total size, or 0.0 2711 . tol - The drop tolerance 2712 . rowp - The row permutation 2713 - colp - The column permutation 2714 2715 Output Parameter: 2716 . B - The permuted, sparsified matrix 2717 2718 Level: advanced 2719 2720 Note: 2721 The default behavior (band = PETSC_DECIDE and frac = 0.0) is to 2722 restrict the half-bandwidth of the resulting matrix to 5% of the 2723 total matrix size. 2724 2725 .keywords: matrix, permute, sparsify 2726 2727 .seealso: MatGetOrdering(), MatPermute() 2728 @*/ 2729 int MatPermuteSparsify(Mat A, int band, double frac, double tol, IS rowp, IS colp, Mat *B) 2730 { 2731 IS irowp, icolp; 2732 int *rows, *cols; 2733 int M, N, locRowStart, locRowEnd; 2734 int nz, newNz; 2735 int *cwork, *cnew; 2736 PetscScalar *vwork, *vnew; 2737 int bw, size; 2738 int row, locRow, newRow, col, newCol; 2739 int ierr; 2740 2741 PetscFunctionBegin; 2742 PetscValidHeaderSpecific(A, MAT_COOKIE); 2743 PetscValidHeaderSpecific(rowp, IS_COOKIE); 2744 PetscValidHeaderSpecific(colp, IS_COOKIE); 2745 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix"); 2746 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 2747 if (!A->ops->permutesparsify) { 2748 ierr = MatGetSize(A, &M, &N); CHKERRQ(ierr); 2749 ierr = MatGetOwnershipRange(A, &locRowStart, &locRowEnd); CHKERRQ(ierr); 2750 ierr = ISGetSize(rowp, &size); CHKERRQ(ierr); 2751 if (size != M) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for row permutation, should be %d", size, M); 2752 ierr = ISGetSize(colp, &size); CHKERRQ(ierr); 2753 if (size != N) SETERRQ2(PETSC_ERR_ARG_WRONG, "Wrong size %d for column permutation, should be %d", size, N); 2754 ierr = ISInvertPermutation(rowp, 0, &irowp); CHKERRQ(ierr); 2755 ierr = ISGetIndices(irowp, &rows); CHKERRQ(ierr); 2756 ierr = ISInvertPermutation(colp, 0, &icolp); CHKERRQ(ierr); 2757 ierr = ISGetIndices(icolp, &cols); CHKERRQ(ierr); 2758 ierr = PetscMalloc(N * sizeof(int), &cnew); CHKERRQ(ierr); 2759 ierr = PetscMalloc(N * sizeof(PetscScalar), &vnew); CHKERRQ(ierr); 2760 2761 /* Setup bandwidth to include */ 2762 if (band == PETSC_DECIDE) { 2763 if (frac <= 0.0) 2764 bw = (int) (M * 0.05); 2765 else 2766 bw = (int) (M * frac); 2767 } else { 2768 if (band <= 0) SETERRQ(PETSC_ERR_ARG_WRONG, "Bandwidth must be a positive integer"); 2769 bw = band; 2770 } 2771 2772 /* Put values into new matrix */ 2773 ierr = MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, B); CHKERRQ(ierr); 2774 for(row = locRowStart, locRow = 0; row < locRowEnd; row++, locRow++) { 2775 ierr = MatGetRow(A, row, &nz, &cwork, &vwork); CHKERRQ(ierr); 2776 newRow = rows[locRow]+locRowStart; 2777 for(col = 0, newNz = 0; col < nz; col++) { 2778 newCol = cols[cwork[col]]; 2779 if ((newCol >= newRow - bw) && (newCol < newRow + bw) && (PetscAbsScalar(vwork[col]) >= tol)) { 2780 cnew[newNz] = newCol; 2781 vnew[newNz] = vwork[col]; 2782 newNz++; 2783 } 2784 } 2785 ierr = MatSetValues(*B, 1, &newRow, newNz, cnew, vnew, INSERT_VALUES); CHKERRQ(ierr); 2786 ierr = MatRestoreRow(A, row, &nz, &cwork, &vwork); CHKERRQ(ierr); 2787 } 2788 ierr = PetscFree(cnew); CHKERRQ(ierr); 2789 ierr = PetscFree(vnew); CHKERRQ(ierr); 2790 ierr = MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2791 ierr = MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2792 ierr = ISRestoreIndices(irowp, &rows); CHKERRQ(ierr); 2793 ierr = ISRestoreIndices(icolp, &cols); CHKERRQ(ierr); 2794 ierr = ISDestroy(irowp); CHKERRQ(ierr); 2795 ierr = ISDestroy(icolp); CHKERRQ(ierr); 2796 } else { 2797 ierr = (*A->ops->permutesparsify)(A, band, frac, tol, rowp, colp, B); CHKERRQ(ierr); 2798 } 2799 PetscFunctionReturn(0); 2800 } 2801 2802 #undef __FUNCT__ 2803 #define __FUNCT__ "MatEqual" 2804 /*@ 2805 MatEqual - Compares two matrices. 2806 2807 Collective on Mat 2808 2809 Input Parameters: 2810 + A - the first matrix 2811 - B - the second matrix 2812 2813 Output Parameter: 2814 . flg - PETSC_TRUE if the matrices are equal; PETSC_FALSE otherwise. 2815 2816 Level: intermediate 2817 2818 Concepts: matrices^equality between 2819 @*/ 2820 int MatEqual(Mat A,Mat B,PetscTruth *flg) 2821 { 2822 int ierr; 2823 2824 PetscFunctionBegin; 2825 PetscValidHeaderSpecific(A,MAT_COOKIE); 2826 PetscValidHeaderSpecific(B,MAT_COOKIE); 2827 PetscValidType(A); 2828 MatPreallocated(A); 2829 PetscValidType(B); 2830 MatPreallocated(B); 2831 PetscValidIntPointer(flg); 2832 PetscCheckSameComm(A,B); 2833 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2834 if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2835 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); 2836 if (!A->ops->equal) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",A->type_name); 2837 ierr = (*A->ops->equal)(A,B,flg);CHKERRQ(ierr); 2838 PetscFunctionReturn(0); 2839 } 2840 2841 #undef __FUNCT__ 2842 #define __FUNCT__ "MatDiagonalScale" 2843 /*@ 2844 MatDiagonalScale - Scales a matrix on the left and right by diagonal 2845 matrices that are stored as vectors. Either of the two scaling 2846 matrices can be PETSC_NULL. 2847 2848 Collective on Mat 2849 2850 Input Parameters: 2851 + mat - the matrix to be scaled 2852 . l - the left scaling vector (or PETSC_NULL) 2853 - r - the right scaling vector (or PETSC_NULL) 2854 2855 Notes: 2856 MatDiagonalScale() computes A = LAR, where 2857 L = a diagonal matrix, R = a diagonal matrix 2858 2859 Level: intermediate 2860 2861 Concepts: matrices^diagonal scaling 2862 Concepts: diagonal scaling of matrices 2863 2864 .seealso: MatScale() 2865 @*/ 2866 int MatDiagonalScale(Mat mat,Vec l,Vec r) 2867 { 2868 int ierr; 2869 2870 PetscFunctionBegin; 2871 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2872 PetscValidType(mat); 2873 MatPreallocated(mat); 2874 if (!mat->ops->diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2875 if (l) {PetscValidHeaderSpecific(l,VEC_COOKIE);PetscCheckSameComm(mat,l);} 2876 if (r) {PetscValidHeaderSpecific(r,VEC_COOKIE);PetscCheckSameComm(mat,r);} 2877 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2878 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2879 2880 ierr = PetscLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 2881 ierr = (*mat->ops->diagonalscale)(mat,l,r);CHKERRQ(ierr); 2882 ierr = PetscLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 2883 PetscFunctionReturn(0); 2884 } 2885 2886 #undef __FUNCT__ 2887 #define __FUNCT__ "MatScale" 2888 /*@ 2889 MatScale - Scales all elements of a matrix by a given number. 2890 2891 Collective on Mat 2892 2893 Input Parameters: 2894 + mat - the matrix to be scaled 2895 - a - the scaling value 2896 2897 Output Parameter: 2898 . mat - the scaled matrix 2899 2900 Level: intermediate 2901 2902 Concepts: matrices^scaling all entries 2903 2904 .seealso: MatDiagonalScale() 2905 @*/ 2906 int MatScale(PetscScalar *a,Mat mat) 2907 { 2908 int ierr; 2909 2910 PetscFunctionBegin; 2911 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2912 PetscValidType(mat); 2913 MatPreallocated(mat); 2914 PetscValidScalarPointer(a); 2915 if (!mat->ops->scale) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2916 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2917 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2918 2919 ierr = PetscLogEventBegin(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 2920 ierr = (*mat->ops->scale)(a,mat);CHKERRQ(ierr); 2921 ierr = PetscLogEventEnd(MAT_Scale,mat,0,0,0);CHKERRQ(ierr); 2922 PetscFunctionReturn(0); 2923 } 2924 2925 #undef __FUNCT__ 2926 #define __FUNCT__ "MatNorm" 2927 /*@ 2928 MatNorm - Calculates various norms of a matrix. 2929 2930 Collective on Mat 2931 2932 Input Parameters: 2933 + mat - the matrix 2934 - type - the type of norm, NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY 2935 2936 Output Parameters: 2937 . nrm - the resulting norm 2938 2939 Level: intermediate 2940 2941 Concepts: matrices^norm 2942 Concepts: norm^of matrix 2943 @*/ 2944 int MatNorm(Mat mat,NormType type,PetscReal *nrm) 2945 { 2946 int ierr; 2947 2948 PetscFunctionBegin; 2949 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2950 PetscValidType(mat); 2951 MatPreallocated(mat); 2952 PetscValidScalarPointer(nrm); 2953 2954 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2955 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2956 if (!mat->ops->norm) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 2957 ierr = (*mat->ops->norm)(mat,type,nrm);CHKERRQ(ierr); 2958 PetscFunctionReturn(0); 2959 } 2960 2961 /* 2962 This variable is used to prevent counting of MatAssemblyBegin() that 2963 are called from within a MatAssemblyEnd(). 2964 */ 2965 static int MatAssemblyEnd_InUse = 0; 2966 #undef __FUNCT__ 2967 #define __FUNCT__ "MatAssemblyBegin" 2968 /*@ 2969 MatAssemblyBegin - Begins assembling the matrix. This routine should 2970 be called after completing all calls to MatSetValues(). 2971 2972 Collective on Mat 2973 2974 Input Parameters: 2975 + mat - the matrix 2976 - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 2977 2978 Notes: 2979 MatSetValues() generally caches the values. The matrix is ready to 2980 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 2981 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 2982 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 2983 using the matrix. 2984 2985 Level: beginner 2986 2987 Concepts: matrices^assembling 2988 2989 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssembled() 2990 @*/ 2991 int MatAssemblyBegin(Mat mat,MatAssemblyType type) 2992 { 2993 int ierr; 2994 2995 PetscFunctionBegin; 2996 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2997 PetscValidType(mat); 2998 MatPreallocated(mat); 2999 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix.\nDid you forget to call MatSetUnfactored()?"); 3000 if (mat->assembled) { 3001 mat->was_assembled = PETSC_TRUE; 3002 mat->assembled = PETSC_FALSE; 3003 } 3004 if (!MatAssemblyEnd_InUse) { 3005 ierr = PetscLogEventBegin(MAT_AssemblyBegin,mat,0,0,0);CHKERRQ(ierr); 3006 if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);} 3007 ierr = PetscLogEventEnd(MAT_AssemblyBegin,mat,0,0,0);CHKERRQ(ierr); 3008 } else { 3009 if (mat->ops->assemblybegin){ierr = (*mat->ops->assemblybegin)(mat,type);CHKERRQ(ierr);} 3010 } 3011 PetscFunctionReturn(0); 3012 } 3013 3014 #undef __FUNCT__ 3015 #define __FUNCT__ "MatAssembed" 3016 /*@ 3017 MatAssembled - Indicates if a matrix has been assembled and is ready for 3018 use; for example, in matrix-vector product. 3019 3020 Collective on Mat 3021 3022 Input Parameter: 3023 . mat - the matrix 3024 3025 Output Parameter: 3026 . assembled - PETSC_TRUE or PETSC_FALSE 3027 3028 Level: advanced 3029 3030 Concepts: matrices^assembled? 3031 3032 .seealso: MatAssemblyEnd(), MatSetValues(), MatAssemblyBegin() 3033 @*/ 3034 int MatAssembled(Mat mat,PetscTruth *assembled) 3035 { 3036 PetscFunctionBegin; 3037 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3038 PetscValidType(mat); 3039 MatPreallocated(mat); 3040 *assembled = mat->assembled; 3041 PetscFunctionReturn(0); 3042 } 3043 3044 #undef __FUNCT__ 3045 #define __FUNCT__ "MatView_Private" 3046 /* 3047 Processes command line options to determine if/how a matrix 3048 is to be viewed. Called by MatAssemblyEnd() and MatLoad(). 3049 */ 3050 int MatView_Private(Mat mat) 3051 { 3052 int ierr; 3053 PetscTruth flg; 3054 3055 PetscFunctionBegin; 3056 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_info",&flg);CHKERRQ(ierr); 3057 if (flg) { 3058 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 3059 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3060 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3061 } 3062 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_info_detailed",&flg);CHKERRQ(ierr); 3063 if (flg) { 3064 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_INFO_LONG);CHKERRQ(ierr); 3065 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3066 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3067 } 3068 ierr = PetscOptionsHasName(mat->prefix,"-mat_view",&flg);CHKERRQ(ierr); 3069 if (flg) { 3070 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3071 } 3072 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_matlab",&flg);CHKERRQ(ierr); 3073 if (flg) { 3074 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(mat->comm),PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); 3075 ierr = MatView(mat,PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3076 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 3077 } 3078 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_draw",&flg);CHKERRQ(ierr); 3079 if (flg) { 3080 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_contour",&flg);CHKERRQ(ierr); 3081 if (flg) { 3082 PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 3083 } 3084 ierr = MatView(mat,PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3085 ierr = PetscViewerFlush(PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3086 if (flg) { 3087 PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));CHKERRQ(ierr); 3088 } 3089 } 3090 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_socket",&flg);CHKERRQ(ierr); 3091 if (flg) { 3092 ierr = MatView(mat,PETSC_VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr); 3093 ierr = PetscViewerFlush(PETSC_VIEWER_SOCKET_(mat->comm));CHKERRQ(ierr); 3094 } 3095 ierr = PetscOptionsHasName(mat->prefix,"-mat_view_binary",&flg);CHKERRQ(ierr); 3096 if (flg) { 3097 ierr = MatView(mat,PETSC_VIEWER_BINARY_(mat->comm));CHKERRQ(ierr); 3098 ierr = PetscViewerFlush(PETSC_VIEWER_BINARY_(mat->comm));CHKERRQ(ierr); 3099 } 3100 PetscFunctionReturn(0); 3101 } 3102 3103 #undef __FUNCT__ 3104 #define __FUNCT__ "MatAssemblyEnd" 3105 /*@ 3106 MatAssemblyEnd - Completes assembling the matrix. This routine should 3107 be called after MatAssemblyBegin(). 3108 3109 Collective on Mat 3110 3111 Input Parameters: 3112 + mat - the matrix 3113 - type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 3114 3115 Options Database Keys: 3116 + -mat_view_info - Prints info on matrix at conclusion of MatEndAssembly() 3117 . -mat_view_info_detailed - Prints more detailed info 3118 . -mat_view - Prints matrix in ASCII format 3119 . -mat_view_matlab - Prints matrix in Matlab format 3120 . -mat_view_draw - PetscDraws nonzero structure of matrix, using MatView() and PetscDrawOpenX(). 3121 . -display <name> - Sets display name (default is host) 3122 - -draw_pause <sec> - Sets number of seconds to pause after display 3123 3124 Notes: 3125 MatSetValues() generally caches the values. The matrix is ready to 3126 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 3127 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 3128 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 3129 using the matrix. 3130 3131 Level: beginner 3132 3133 .seealso: MatAssemblyBegin(), MatSetValues(), PetscDrawOpenX(), MatView(), MatAssembled() 3134 @*/ 3135 int MatAssemblyEnd(Mat mat,MatAssemblyType type) 3136 { 3137 int ierr; 3138 static int inassm = 0; 3139 3140 PetscFunctionBegin; 3141 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3142 PetscValidType(mat); 3143 MatPreallocated(mat); 3144 3145 inassm++; 3146 MatAssemblyEnd_InUse++; 3147 if (MatAssemblyEnd_InUse == 1) { /* Do the logging only the first time through */ 3148 ierr = PetscLogEventBegin(MAT_AssemblyEnd,mat,0,0,0);CHKERRQ(ierr); 3149 if (mat->ops->assemblyend) { 3150 ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr); 3151 } 3152 ierr = PetscLogEventEnd(MAT_AssemblyEnd,mat,0,0,0);CHKERRQ(ierr); 3153 } else { 3154 if (mat->ops->assemblyend) { 3155 ierr = (*mat->ops->assemblyend)(mat,type);CHKERRQ(ierr); 3156 } 3157 } 3158 3159 /* Flush assembly is not a true assembly */ 3160 if (type != MAT_FLUSH_ASSEMBLY) { 3161 mat->assembled = PETSC_TRUE; mat->num_ass++; 3162 } 3163 mat->insertmode = NOT_SET_VALUES; 3164 MatAssemblyEnd_InUse--; 3165 3166 if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) { 3167 ierr = MatView_Private(mat);CHKERRQ(ierr); 3168 } 3169 inassm--; 3170 PetscFunctionReturn(0); 3171 } 3172 3173 3174 #undef __FUNCT__ 3175 #define __FUNCT__ "MatCompress" 3176 /*@ 3177 MatCompress - Tries to store the matrix in as little space as 3178 possible. May fail if memory is already fully used, since it 3179 tries to allocate new space. 3180 3181 Collective on Mat 3182 3183 Input Parameters: 3184 . mat - the matrix 3185 3186 Level: advanced 3187 3188 @*/ 3189 int MatCompress(Mat mat) 3190 { 3191 int ierr; 3192 3193 PetscFunctionBegin; 3194 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3195 PetscValidType(mat); 3196 MatPreallocated(mat); 3197 if (mat->ops->compress) {ierr = (*mat->ops->compress)(mat);CHKERRQ(ierr);} 3198 PetscFunctionReturn(0); 3199 } 3200 3201 #undef __FUNCT__ 3202 #define __FUNCT__ "MatSetOption" 3203 /*@ 3204 MatSetOption - Sets a parameter option for a matrix. Some options 3205 may be specific to certain storage formats. Some options 3206 determine how values will be inserted (or added). Sorted, 3207 row-oriented input will generally assemble the fastest. The default 3208 is row-oriented, nonsorted input. 3209 3210 Collective on Mat 3211 3212 Input Parameters: 3213 + mat - the matrix 3214 - option - the option, one of those listed below (and possibly others), 3215 e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERR 3216 3217 Options Describing Matrix Structure: 3218 + MAT_SYMMETRIC - symmetric in terms of both structure and value 3219 - MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure 3220 3221 Options For Use with MatSetValues(): 3222 Insert a logically dense subblock, which can be 3223 + MAT_ROW_ORIENTED - row-oriented 3224 . MAT_COLUMN_ORIENTED - column-oriented 3225 . MAT_ROWS_SORTED - sorted by row 3226 . MAT_ROWS_UNSORTED - not sorted by row 3227 . MAT_COLUMNS_SORTED - sorted by column 3228 - MAT_COLUMNS_UNSORTED - not sorted by column 3229 3230 Not these options reflect the data you pass in with MatSetValues(); it has 3231 nothing to do with how the data is stored internally in the matrix 3232 data structure. 3233 3234 When (re)assembling a matrix, we can restrict the input for 3235 efficiency/debugging purposes. These options include 3236 + MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be 3237 allowed if they generate a new nonzero 3238 . MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed 3239 . MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if 3240 they generate a nonzero in a new diagonal (for block diagonal format only) 3241 . MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only) 3242 . MAT_IGNORE_OFF_PROC_ENTRIES - drops off-processor entries 3243 . MAT_NEW_NONZERO_LOCATION_ERR - generates an error for new matrix entry 3244 - MAT_USE_HASH_TABLE - uses a hash table to speed up matrix assembly 3245 3246 Notes: 3247 Some options are relevant only for particular matrix types and 3248 are thus ignored by others. Other options are not supported by 3249 certain matrix types and will generate an error message if set. 3250 3251 If using a Fortran 77 module to compute a matrix, one may need to 3252 use the column-oriented option (or convert to the row-oriented 3253 format). 3254 3255 MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion 3256 that would generate a new entry in the nonzero structure is instead 3257 ignored. Thus, if memory has not alredy been allocated for this particular 3258 data, then the insertion is ignored. For dense matrices, in which 3259 the entire array is allocated, no entries are ever ignored. 3260 3261 MAT_NEW_NONZERO_LOCATION_ERR indicates that any add or insertion 3262 that would generate a new entry in the nonzero structure instead produces 3263 an error. (Currently supported for AIJ and BAIJ formats only.) 3264 This is a useful flag when using SAME_NONZERO_PATTERN in calling 3265 SLESSetOperators() to ensure that the nonzero pattern truely does 3266 remain unchanged. 3267 3268 MAT_NEW_NONZERO_ALLOCATION_ERR indicates that any add or insertion 3269 that would generate a new entry that has not been preallocated will 3270 instead produce an error. (Currently supported for AIJ and BAIJ formats 3271 only.) This is a useful flag when debugging matrix memory preallocation. 3272 3273 MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for 3274 other processors should be dropped, rather than stashed. 3275 This is useful if you know that the "owning" processor is also 3276 always generating the correct matrix entries, so that PETSc need 3277 not transfer duplicate entries generated on another processor. 3278 3279 MAT_USE_HASH_TABLE indicates that a hash table be used to improve the 3280 searches during matrix assembly. When this flag is set, the hash table 3281 is created during the first Matrix Assembly. This hash table is 3282 used the next time through, during MatSetVaules()/MatSetVaulesBlocked() 3283 to improve the searching of indices. MAT_NO_NEW_NONZERO_LOCATIONS flag 3284 should be used with MAT_USE_HASH_TABLE flag. This option is currently 3285 supported by MATMPIBAIJ format only. 3286 3287 MAT_KEEP_ZEROED_ROWS indicates when MatZeroRows() is called the zeroed entries 3288 are kept in the nonzero structure 3289 3290 MAT_IGNORE_ZERO_ENTRIES - when using ADD_VALUES for AIJ matrices this will stop 3291 zero values from creating a zero location in the matrix 3292 3293 MAT_USE_INODES - indicates using inode version of the code - works with AIJ and 3294 ROWBS matrix types 3295 3296 MAT_DO_NOT_USE_INODES - indicates not using inode version of the code - works 3297 with AIJ and ROWBS matrix types 3298 3299 Level: intermediate 3300 3301 Concepts: matrices^setting options 3302 3303 @*/ 3304 int MatSetOption(Mat mat,MatOption op) 3305 { 3306 int ierr; 3307 3308 PetscFunctionBegin; 3309 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3310 PetscValidType(mat); 3311 MatPreallocated(mat); 3312 switch (op) { 3313 case MAT_SYMMETRIC: 3314 mat->symmetric = PETSC_TRUE; 3315 mat->structurally_symmetric = PETSC_TRUE; 3316 break; 3317 case MAT_STRUCTURALLY_SYMMETRIC: 3318 mat->structurally_symmetric = PETSC_TRUE; 3319 break; 3320 default: 3321 if (mat->ops->setoption) {ierr = (*mat->ops->setoption)(mat,op);CHKERRQ(ierr);} 3322 break; 3323 } 3324 PetscFunctionReturn(0); 3325 } 3326 3327 #undef __FUNCT__ 3328 #define __FUNCT__ "MatZeroEntries" 3329 /*@ 3330 MatZeroEntries - Zeros all entries of a matrix. For sparse matrices 3331 this routine retains the old nonzero structure. 3332 3333 Collective on Mat 3334 3335 Input Parameters: 3336 . mat - the matrix 3337 3338 Level: intermediate 3339 3340 Concepts: matrices^zeroing 3341 3342 .seealso: MatZeroRows() 3343 @*/ 3344 int MatZeroEntries(Mat mat) 3345 { 3346 int ierr; 3347 3348 PetscFunctionBegin; 3349 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3350 PetscValidType(mat); 3351 MatPreallocated(mat); 3352 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3353 if (!mat->ops->zeroentries) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3354 3355 ierr = PetscLogEventBegin(MAT_ZeroEntries,mat,0,0,0);CHKERRQ(ierr); 3356 ierr = (*mat->ops->zeroentries)(mat);CHKERRQ(ierr); 3357 ierr = PetscLogEventEnd(MAT_ZeroEntries,mat,0,0,0);CHKERRQ(ierr); 3358 PetscFunctionReturn(0); 3359 } 3360 3361 #undef __FUNCT__ 3362 #define __FUNCT__ "MatZeroRows" 3363 /*@C 3364 MatZeroRows - Zeros all entries (except possibly the main diagonal) 3365 of a set of rows of a matrix. 3366 3367 Collective on Mat 3368 3369 Input Parameters: 3370 + mat - the matrix 3371 . is - index set of rows to remove 3372 - diag - pointer to value put in all diagonals of eliminated rows. 3373 Note that diag is not a pointer to an array, but merely a 3374 pointer to a single value. 3375 3376 Notes: 3377 For the AIJ and BAIJ matrix formats this removes the old nonzero structure, 3378 but does not release memory. For the dense and block diagonal 3379 formats this does not alter the nonzero structure. 3380 3381 If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure 3382 of the matrix is not changed (even for AIJ and BAIJ matrices) the values are 3383 merely zeroed. 3384 3385 The user can set a value in the diagonal entry (or for the AIJ and 3386 row formats can optionally remove the main diagonal entry from the 3387 nonzero structure as well, by passing a null pointer (PETSC_NULL 3388 in C or PETSC_NULL_SCALAR in Fortran) as the final argument). 3389 3390 For the parallel case, all processes that share the matrix (i.e., 3391 those in the communicator used for matrix creation) MUST call this 3392 routine, regardless of whether any rows being zeroed are owned by 3393 them. 3394 3395 3396 Level: intermediate 3397 3398 Concepts: matrices^zeroing rows 3399 3400 .seealso: MatZeroEntries(), MatZeroRowsLocal(), MatSetOption() 3401 @*/ 3402 int MatZeroRows(Mat mat,IS is,PetscScalar *diag) 3403 { 3404 int ierr; 3405 3406 PetscFunctionBegin; 3407 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3408 PetscValidType(mat); 3409 MatPreallocated(mat); 3410 PetscValidHeaderSpecific(is,IS_COOKIE); 3411 if (diag) PetscValidScalarPointer(diag); 3412 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3413 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3414 if (!mat->ops->zerorows) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3415 3416 ierr = (*mat->ops->zerorows)(mat,is,diag);CHKERRQ(ierr); 3417 ierr = MatView_Private(mat);CHKERRQ(ierr); 3418 PetscFunctionReturn(0); 3419 } 3420 3421 #undef __FUNCT__ 3422 #define __FUNCT__ "MatZeroRowsLocal" 3423 /*@C 3424 MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal) 3425 of a set of rows of a matrix; using local numbering of rows. 3426 3427 Collective on Mat 3428 3429 Input Parameters: 3430 + mat - the matrix 3431 . is - index set of rows to remove 3432 - diag - pointer to value put in all diagonals of eliminated rows. 3433 Note that diag is not a pointer to an array, but merely a 3434 pointer to a single value. 3435 3436 Notes: 3437 Before calling MatZeroRowsLocal(), the user must first set the 3438 local-to-global mapping by calling MatSetLocalToGlobalMapping(). 3439 3440 For the AIJ matrix formats this removes the old nonzero structure, 3441 but does not release memory. For the dense and block diagonal 3442 formats this does not alter the nonzero structure. 3443 3444 If the option MatSetOption(mat,MAT_KEEP_ZEROED_ROWS) the nonzero structure 3445 of the matrix is not changed (even for AIJ and BAIJ matrices) the values are 3446 merely zeroed. 3447 3448 The user can set a value in the diagonal entry (or for the AIJ and 3449 row formats can optionally remove the main diagonal entry from the 3450 nonzero structure as well, by passing a null pointer (PETSC_NULL 3451 in C or PETSC_NULL_SCALAR in Fortran) as the final argument). 3452 3453 Level: intermediate 3454 3455 Concepts: matrices^zeroing 3456 3457 .seealso: MatZeroEntries(), MatZeroRows(), MatSetLocalToGlobalMapping 3458 @*/ 3459 int MatZeroRowsLocal(Mat mat,IS is,PetscScalar *diag) 3460 { 3461 int ierr; 3462 IS newis; 3463 3464 PetscFunctionBegin; 3465 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3466 PetscValidType(mat); 3467 MatPreallocated(mat); 3468 PetscValidHeaderSpecific(is,IS_COOKIE); 3469 if (diag) PetscValidScalarPointer(diag); 3470 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3471 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3472 3473 if (mat->ops->zerorowslocal) { 3474 ierr = (*mat->ops->zerorowslocal)(mat,is,diag);CHKERRQ(ierr); 3475 } else { 3476 if (!mat->mapping) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Need to provide local to global mapping to matrix first"); 3477 ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis);CHKERRQ(ierr); 3478 ierr = (*mat->ops->zerorows)(mat,newis,diag);CHKERRQ(ierr); 3479 ierr = ISDestroy(newis);CHKERRQ(ierr); 3480 } 3481 PetscFunctionReturn(0); 3482 } 3483 3484 #undef __FUNCT__ 3485 #define __FUNCT__ "MatGetSize" 3486 /*@ 3487 MatGetSize - Returns the numbers of rows and columns in a matrix. 3488 3489 Not Collective 3490 3491 Input Parameter: 3492 . mat - the matrix 3493 3494 Output Parameters: 3495 + m - the number of global rows 3496 - n - the number of global columns 3497 3498 Level: beginner 3499 3500 Concepts: matrices^size 3501 3502 .seealso: MatGetLocalSize() 3503 @*/ 3504 int MatGetSize(Mat mat,int *m,int* n) 3505 { 3506 PetscFunctionBegin; 3507 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3508 if (m) *m = mat->M; 3509 if (n) *n = mat->N; 3510 PetscFunctionReturn(0); 3511 } 3512 3513 #undef __FUNCT__ 3514 #define __FUNCT__ "MatGetLocalSize" 3515 /*@ 3516 MatGetLocalSize - Returns the number of rows and columns in a matrix 3517 stored locally. This information may be implementation dependent, so 3518 use with care. 3519 3520 Not Collective 3521 3522 Input Parameters: 3523 . mat - the matrix 3524 3525 Output Parameters: 3526 + m - the number of local rows 3527 - n - the number of local columns 3528 3529 Level: beginner 3530 3531 Concepts: matrices^local size 3532 3533 .seealso: MatGetSize() 3534 @*/ 3535 int MatGetLocalSize(Mat mat,int *m,int* n) 3536 { 3537 PetscFunctionBegin; 3538 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3539 if (m) *m = mat->m; 3540 if (n) *n = mat->n; 3541 PetscFunctionReturn(0); 3542 } 3543 3544 #undef __FUNCT__ 3545 #define __FUNCT__ "MatGetOwnershipRange" 3546 /*@ 3547 MatGetOwnershipRange - Returns the range of matrix rows owned by 3548 this processor, assuming that the matrix is laid out with the first 3549 n1 rows on the first processor, the next n2 rows on the second, etc. 3550 For certain parallel layouts this range may not be well defined. 3551 3552 Not Collective 3553 3554 Input Parameters: 3555 . mat - the matrix 3556 3557 Output Parameters: 3558 + m - the global index of the first local row 3559 - n - one more than the global index of the last local row 3560 3561 Level: beginner 3562 3563 Concepts: matrices^row ownership 3564 @*/ 3565 int MatGetOwnershipRange(Mat mat,int *m,int* n) 3566 { 3567 int ierr; 3568 3569 PetscFunctionBegin; 3570 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3571 PetscValidType(mat); 3572 MatPreallocated(mat); 3573 if (m) PetscValidIntPointer(m); 3574 if (n) PetscValidIntPointer(n); 3575 ierr = PetscMapGetLocalRange(mat->rmap,m,n);CHKERRQ(ierr); 3576 PetscFunctionReturn(0); 3577 } 3578 3579 #undef __FUNCT__ 3580 #define __FUNCT__ "MatILUFactorSymbolic" 3581 /*@ 3582 MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix. 3583 Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric() 3584 to complete the factorization. 3585 3586 Collective on Mat 3587 3588 Input Parameters: 3589 + mat - the matrix 3590 . row - row permutation 3591 . column - column permutation 3592 - info - structure containing 3593 $ levels - number of levels of fill. 3594 $ expected fill - as ratio of original fill. 3595 $ 1 or 0 - indicating force fill on diagonal (improves robustness for matrices 3596 missing diagonal entries) 3597 3598 Output Parameters: 3599 . fact - new matrix that has been symbolically factored 3600 3601 Notes: 3602 See the users manual for additional information about 3603 choosing the fill factor for better efficiency. 3604 3605 Most users should employ the simplified SLES interface for linear solvers 3606 instead of working directly with matrix algebra routines such as this. 3607 See, e.g., SLESCreate(). 3608 3609 Level: developer 3610 3611 Concepts: matrices^symbolic LU factorization 3612 Concepts: matrices^factorization 3613 Concepts: LU^symbolic factorization 3614 3615 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 3616 MatGetOrdering(), MatILUInfo 3617 3618 @*/ 3619 int MatILUFactorSymbolic(Mat mat,IS row,IS col,MatILUInfo *info,Mat *fact) 3620 { 3621 int ierr; 3622 3623 PetscFunctionBegin; 3624 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3625 PetscValidType(mat); 3626 MatPreallocated(mat); 3627 PetscValidPointer(fact); 3628 PetscValidHeaderSpecific(row,IS_COOKIE); 3629 PetscValidHeaderSpecific(col,IS_COOKIE); 3630 if (info && info->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Levels of fill negative %d",(int)info->levels); 3631 if (info && info->fill < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",info->fill); 3632 if (!mat->ops->ilufactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ILU",mat->type_name); 3633 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3634 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3635 3636 ierr = PetscLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 3637 ierr = (*mat->ops->ilufactorsymbolic)(mat,row,col,info,fact);CHKERRQ(ierr); 3638 ierr = PetscLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0);CHKERRQ(ierr); 3639 PetscFunctionReturn(0); 3640 } 3641 3642 #undef __FUNCT__ 3643 #define __FUNCT__ "MatICCFactorSymbolic" 3644 /*@ 3645 MatICCFactorSymbolic - Performs symbolic incomplete 3646 Cholesky factorization for a symmetric matrix. Use 3647 MatCholeskyFactorNumeric() to complete the factorization. 3648 3649 Collective on Mat 3650 3651 Input Parameters: 3652 + mat - the matrix 3653 . perm - row and column permutation 3654 . fill - levels of fill 3655 - f - expected fill as ratio of original fill 3656 3657 Output Parameter: 3658 . fact - the factored matrix 3659 3660 Notes: 3661 Currently only no-fill factorization is supported. 3662 3663 Most users should employ the simplified SLES interface for linear solvers 3664 instead of working directly with matrix algebra routines such as this. 3665 See, e.g., SLESCreate(). 3666 3667 Level: developer 3668 3669 Concepts: matrices^symbolic incomplete Cholesky factorization 3670 Concepts: matrices^factorization 3671 Concepts: Cholsky^symbolic factorization 3672 3673 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor() 3674 @*/ 3675 int MatICCFactorSymbolic(Mat mat,IS perm,PetscReal f,int fill,Mat *fact) 3676 { 3677 int ierr; 3678 3679 PetscFunctionBegin; 3680 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3681 PetscValidType(mat); 3682 MatPreallocated(mat); 3683 PetscValidPointer(fact); 3684 PetscValidHeaderSpecific(perm,IS_COOKIE); 3685 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3686 if (fill < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Fill negative %d",fill); 3687 if (f < 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Expected fill less than 1.0 %g",f); 3688 if (!mat->ops->iccfactorsymbolic) SETERRQ1(PETSC_ERR_SUP,"Matrix type %s symbolic ICC",mat->type_name); 3689 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3690 3691 ierr = PetscLogEventBegin(MAT_ICCFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 3692 ierr = (*mat->ops->iccfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr); 3693 ierr = PetscLogEventEnd(MAT_ICCFactorSymbolic,mat,perm,0,0);CHKERRQ(ierr); 3694 PetscFunctionReturn(0); 3695 } 3696 3697 #undef __FUNCT__ 3698 #define __FUNCT__ "MatGetArray" 3699 /*@C 3700 MatGetArray - Returns a pointer to the element values in the matrix. 3701 The result of this routine is dependent on the underlying matrix data 3702 structure, and may not even work for certain matrix types. You MUST 3703 call MatRestoreArray() when you no longer need to access the array. 3704 3705 Not Collective 3706 3707 Input Parameter: 3708 . mat - the matrix 3709 3710 Output Parameter: 3711 . v - the location of the values 3712 3713 3714 Fortran Note: 3715 This routine is used differently from Fortran, e.g., 3716 .vb 3717 Mat mat 3718 PetscScalar mat_array(1) 3719 PetscOffset i_mat 3720 int ierr 3721 call MatGetArray(mat,mat_array,i_mat,ierr) 3722 3723 C Access first local entry in matrix; note that array is 3724 C treated as one dimensional 3725 value = mat_array(i_mat + 1) 3726 3727 [... other code ...] 3728 call MatRestoreArray(mat,mat_array,i_mat,ierr) 3729 .ve 3730 3731 See the Fortran chapter of the users manual and 3732 petsc/src/mat/examples/tests for details. 3733 3734 Level: advanced 3735 3736 Concepts: matrices^access array 3737 3738 .seealso: MatRestoreArray(), MatGetArrayF90() 3739 @*/ 3740 int MatGetArray(Mat mat,PetscScalar **v) 3741 { 3742 int ierr; 3743 3744 PetscFunctionBegin; 3745 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3746 PetscValidType(mat); 3747 MatPreallocated(mat); 3748 PetscValidPointer(v); 3749 if (!mat->ops->getarray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3750 ierr = (*mat->ops->getarray)(mat,v);CHKERRQ(ierr); 3751 PetscFunctionReturn(0); 3752 } 3753 3754 #undef __FUNCT__ 3755 #define __FUNCT__ "MatRestoreArray" 3756 /*@C 3757 MatRestoreArray - Restores the matrix after MatGetArray() has been called. 3758 3759 Not Collective 3760 3761 Input Parameter: 3762 + mat - the matrix 3763 - v - the location of the values 3764 3765 Fortran Note: 3766 This routine is used differently from Fortran, e.g., 3767 .vb 3768 Mat mat 3769 PetscScalar mat_array(1) 3770 PetscOffset i_mat 3771 int ierr 3772 call MatGetArray(mat,mat_array,i_mat,ierr) 3773 3774 C Access first local entry in matrix; note that array is 3775 C treated as one dimensional 3776 value = mat_array(i_mat + 1) 3777 3778 [... other code ...] 3779 call MatRestoreArray(mat,mat_array,i_mat,ierr) 3780 .ve 3781 3782 See the Fortran chapter of the users manual and 3783 petsc/src/mat/examples/tests for details 3784 3785 Level: advanced 3786 3787 .seealso: MatGetArray(), MatRestoreArrayF90() 3788 @*/ 3789 int MatRestoreArray(Mat mat,PetscScalar **v) 3790 { 3791 int ierr; 3792 3793 PetscFunctionBegin; 3794 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3795 PetscValidType(mat); 3796 MatPreallocated(mat); 3797 PetscValidPointer(v); 3798 if (!mat->ops->restorearray) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3799 ierr = (*mat->ops->restorearray)(mat,v);CHKERRQ(ierr); 3800 PetscFunctionReturn(0); 3801 } 3802 3803 #undef __FUNCT__ 3804 #define __FUNCT__ "MatGetSubMatrices" 3805 /*@C 3806 MatGetSubMatrices - Extracts several submatrices from a matrix. If submat 3807 points to an array of valid matrices, they may be reused to store the new 3808 submatrices. 3809 3810 Collective on Mat 3811 3812 Input Parameters: 3813 + mat - the matrix 3814 . n - the number of submatrixes to be extracted (on this processor, may be zero) 3815 . irow, icol - index sets of rows and columns to extract 3816 - scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 3817 3818 Output Parameter: 3819 . submat - the array of submatrices 3820 3821 Notes: 3822 MatGetSubMatrices() can extract only sequential submatrices 3823 (from both sequential and parallel matrices). Use MatGetSubMatrix() 3824 to extract a parallel submatrix. 3825 3826 When extracting submatrices from a parallel matrix, each processor can 3827 form a different submatrix by setting the rows and columns of its 3828 individual index sets according to the local submatrix desired. 3829 3830 When finished using the submatrices, the user should destroy 3831 them with MatDestroyMatrices(). 3832 3833 MAT_REUSE_MATRIX can only be used when the nonzero structure of the 3834 original matrix has not changed from that last call to MatGetSubMatrices(). 3835 3836 This routine creates the matrices submat; you should NOT create them before 3837 calling it. 3838 3839 Fortran Note: 3840 The Fortran interface is slightly different from that given below; it 3841 requires one to pass in as submat a Mat (integer) array of size at least m. 3842 3843 Level: advanced 3844 3845 Concepts: matrices^accessing submatrices 3846 Concepts: submatrices 3847 3848 .seealso: MatDestroyMatrices(), MatGetSubMatrix(), MatGetRow(), MatGetDiagonal() 3849 @*/ 3850 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatReuse scall,Mat **submat) 3851 { 3852 int ierr; 3853 3854 PetscFunctionBegin; 3855 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3856 PetscValidType(mat); 3857 MatPreallocated(mat); 3858 if (!mat->ops->getsubmatrices) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3859 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3860 3861 ierr = PetscLogEventBegin(MAT_GetSubMatrices,mat,0,0,0);CHKERRQ(ierr); 3862 ierr = (*mat->ops->getsubmatrices)(mat,n,irow,icol,scall,submat);CHKERRQ(ierr); 3863 ierr = PetscLogEventEnd(MAT_GetSubMatrices,mat,0,0,0);CHKERRQ(ierr); 3864 PetscFunctionReturn(0); 3865 } 3866 3867 #undef __FUNCT__ 3868 #define __FUNCT__ "MatDestroyMatrices" 3869 /*@C 3870 MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices(). 3871 3872 Collective on Mat 3873 3874 Input Parameters: 3875 + n - the number of local matrices 3876 - mat - the matrices 3877 3878 Level: advanced 3879 3880 Notes: Frees not only the matrices, but also the array that contains the matrices 3881 3882 .seealso: MatGetSubMatrices() 3883 @*/ 3884 int MatDestroyMatrices(int n,Mat **mat) 3885 { 3886 int ierr,i; 3887 3888 PetscFunctionBegin; 3889 if (n < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of matrices %d",n); 3890 PetscValidPointer(mat); 3891 for (i=0; i<n; i++) { 3892 ierr = MatDestroy((*mat)[i]);CHKERRQ(ierr); 3893 } 3894 /* memory is allocated even if n = 0 */ 3895 ierr = PetscFree(*mat);CHKERRQ(ierr); 3896 PetscFunctionReturn(0); 3897 } 3898 3899 #undef __FUNCT__ 3900 #define __FUNCT__ "MatIncreaseOverlap" 3901 /*@ 3902 MatIncreaseOverlap - Given a set of submatrices indicated by index sets, 3903 replaces the index sets by larger ones that represent submatrices with 3904 additional overlap. 3905 3906 Collective on Mat 3907 3908 Input Parameters: 3909 + mat - the matrix 3910 . n - the number of index sets 3911 . is - the array of pointers to index sets 3912 - ov - the additional overlap requested 3913 3914 Level: developer 3915 3916 Concepts: overlap 3917 Concepts: ASM^computing overlap 3918 3919 .seealso: MatGetSubMatrices() 3920 @*/ 3921 int MatIncreaseOverlap(Mat mat,int n,IS *is,int ov) 3922 { 3923 int ierr; 3924 3925 PetscFunctionBegin; 3926 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3927 PetscValidType(mat); 3928 MatPreallocated(mat); 3929 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3930 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3931 3932 if (!ov) PetscFunctionReturn(0); 3933 if (!mat->ops->increaseoverlap) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 3934 ierr = PetscLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0);CHKERRQ(ierr); 3935 ierr = (*mat->ops->increaseoverlap)(mat,n,is,ov);CHKERRQ(ierr); 3936 ierr = PetscLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0);CHKERRQ(ierr); 3937 PetscFunctionReturn(0); 3938 } 3939 3940 #undef __FUNCT__ 3941 #define __FUNCT__ "MatPrintHelp" 3942 /*@ 3943 MatPrintHelp - Prints all the options for the matrix. 3944 3945 Collective on Mat 3946 3947 Input Parameter: 3948 . mat - the matrix 3949 3950 Options Database Keys: 3951 + -help - Prints matrix options 3952 - -h - Prints matrix options 3953 3954 Level: developer 3955 3956 .seealso: MatCreate(), MatCreateXXX() 3957 @*/ 3958 int MatPrintHelp(Mat mat) 3959 { 3960 static PetscTruth called = PETSC_FALSE; 3961 int ierr; 3962 MPI_Comm comm; 3963 3964 PetscFunctionBegin; 3965 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3966 PetscValidType(mat); 3967 MatPreallocated(mat); 3968 3969 comm = mat->comm; 3970 if (!called) { 3971 ierr = (*PetscHelpPrintf)(comm,"General matrix options:\n");CHKERRQ(ierr); 3972 ierr = (*PetscHelpPrintf)(comm," -mat_view_info: view basic matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr); 3973 ierr = (*PetscHelpPrintf)(comm," -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n");CHKERRQ(ierr); 3974 ierr = (*PetscHelpPrintf)(comm," -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n");CHKERRQ(ierr); 3975 ierr = (*PetscHelpPrintf)(comm," -draw_pause <sec>: set seconds of display pause\n");CHKERRQ(ierr); 3976 ierr = (*PetscHelpPrintf)(comm," -display <name>: set alternate display\n");CHKERRQ(ierr); 3977 called = PETSC_TRUE; 3978 } 3979 if (mat->ops->printhelp) { 3980 ierr = (*mat->ops->printhelp)(mat);CHKERRQ(ierr); 3981 } 3982 PetscFunctionReturn(0); 3983 } 3984 3985 #undef __FUNCT__ 3986 #define __FUNCT__ "MatGetBlockSize" 3987 /*@ 3988 MatGetBlockSize - Returns the matrix block size; useful especially for the 3989 block row and block diagonal formats. 3990 3991 Not Collective 3992 3993 Input Parameter: 3994 . mat - the matrix 3995 3996 Output Parameter: 3997 . bs - block size 3998 3999 Notes: 4000 Block diagonal formats are MATSEQBDIAG, MATMPIBDIAG. 4001 Block row formats are MATSEQBAIJ, MATMPIBAIJ 4002 4003 Level: intermediate 4004 4005 Concepts: matrices^block size 4006 4007 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag() 4008 @*/ 4009 int MatGetBlockSize(Mat mat,int *bs) 4010 { 4011 int ierr; 4012 4013 PetscFunctionBegin; 4014 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4015 PetscValidType(mat); 4016 MatPreallocated(mat); 4017 PetscValidIntPointer(bs); 4018 if (!mat->ops->getblocksize) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4019 ierr = (*mat->ops->getblocksize)(mat,bs);CHKERRQ(ierr); 4020 PetscFunctionReturn(0); 4021 } 4022 4023 #undef __FUNCT__ 4024 #define __FUNCT__ "MatGetRowIJ" 4025 /*@C 4026 MatGetRowIJ - Returns the compressed row storage i and j indices for sequential matrices. 4027 4028 Collective on Mat 4029 4030 Input Parameters: 4031 + mat - the matrix 4032 . shift - 0 or 1 indicating we want the indices starting at 0 or 1 4033 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4034 symmetrized 4035 4036 Output Parameters: 4037 + n - number of rows in the (possibly compressed) matrix 4038 . ia - the row pointers 4039 . ja - the column indices 4040 - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned 4041 4042 Level: developer 4043 4044 .seealso: MatGetColumnIJ(), MatRestoreRowIJ() 4045 @*/ 4046 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 4047 { 4048 int ierr; 4049 4050 PetscFunctionBegin; 4051 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4052 PetscValidType(mat); 4053 MatPreallocated(mat); 4054 if (ia) PetscValidIntPointer(ia); 4055 if (ja) PetscValidIntPointer(ja); 4056 PetscValidIntPointer(done); 4057 if (!mat->ops->getrowij) *done = PETSC_FALSE; 4058 else { 4059 *done = PETSC_TRUE; 4060 ierr = (*mat->ops->getrowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4061 } 4062 PetscFunctionReturn(0); 4063 } 4064 4065 #undef __FUNCT__ 4066 #define __FUNCT__ "MatGetColumnIJ" 4067 /*@C 4068 MatGetColumnIJ - Returns the compressed column storage i and j indices for sequential matrices. 4069 4070 Collective on Mat 4071 4072 Input Parameters: 4073 + mat - the matrix 4074 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4075 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4076 symmetrized 4077 4078 Output Parameters: 4079 + n - number of columns in the (possibly compressed) matrix 4080 . ia - the column pointers 4081 . ja - the row indices 4082 - done - PETSC_TRUE or PETSC_FALSE, indicating whether the values have been returned 4083 4084 Level: developer 4085 4086 .seealso: MatGetRowIJ(), MatRestoreColumnIJ() 4087 @*/ 4088 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 4089 { 4090 int ierr; 4091 4092 PetscFunctionBegin; 4093 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4094 PetscValidType(mat); 4095 MatPreallocated(mat); 4096 if (ia) PetscValidIntPointer(ia); 4097 if (ja) PetscValidIntPointer(ja); 4098 PetscValidIntPointer(done); 4099 4100 if (!mat->ops->getcolumnij) *done = PETSC_FALSE; 4101 else { 4102 *done = PETSC_TRUE; 4103 ierr = (*mat->ops->getcolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4104 } 4105 PetscFunctionReturn(0); 4106 } 4107 4108 #undef __FUNCT__ 4109 #define __FUNCT__ "MatRestoreRowIJ" 4110 /*@C 4111 MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with 4112 MatGetRowIJ(). 4113 4114 Collective on Mat 4115 4116 Input Parameters: 4117 + mat - the matrix 4118 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4119 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4120 symmetrized 4121 4122 Output Parameters: 4123 + n - size of (possibly compressed) matrix 4124 . ia - the row pointers 4125 . ja - the column indices 4126 - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 4127 4128 Level: developer 4129 4130 .seealso: MatGetRowIJ(), MatRestoreColumnIJ() 4131 @*/ 4132 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 4133 { 4134 int ierr; 4135 4136 PetscFunctionBegin; 4137 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4138 PetscValidType(mat); 4139 MatPreallocated(mat); 4140 if (ia) PetscValidIntPointer(ia); 4141 if (ja) PetscValidIntPointer(ja); 4142 PetscValidIntPointer(done); 4143 4144 if (!mat->ops->restorerowij) *done = PETSC_FALSE; 4145 else { 4146 *done = PETSC_TRUE; 4147 ierr = (*mat->ops->restorerowij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4148 } 4149 PetscFunctionReturn(0); 4150 } 4151 4152 #undef __FUNCT__ 4153 #define __FUNCT__ "MatRestoreColumnIJ" 4154 /*@C 4155 MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with 4156 MatGetColumnIJ(). 4157 4158 Collective on Mat 4159 4160 Input Parameters: 4161 + mat - the matrix 4162 . shift - 1 or zero indicating we want the indices starting at 0 or 1 4163 - symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 4164 symmetrized 4165 4166 Output Parameters: 4167 + n - size of (possibly compressed) matrix 4168 . ia - the column pointers 4169 . ja - the row indices 4170 - done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 4171 4172 Level: developer 4173 4174 .seealso: MatGetColumnIJ(), MatRestoreRowIJ() 4175 @*/ 4176 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 4177 { 4178 int ierr; 4179 4180 PetscFunctionBegin; 4181 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4182 PetscValidType(mat); 4183 MatPreallocated(mat); 4184 if (ia) PetscValidIntPointer(ia); 4185 if (ja) PetscValidIntPointer(ja); 4186 PetscValidIntPointer(done); 4187 4188 if (!mat->ops->restorecolumnij) *done = PETSC_FALSE; 4189 else { 4190 *done = PETSC_TRUE; 4191 ierr = (*mat->ops->restorecolumnij)(mat,shift,symmetric,n,ia,ja,done);CHKERRQ(ierr); 4192 } 4193 PetscFunctionReturn(0); 4194 } 4195 4196 #undef __FUNCT__ 4197 #define __FUNCT__ "MatColoringPatch" 4198 /*@C 4199 MatColoringPatch -Used inside matrix coloring routines that 4200 use MatGetRowIJ() and/or MatGetColumnIJ(). 4201 4202 Collective on Mat 4203 4204 Input Parameters: 4205 + mat - the matrix 4206 . n - number of colors 4207 - colorarray - array indicating color for each column 4208 4209 Output Parameters: 4210 . iscoloring - coloring generated using colorarray information 4211 4212 Level: developer 4213 4214 .seealso: MatGetRowIJ(), MatGetColumnIJ() 4215 4216 @*/ 4217 int MatColoringPatch(Mat mat,int n,int ncolors,int *colorarray,ISColoring *iscoloring) 4218 { 4219 int ierr; 4220 4221 PetscFunctionBegin; 4222 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4223 PetscValidType(mat); 4224 MatPreallocated(mat); 4225 PetscValidIntPointer(colorarray); 4226 4227 if (!mat->ops->coloringpatch){ 4228 ierr = ISColoringCreate(mat->comm,n,colorarray,iscoloring);CHKERRQ(ierr); 4229 } else { 4230 ierr = (*mat->ops->coloringpatch)(mat,n,ncolors,colorarray,iscoloring);CHKERRQ(ierr); 4231 } 4232 PetscFunctionReturn(0); 4233 } 4234 4235 4236 #undef __FUNCT__ 4237 #define __FUNCT__ "MatSetUnfactored" 4238 /*@ 4239 MatSetUnfactored - Resets a factored matrix to be treated as unfactored. 4240 4241 Collective on Mat 4242 4243 Input Parameter: 4244 . mat - the factored matrix to be reset 4245 4246 Notes: 4247 This routine should be used only with factored matrices formed by in-place 4248 factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE 4249 format). This option can save memory, for example, when solving nonlinear 4250 systems with a matrix-free Newton-Krylov method and a matrix-based, in-place 4251 ILU(0) preconditioner. 4252 4253 Note that one can specify in-place ILU(0) factorization by calling 4254 .vb 4255 PCType(pc,PCILU); 4256 PCILUSeUseInPlace(pc); 4257 .ve 4258 or by using the options -pc_type ilu -pc_ilu_in_place 4259 4260 In-place factorization ILU(0) can also be used as a local 4261 solver for the blocks within the block Jacobi or additive Schwarz 4262 methods (runtime option: -sub_pc_ilu_in_place). See the discussion 4263 of these preconditioners in the users manual for details on setting 4264 local solver options. 4265 4266 Most users should employ the simplified SLES interface for linear solvers 4267 instead of working directly with matrix algebra routines such as this. 4268 See, e.g., SLESCreate(). 4269 4270 Level: developer 4271 4272 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace() 4273 4274 Concepts: matrices^unfactored 4275 4276 @*/ 4277 int MatSetUnfactored(Mat mat) 4278 { 4279 int ierr; 4280 4281 PetscFunctionBegin; 4282 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4283 PetscValidType(mat); 4284 MatPreallocated(mat); 4285 mat->factor = 0; 4286 if (!mat->ops->setunfactored) PetscFunctionReturn(0); 4287 ierr = (*mat->ops->setunfactored)(mat);CHKERRQ(ierr); 4288 PetscFunctionReturn(0); 4289 } 4290 4291 /*MC 4292 MatGetArrayF90 - Accesses a matrix array from Fortran90. 4293 4294 Synopsis: 4295 MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 4296 4297 Not collective 4298 4299 Input Parameter: 4300 . x - matrix 4301 4302 Output Parameters: 4303 + xx_v - the Fortran90 pointer to the array 4304 - ierr - error code 4305 4306 Example of Usage: 4307 .vb 4308 PetscScalar, pointer xx_v(:) 4309 .... 4310 call MatGetArrayF90(x,xx_v,ierr) 4311 a = xx_v(3) 4312 call MatRestoreArrayF90(x,xx_v,ierr) 4313 .ve 4314 4315 Notes: 4316 Not yet supported for all F90 compilers 4317 4318 Level: advanced 4319 4320 .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray() 4321 4322 Concepts: matrices^accessing array 4323 4324 M*/ 4325 4326 /*MC 4327 MatRestoreArrayF90 - Restores a matrix array that has been 4328 accessed with MatGetArrayF90(). 4329 4330 Synopsis: 4331 MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 4332 4333 Not collective 4334 4335 Input Parameters: 4336 + x - matrix 4337 - xx_v - the Fortran90 pointer to the array 4338 4339 Output Parameter: 4340 . ierr - error code 4341 4342 Example of Usage: 4343 .vb 4344 PetscScalar, pointer xx_v(:) 4345 .... 4346 call MatGetArrayF90(x,xx_v,ierr) 4347 a = xx_v(3) 4348 call MatRestoreArrayF90(x,xx_v,ierr) 4349 .ve 4350 4351 Notes: 4352 Not yet supported for all F90 compilers 4353 4354 Level: advanced 4355 4356 .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray() 4357 4358 M*/ 4359 4360 4361 #undef __FUNCT__ 4362 #define __FUNCT__ "MatGetSubMatrix" 4363 /*@ 4364 MatGetSubMatrix - Gets a single submatrix on the same number of processors 4365 as the original matrix. 4366 4367 Collective on Mat 4368 4369 Input Parameters: 4370 + mat - the original matrix 4371 . isrow - rows this processor should obtain 4372 . iscol - columns for all processors you wish to keep 4373 . csize - number of columns "local" to this processor (does nothing for sequential 4374 matrices). This should match the result from VecGetLocalSize(x,...) if you 4375 plan to use the matrix in a A*x; alternatively, you can use PETSC_DECIDE 4376 - cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4377 4378 Output Parameter: 4379 . newmat - the new submatrix, of the same type as the old 4380 4381 Level: advanced 4382 4383 Notes: the iscol argument MUST be the same on each processor. 4384 4385 The first time this is called you should use a cll of MAT_INITIAL_MATRIX, 4386 the MatGetSubMatrix() routine will create the newmat for you. Any additional calls 4387 to this routine with a mat of the same nonzero structure will reuse the matrix 4388 generated the first time. 4389 4390 Concepts: matrices^submatrices 4391 4392 .seealso: MatGetSubMatrices() 4393 @*/ 4394 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,Mat *newmat) 4395 { 4396 int ierr, size; 4397 Mat *local; 4398 4399 PetscFunctionBegin; 4400 PetscValidType(mat); 4401 MatPreallocated(mat); 4402 ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr); 4403 4404 /* if original matrix is on just one processor then use submatrix generated */ 4405 if (!mat->ops->getsubmatrix && size == 1 && cll == MAT_REUSE_MATRIX) { 4406 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr); 4407 PetscFunctionReturn(0); 4408 } else if (!mat->ops->getsubmatrix && size == 1) { 4409 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr); 4410 *newmat = *local; 4411 ierr = PetscFree(local);CHKERRQ(ierr); 4412 PetscFunctionReturn(0); 4413 } 4414 4415 if (!mat->ops->getsubmatrix) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4416 ierr = (*mat->ops->getsubmatrix)(mat,isrow,iscol,csize,cll,newmat);CHKERRQ(ierr); 4417 PetscFunctionReturn(0); 4418 } 4419 4420 #undef __FUNCT__ 4421 #define __FUNCT__ "MatGetPetscMaps" 4422 /*@C 4423 MatGetPetscMaps - Returns the maps associated with the matrix. 4424 4425 Not Collective 4426 4427 Input Parameter: 4428 . mat - the matrix 4429 4430 Output Parameters: 4431 + rmap - the row (right) map 4432 - cmap - the column (left) map 4433 4434 Level: developer 4435 4436 Concepts: maps^getting from matrix 4437 4438 @*/ 4439 int MatGetPetscMaps(Mat mat,PetscMap *rmap,PetscMap *cmap) 4440 { 4441 int ierr; 4442 4443 PetscFunctionBegin; 4444 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4445 PetscValidType(mat); 4446 MatPreallocated(mat); 4447 ierr = (*mat->ops->getmaps)(mat,rmap,cmap);CHKERRQ(ierr); 4448 PetscFunctionReturn(0); 4449 } 4450 4451 /* 4452 Version that works for all PETSc matrices 4453 */ 4454 #undef __FUNCT__ 4455 #define __FUNCT__ "MatGetPetscMaps_Petsc" 4456 int MatGetPetscMaps_Petsc(Mat mat,PetscMap *rmap,PetscMap *cmap) 4457 { 4458 PetscFunctionBegin; 4459 if (rmap) *rmap = mat->rmap; 4460 if (cmap) *cmap = mat->cmap; 4461 PetscFunctionReturn(0); 4462 } 4463 4464 #undef __FUNCT__ 4465 #define __FUNCT__ "MatSetStashInitialSize" 4466 /*@ 4467 MatSetStashInitialSize - sets the sizes of the matrix stash, that is 4468 used during the assembly process to store values that belong to 4469 other processors. 4470 4471 Not Collective 4472 4473 Input Parameters: 4474 + mat - the matrix 4475 . size - the initial size of the stash. 4476 - bsize - the initial size of the block-stash(if used). 4477 4478 Options Database Keys: 4479 + -matstash_initial_size <size> or <size0,size1,...sizep-1> 4480 - -matstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1> 4481 4482 Level: intermediate 4483 4484 Notes: 4485 The block-stash is used for values set with VecSetValuesBlocked() while 4486 the stash is used for values set with VecSetValues() 4487 4488 Run with the option -log_info and look for output of the form 4489 MatAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs. 4490 to determine the appropriate value, MM, to use for size and 4491 MatAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs. 4492 to determine the value, BMM to use for bsize 4493 4494 Concepts: stash^setting matrix size 4495 Concepts: matrices^stash 4496 4497 @*/ 4498 int MatSetStashInitialSize(Mat mat,int size, int bsize) 4499 { 4500 int ierr; 4501 4502 PetscFunctionBegin; 4503 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4504 PetscValidType(mat); 4505 MatPreallocated(mat); 4506 ierr = MatStashSetInitialSize_Private(&mat->stash,size);CHKERRQ(ierr); 4507 ierr = MatStashSetInitialSize_Private(&mat->bstash,bsize);CHKERRQ(ierr); 4508 PetscFunctionReturn(0); 4509 } 4510 4511 #undef __FUNCT__ 4512 #define __FUNCT__ "MatInterpolateAdd" 4513 /*@ 4514 MatInterpolateAdd - w = y + A*x or A'*x depending on the shape of 4515 the matrix 4516 4517 Collective on Mat 4518 4519 Input Parameters: 4520 + mat - the matrix 4521 . x,y - the vectors 4522 - w - where the result is stored 4523 4524 Level: intermediate 4525 4526 Notes: 4527 w may be the same vector as y. 4528 4529 This allows one to use either the restriction or interpolation (its transpose) 4530 matrix to do the interpolation 4531 4532 Concepts: interpolation 4533 4534 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict() 4535 4536 @*/ 4537 int MatInterpolateAdd(Mat A,Vec x,Vec y,Vec w) 4538 { 4539 int M,N,ierr; 4540 4541 PetscFunctionBegin; 4542 PetscValidType(A); 4543 MatPreallocated(A); 4544 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4545 if (N > M) { 4546 ierr = MatMultTransposeAdd(A,x,y,w);CHKERRQ(ierr); 4547 } else { 4548 ierr = MatMultAdd(A,x,y,w);CHKERRQ(ierr); 4549 } 4550 PetscFunctionReturn(0); 4551 } 4552 4553 #undef __FUNCT__ 4554 #define __FUNCT__ "MatInterpolate" 4555 /*@ 4556 MatInterpolate - y = A*x or A'*x depending on the shape of 4557 the matrix 4558 4559 Collective on Mat 4560 4561 Input Parameters: 4562 + mat - the matrix 4563 - x,y - the vectors 4564 4565 Level: intermediate 4566 4567 Notes: 4568 This allows one to use either the restriction or interpolation (its transpose) 4569 matrix to do the interpolation 4570 4571 Concepts: matrices^interpolation 4572 4573 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatRestrict() 4574 4575 @*/ 4576 int MatInterpolate(Mat A,Vec x,Vec y) 4577 { 4578 int M,N,ierr; 4579 4580 PetscFunctionBegin; 4581 PetscValidType(A); 4582 MatPreallocated(A); 4583 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4584 if (N > M) { 4585 ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr); 4586 } else { 4587 ierr = MatMult(A,x,y);CHKERRQ(ierr); 4588 } 4589 PetscFunctionReturn(0); 4590 } 4591 4592 #undef __FUNCT__ 4593 #define __FUNCT__ "MatRestrict" 4594 /*@ 4595 MatRestrict - y = A*x or A'*x 4596 4597 Collective on Mat 4598 4599 Input Parameters: 4600 + mat - the matrix 4601 - x,y - the vectors 4602 4603 Level: intermediate 4604 4605 Notes: 4606 This allows one to use either the restriction or interpolation (its transpose) 4607 matrix to do the restriction 4608 4609 Concepts: matrices^restriction 4610 4611 .seealso: MatMultAdd(), MatMultTransposeAdd(), MatInterpolate() 4612 4613 @*/ 4614 int MatRestrict(Mat A,Vec x,Vec y) 4615 { 4616 int M,N,ierr; 4617 4618 PetscFunctionBegin; 4619 PetscValidType(A); 4620 MatPreallocated(A); 4621 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 4622 if (N > M) { 4623 ierr = MatMult(A,x,y);CHKERRQ(ierr); 4624 } else { 4625 ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr); 4626 } 4627 PetscFunctionReturn(0); 4628 } 4629 4630 #undef __FUNCT__ 4631 #define __FUNCT__ "MatNullSpaceAttach" 4632 /*@C 4633 MatNullSpaceAttach - attaches a null space to a matrix. 4634 This null space will be removed from the resulting vector whenever 4635 MatMult() is called 4636 4637 Collective on Mat 4638 4639 Input Parameters: 4640 + mat - the matrix 4641 - nullsp - the null space object 4642 4643 Level: developer 4644 4645 Notes: 4646 Overwrites any previous null space that may have been attached 4647 4648 Concepts: null space^attaching to matrix 4649 4650 .seealso: MatCreate(), MatNullSpaceCreate() 4651 @*/ 4652 int MatNullSpaceAttach(Mat mat,MatNullSpace nullsp) 4653 { 4654 int ierr; 4655 4656 PetscFunctionBegin; 4657 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4658 PetscValidType(mat); 4659 MatPreallocated(mat); 4660 PetscValidHeaderSpecific(nullsp,MAT_NULLSPACE_COOKIE); 4661 4662 if (mat->nullsp) { 4663 ierr = MatNullSpaceDestroy(mat->nullsp);CHKERRQ(ierr); 4664 } 4665 mat->nullsp = nullsp; 4666 ierr = PetscObjectReference((PetscObject)nullsp);CHKERRQ(ierr); 4667 PetscFunctionReturn(0); 4668 } 4669 4670 #undef __FUNCT__ 4671 #define __FUNCT__ "MatICCFactor" 4672 /*@ 4673 MatICCFactor - Performs in-place incomplete Cholesky factorization of matrix. 4674 4675 Collective on Mat 4676 4677 Input Parameters: 4678 + mat - the matrix 4679 . row - row/column permutation 4680 . fill - expected fill factor >= 1.0 4681 - level - level of fill, for ICC(k) 4682 4683 Notes: 4684 Probably really in-place only when level of fill is zero, otherwise allocates 4685 new space to store factored matrix and deletes previous memory. 4686 4687 Most users should employ the simplified SLES interface for linear solvers 4688 instead of working directly with matrix algebra routines such as this. 4689 See, e.g., SLESCreate(). 4690 4691 Level: developer 4692 4693 Concepts: matrices^incomplete Cholesky factorization 4694 Concepts: Cholesky factorization 4695 4696 .seealso: MatICCFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 4697 @*/ 4698 int MatICCFactor(Mat mat,IS row,PetscReal fill,int level) 4699 { 4700 int ierr; 4701 4702 PetscFunctionBegin; 4703 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4704 PetscValidType(mat); 4705 MatPreallocated(mat); 4706 if (mat->M != mat->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square"); 4707 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4708 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4709 if (!mat->ops->iccfactor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4710 ierr = (*mat->ops->iccfactor)(mat,row,fill,level);CHKERRQ(ierr); 4711 PetscFunctionReturn(0); 4712 } 4713 4714 #undef __FUNCT__ 4715 #define __FUNCT__ "MatSetValuesAdic" 4716 /*@ 4717 MatSetValuesAdic - Sets values computed with ADIC automatic differentiation into a matrix. 4718 4719 Not Collective 4720 4721 Input Parameters: 4722 + mat - the matrix 4723 - v - the values compute with ADIC 4724 4725 Level: developer 4726 4727 Notes: 4728 Must call MatSetColoring() before using this routine. Also this matrix must already 4729 have its nonzero pattern determined. 4730 4731 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4732 MatSetValues(), MatSetColoring(), MatSetValuesAdifor() 4733 @*/ 4734 int MatSetValuesAdic(Mat mat,void *v) 4735 { 4736 int ierr; 4737 4738 PetscFunctionBegin; 4739 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4740 PetscValidType(mat); 4741 4742 if (!mat->assembled) { 4743 SETERRQ(1,"Matrix must be already assembled"); 4744 } 4745 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4746 if (!mat->ops->setvaluesadic) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4747 ierr = (*mat->ops->setvaluesadic)(mat,v);CHKERRQ(ierr); 4748 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4749 ierr = MatView_Private(mat);CHKERRQ(ierr); 4750 PetscFunctionReturn(0); 4751 } 4752 4753 4754 #undef __FUNCT__ 4755 #define __FUNCT__ "MatSetColoring" 4756 /*@ 4757 MatSetColoring - Sets a coloring used by calls to MatSetValuesAdic() 4758 4759 Not Collective 4760 4761 Input Parameters: 4762 + mat - the matrix 4763 - coloring - the coloring 4764 4765 Level: developer 4766 4767 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4768 MatSetValues(), MatSetValuesAdic() 4769 @*/ 4770 int MatSetColoring(Mat mat,ISColoring coloring) 4771 { 4772 int ierr; 4773 4774 PetscFunctionBegin; 4775 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4776 PetscValidType(mat); 4777 4778 if (!mat->assembled) { 4779 SETERRQ(1,"Matrix must be already assembled"); 4780 } 4781 if (!mat->ops->setcoloring) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4782 ierr = (*mat->ops->setcoloring)(mat,coloring);CHKERRQ(ierr); 4783 PetscFunctionReturn(0); 4784 } 4785 4786 #undef __FUNCT__ 4787 #define __FUNCT__ "MatSetValuesAdifor" 4788 /*@ 4789 MatSetValuesAdifor - Sets values computed with automatic differentiation into a matrix. 4790 4791 Not Collective 4792 4793 Input Parameters: 4794 + mat - the matrix 4795 . nl - leading dimension of v 4796 - v - the values compute with ADIFOR 4797 4798 Level: developer 4799 4800 Notes: 4801 Must call MatSetColoring() before using this routine. Also this matrix must already 4802 have its nonzero pattern determined. 4803 4804 .seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked(), MatSetValuesLocal(), 4805 MatSetValues(), MatSetColoring() 4806 @*/ 4807 int MatSetValuesAdifor(Mat mat,int nl,void *v) 4808 { 4809 int ierr; 4810 4811 PetscFunctionBegin; 4812 PetscValidHeaderSpecific(mat,MAT_COOKIE); 4813 PetscValidType(mat); 4814 4815 if (!mat->assembled) { 4816 SETERRQ(1,"Matrix must be already assembled"); 4817 } 4818 ierr = PetscLogEventBegin(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4819 if (!mat->ops->setvaluesadifor) SETERRQ1(PETSC_ERR_SUP,"Mat type %s",mat->type_name); 4820 ierr = (*mat->ops->setvaluesadifor)(mat,nl,v);CHKERRQ(ierr); 4821 ierr = PetscLogEventEnd(MAT_SetValues,mat,0,0,0);CHKERRQ(ierr); 4822 PetscFunctionReturn(0); 4823 } 4824