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