1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: matrix.c,v 1.266 1997/10/22 15:07:53 bsmith Exp bsmith $"; 3 #endif 4 5 /* 6 This is where the abstract matrix operations are defined 7 */ 8 9 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 10 #include "src/vec/vecimpl.h" 11 #include "pinclude/pviewer.h" 12 13 #undef __FUNC__ 14 #define __FUNC__ "MatGetRow" 15 /*@C 16 MatGetRow - Gets a row of a matrix. You MUST call MatRestoreRow() 17 for each row that you get to ensure that your application does 18 not bleed memory. 19 20 Input Parameters: 21 . mat - the matrix 22 . row - the row to get 23 24 Output Parameters: 25 . ncols - the number of nonzeros in the row 26 . cols - if nonzero, the column numbers 27 . vals - if nonzero, the values 28 29 Notes: 30 This routine is provided for people who need to have direct access 31 to the structure of a matrix. We hope that we provide enough 32 high-level matrix routines that few users will need it. 33 34 For better efficiency, set cols and/or vals to PETSC_NULL if you do 35 not wish to extract these quantities. 36 37 The user can only examine the values extracted with MatGetRow(); 38 the values cannot be altered. To change the matrix entries, one 39 must use MatSetValues(). 40 41 Fortran Notes: 42 $ The calling sequence from Fortran is 43 $ 44 $ MatGetRow(matrix,row,ncols,cols,values,ierr) 45 $ Mat matrix (input) 46 $ integer row (input) 47 $ integer ncols (output) 48 $ integer cols(maxcols) (output) 49 $ double precision (or double complex) values(maxcols) output 50 $ Where maxcols is larger then or equal to the maximum nonzeros in 51 $ any row of the matrix. 52 53 Caution: 54 Do not try to change the contents of the output arrays (cols and vals). 55 In some cases, this may corrupt the matrix. 56 57 .keywords: matrix, row, get, extract 58 59 .seealso: MatRestoreRow(), MatSetValues() 60 @*/ 61 int MatGetRow(Mat mat,int row,int *ncols,int **cols,Scalar **vals) 62 { 63 int ierr; 64 65 PetscFunctionBegin; 66 PetscValidHeaderSpecific(mat,MAT_COOKIE); 67 PetscValidIntPointer(ncols); 68 if (!mat->assembled) SETERRQ(1,0,"Not for unassembled matrix"); 69 if (mat->factor) SETERRQ(1,0,"Not for factored matrix"); 70 if (!mat->ops.getrow) SETERRQ(PETSC_ERR_SUP,0,""); 71 PLogEventBegin(MAT_GetRow,mat,0,0,0); 72 ierr = (*mat->ops.getrow)(mat,row,ncols,cols,vals); CHKERRQ(ierr); 73 PLogEventEnd(MAT_GetRow,mat,0,0,0); 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNC__ 78 #define __FUNC__ "MatRestoreRow" 79 /*@C 80 MatRestoreRow - Frees any temporary space allocated by MatGetRow(). 81 82 Input Parameters: 83 . mat - the matrix 84 . row - the row to get 85 . ncols, cols - the number of nonzeros and their columns 86 . vals - if nonzero the column values 87 88 Fortran Notes: 89 $ The calling sequence from Fortran is 90 $ 91 $ MatRestoreRow(matrix,row,ncols,cols,values,ierr) 92 $ Mat matrix (input) 93 $ integer row (input) 94 $ integer ncols (output) 95 $ integer cols(maxcols) (output) 96 $ double precision (or double complex) values(maxcols) output 97 $ Where maxcols is larger then or equal to the maximum nonzeros in 98 $ any row of the matrix. In Fortran you must call MatRestoreRow() 99 $ after a call to MatGetRow() before you make another call to 100 $ MatGetRow() 101 102 .keywords: matrix, row, restore 103 104 .seealso: MatGetRow() 105 @*/ 106 int MatRestoreRow(Mat mat,int row,int *ncols,int **cols,Scalar **vals) 107 { 108 int ierr; 109 110 PetscFunctionBegin; 111 PetscValidHeaderSpecific(mat,MAT_COOKIE); 112 PetscValidIntPointer(ncols); 113 if (!mat->assembled) SETERRQ(1,0,"Not for unassembled matrix"); 114 if (!mat->ops.restorerow) PetscFunctionReturn(0); 115 ierr = (*mat->ops.restorerow)(mat,row,ncols,cols,vals);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNC__ 120 #define __FUNC__ "MatView" 121 /*@C 122 MatView - Visualizes a matrix object. 123 124 Input Parameters: 125 . mat - the matrix 126 . ptr - visualization context 127 128 Notes: 129 The available visualization contexts include 130 $ VIEWER_STDOUT_SELF - standard output (default) 131 $ VIEWER_STDOUT_WORLD - synchronized standard 132 $ output where only the first processor opens 133 $ the file. All other processors send their 134 $ data to the first processor to print. 135 $ VIEWER_DRAWX_WORLD - graphical display of nonzero structure 136 137 The user can open alternative vistualization contexts with 138 $ ViewerFileOpenASCII() - output to a specified file 139 $ ViewerFileOpenBinary() - output in binary to a 140 $ specified file; corresponding input uses MatLoad() 141 $ ViewerDrawOpenX() - output nonzero matrix structure to 142 $ an X window display 143 $ ViewerMatlabOpen() - output matrix to Matlab viewer. 144 $ Currently only the sequential dense and AIJ 145 $ matrix types support the Matlab viewer. 146 147 The user can call ViewerSetFormat() to specify the output 148 format of ASCII printed objects (when using VIEWER_STDOUT_SELF, 149 VIEWER_STDOUT_WORLD and ViewerFileOpenASCII). Available formats include 150 $ VIEWER_FORMAT_ASCII_DEFAULT - default, prints matrix contents 151 $ VIEWER_FORMAT_ASCII_MATLAB - Matlab format 152 $ VIEWER_FORMAT_ASCII_IMPL - implementation-specific format 153 $ (which is in many cases the same as the default) 154 $ VIEWER_FORMAT_ASCII_INFO - basic information about the matrix 155 $ size and structure (not the matrix entries) 156 $ VIEWER_FORMAT_ASCII_INFO_LONG - more detailed information about the 157 $ matrix structure 158 159 .keywords: matrix, view, visualize, output, print, write, draw 160 161 .seealso: ViewerSetFormat(), ViewerFileOpenASCII(), ViewerDrawOpenX(), 162 ViewerMatlabOpen(), ViewerFileOpenBinary(), MatLoad() 163 @*/ 164 int MatView(Mat mat,Viewer viewer) 165 { 166 int format, ierr, rows, cols; 167 FILE *fd; 168 char *cstr; 169 ViewerType vtype; 170 MPI_Comm comm = mat->comm; 171 172 PetscFunctionBegin; 173 PetscValidHeaderSpecific(mat,MAT_COOKIE); 174 if (viewer) PetscValidHeaderSpecific(viewer,VIEWER_COOKIE); 175 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 176 177 if (!viewer) { 178 viewer = VIEWER_STDOUT_SELF; 179 } 180 181 ierr = ViewerGetType(viewer,&vtype); 182 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 183 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 184 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 185 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 186 PetscFPrintf(comm,fd,"Matrix Object:\n"); 187 ierr = MatGetType(mat,PETSC_NULL,&cstr); CHKERRQ(ierr); 188 ierr = MatGetSize(mat,&rows,&cols); CHKERRQ(ierr); 189 PetscFPrintf(comm,fd," type=%s, rows=%d, cols=%d\n",cstr,rows,cols); 190 if (mat->ops.getinfo) { 191 MatInfo info; 192 ierr = MatGetInfo(mat,MAT_GLOBAL_SUM,&info); CHKERRQ(ierr); 193 PetscFPrintf(comm,fd," total: nonzeros=%d, allocated nonzeros=%d\n", 194 (int)info.nz_used,(int)info.nz_allocated); 195 } 196 } 197 } 198 if (mat->view) {ierr = (*mat->view)((PetscObject)mat,viewer); CHKERRQ(ierr);} 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNC__ 203 #define __FUNC__ "MatDestroy" 204 /*@C 205 MatDestroy - Frees space taken by a matrix. 206 207 Input Parameter: 208 . mat - the matrix 209 210 .keywords: matrix, destroy 211 @*/ 212 int MatDestroy(Mat mat) 213 { 214 int ierr; 215 216 PetscFunctionBegin; 217 PetscValidHeaderSpecific(mat,MAT_COOKIE); 218 if (--mat->refct > 0) PetscFunctionReturn(0); 219 220 if (mat->mapping) { 221 ierr = ISLocalToGlobalMappingDestroy(mat->mapping); CHKERRQ(ierr); 222 } 223 if (mat->bmapping) { 224 ierr = ISLocalToGlobalMappingDestroy(mat->bmapping); CHKERRQ(ierr); 225 } 226 ierr = (*mat->destroy)((PetscObject)mat); CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 #undef __FUNC__ 231 #define __FUNC__ "MatValid" 232 /*@ 233 MatValid - Checks whether a matrix object is valid. 234 235 Input Parameter: 236 . m - the matrix to check 237 238 Output Parameter: 239 flg - flag indicating matrix status, eMAT_IGNORE_OFF_PROC_ENTRIESither 240 $ PETSC_TRUE if matrix is valid; 241 $ PETSC_FALSE otherwise. 242 243 .keywords: matrix, valid 244 @*/ 245 int MatValid(Mat m,PetscTruth *flg) 246 { 247 PetscFunctionBegin; 248 PetscValidIntPointer(flg); 249 if (!m) *flg = PETSC_FALSE; 250 else if (m->cookie != MAT_COOKIE) *flg = PETSC_FALSE; 251 else *flg = PETSC_TRUE; 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNC__ 256 #define __FUNC__ "MatSetValues" 257 /*@ 258 MatSetValues - Inserts or adds a block of values into a matrix. 259 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 260 MUST be called after all calls to MatSetValues() have been completed. 261 262 Input Parameters: 263 . mat - the matrix 264 . v - a logically two-dimensional array of values 265 . m, idxm - the number of rows and their global indices 266 . n, idxn - the number of columns and their global indices 267 . addv - either ADD_VALUES or INSERT_VALUES, where 268 $ ADD_VALUES - adds values to any existing entries 269 $ INSERT_VALUES - replaces existing entries with new values 270 271 Notes: 272 By default the values, v, are row-oriented and unsorted. 273 See MatSetOptions() for other options. 274 275 Calls to MatSetValues() with the INSERT_VALUES and ADD_VALUES 276 options cannot be mixed without intervening calls to the assembly 277 routines. 278 279 MatSetValues() uses 0-based row and column numbers in Fortran 280 as well as in C. 281 282 Efficiency Alert: 283 The routine MatSetValuesBlocked() may offer much better efficiency 284 for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ). 285 286 .keywords: matrix, insert, add, set, values 287 288 .seealso: MatSetOptions(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesBlocked() 289 @*/ 290 int MatSetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv) 291 { 292 int ierr; 293 294 PetscFunctionBegin; 295 if (!m || !n) PetscFunctionReturn(0); /* no values to insert */ 296 PetscValidHeaderSpecific(mat,MAT_COOKIE); 297 PetscValidIntPointer(idxm); 298 PetscValidIntPointer(idxn); 299 PetscValidScalarPointer(v); 300 if (mat->insertmode == NOT_SET_VALUES) { 301 mat->insertmode = addv; 302 } 303 #if defined(USE_PETSC_BOPT_g) 304 else if (mat->insertmode != addv) { 305 SETERRQ(1,1,"Cannot mix add values and insert values"); 306 } 307 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 308 #endif 309 310 if (mat->assembled) { 311 mat->was_assembled = PETSC_TRUE; 312 mat->assembled = PETSC_FALSE; 313 } 314 PLogEventBegin(MAT_SetValues,mat,0,0,0); 315 ierr = (*mat->ops.setvalues)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr); 316 PLogEventEnd(MAT_SetValues,mat,0,0,0); 317 PetscFunctionReturn(0); 318 } 319 320 #undef __FUNC__ 321 #define __FUNC__ "MatSetValuesBlocked" 322 /*@ 323 MatSetValuesBlocked - Inserts or adds a block of values into a matrix. 324 325 Input Parameters: 326 . mat - the matrix 327 . v - a logically two-dimensional array of values 328 . m, idxm - the number of block rows and their global block indices 329 . n, idxn - the number of block columns and their global block indices 330 . addv - either ADD_VALUES or INSERT_VALUES, where 331 $ ADD_VALUES - adds values to any existing entries 332 $ INSERT_VALUES - replaces existing entries with new values 333 334 Notes: 335 By default the values, v, are row-oriented and unsorted. So the layout of 336 v is the same as for MatSetValues(). See MatSetOptions() for other options. 337 338 Calls to MatSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES 339 options cannot be mixed without intervening calls to the assembly 340 routines. 341 342 MatSetValuesBlocked() uses 0-based row and column numbers in Fortran 343 as well as in C. 344 345 Each time an entry is set within a sparse matrix via MatSetValues(), 346 internal searching must be done to determine where to place the the 347 data in the matrix storage space. By instead inserting blocks of 348 entries via MatSetValuesBlocked(), the overhead of matrix assembly is 349 reduced. 350 351 Restrictions: 352 MatSetValuesBlocked() is currently supported only for the block AIJ 353 matrix format (MATSEQBAIJ and MATMPIBAIJ, which are created via 354 MatCreateSeqBAIJ() and MatCreateMPIBAIJ()). 355 356 .keywords: matrix, insert, add, set, values 357 358 .seealso: MatSetOptions(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues() 359 @*/ 360 int MatSetValuesBlocked(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v,InsertMode addv) 361 { 362 int ierr; 363 364 PetscFunctionBegin; 365 if (!m || !n) PetscFunctionReturn(0); /* no values to insert */ 366 PetscValidHeaderSpecific(mat,MAT_COOKIE); 367 PetscValidIntPointer(idxm); 368 PetscValidIntPointer(idxn); 369 PetscValidScalarPointer(v); 370 if (mat->insertmode == NOT_SET_VALUES) { 371 mat->insertmode = addv; 372 } 373 #if defined(USE_PETSC_BOPT_g) 374 else if (mat->insertmode != addv) { 375 SETERRQ(1,1,"Cannot mix add values and insert values"); 376 } 377 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 378 #endif 379 380 if (mat->assembled) { 381 mat->was_assembled = PETSC_TRUE; 382 mat->assembled = PETSC_FALSE; 383 } 384 PLogEventBegin(MAT_SetValues,mat,0,0,0); 385 ierr = (*mat->ops.setvaluesblocked)(mat,m,idxm,n,idxn,v,addv);CHKERRQ(ierr); 386 PLogEventEnd(MAT_SetValues,mat,0,0,0); 387 PetscFunctionReturn(0); 388 } 389 390 /*MC 391 MatSetValue - Set a single entry into a matrix. 392 393 Input Parameters: 394 . m - the matrix 395 . row - the row location of the entry 396 . col - the column location of the entry 397 . value - the value to insert 398 . mode - either INSERT_VALUES or ADD_VALUES 399 400 Synopsis: 401 void MatSetValue(Mat m,int row,int col,Scalar value,InsertMode mode); 402 403 Notes: For efficiency one should use MatSetValues() and set 404 several or many values simultaneously. 405 406 .seealso: MatSetValues() 407 M*/ 408 409 #undef __FUNC__ 410 #define __FUNC__ "MatGetValues" 411 /*@ 412 MatGetValues - Gets a block of values from a matrix. 413 414 Input Parameters: 415 . mat - the matrix 416 . v - a logically two-dimensional array for storing the values 417 . m, idxm - the number of rows and their global indices 418 . n, idxn - the number of columns and their global indices 419 420 Notes: 421 The user must allocate space (m*n Scalars) for the values, v. 422 The values, v, are then returned in a row-oriented format, 423 analogous to that used by default in MatSetValues(). 424 425 MatGetValues() uses 0-based row and column numbers in 426 Fortran as well as in C. 427 428 MatGetValues() requires that the matrix has been assembled 429 with MatAssemblyBegin()/MatAssemblyEnd(). Thus, calls to 430 MatSetValues() and MatGetValues() CANNOT be made in succession 431 without intermediate matrix assembly. 432 433 .keywords: matrix, get, values 434 435 .seealso: MatGetRow(), MatGetSubmatrices(), MatSetValues() 436 @*/ 437 int MatGetValues(Mat mat,int m,int *idxm,int n,int *idxn,Scalar *v) 438 { 439 int ierr; 440 441 PetscFunctionBegin; 442 PetscValidHeaderSpecific(mat,MAT_COOKIE); 443 PetscValidIntPointer(idxm); 444 PetscValidIntPointer(idxn); 445 PetscValidScalarPointer(v); 446 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 447 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 448 if (!mat->ops.getvalues) SETERRQ(PETSC_ERR_SUP,0,""); 449 450 PLogEventBegin(MAT_GetValues,mat,0,0,0); 451 ierr = (*mat->ops.getvalues)(mat,m,idxm,n,idxn,v); CHKERRQ(ierr); 452 PLogEventEnd(MAT_GetValues,mat,0,0,0); 453 PetscFunctionReturn(0); 454 } 455 456 #undef __FUNC__ 457 #define __FUNC__ "MatSetLocalToGlobalMapping" 458 /*@ 459 MatSetLocalToGlobalMapping - Sets a local-to-global numbering for use by 460 the routine MatSetValuesLocal() to allow users to insert matrix entries 461 using a local (per-processor) numbering. 462 463 Input Parameters: 464 . x - the matrix 465 . mapping - mapping created with ISLocalToGlobalMappingCreate() 466 or ISLocalToGlobalMappingCreateIS() 467 468 .keywords: matrix, set, values, local, global, mapping 469 470 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesLocal() 471 @*/ 472 int MatSetLocalToGlobalMapping(Mat x,ISLocalToGlobalMapping mapping) 473 { 474 PetscFunctionBegin; 475 PetscValidHeaderSpecific(x,MAT_COOKIE); 476 PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE); 477 478 if (x->mapping) { 479 SETERRQ(1,0,"Mapping already set for matrix"); 480 } 481 482 x->mapping = mapping; 483 PetscObjectReference((PetscObject)mapping); 484 PetscFunctionReturn(0); 485 } 486 487 #undef __FUNC__ 488 #define __FUNC__ "MatSetLocalToGlobalMappingBlocked" 489 /*@ 490 MatSetLocalToGlobalMappingBlocked - Sets a local-to-global numbering for use 491 by the routine MatSetValuesBlockedLocal() to allow users to insert matrix 492 entries using a local (per-processor) numbering. 493 494 Input Parameters: 495 . x - the matrix 496 . mapping - mapping created with ISLocalToGlobalMappingCreate() or 497 ISLocalToGlobalMappingCreateIS() 498 499 .keywords: matrix, set, values, local ordering 500 501 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetValuesBlockedLocal(), 502 MatSetValuesBlocked(), MatSetValuesLocal() 503 @*/ 504 int MatSetLocalToGlobalMappingBlocked(Mat x,ISLocalToGlobalMapping mapping) 505 { 506 PetscFunctionBegin; 507 PetscValidHeaderSpecific(x,MAT_COOKIE); 508 PetscValidHeaderSpecific(mapping,IS_LTOGM_COOKIE); 509 510 if (x->bmapping) { 511 SETERRQ(1,0,"Mapping already set for matrix"); 512 } 513 514 x->bmapping = mapping; 515 PetscObjectReference((PetscObject)mapping); 516 PetscFunctionReturn(0); 517 } 518 519 #undef __FUNC__ 520 #define __FUNC__ "MatSetValuesLocal" 521 /*@ 522 MatSetValuesLocal - Inserts or adds values into certain locations of a matrix, 523 using a local ordering of the nodes. 524 525 Input Parameters: 526 . x - the matrix 527 . nrow, irow - number of rows and their local indices 528 . ncol, icol - number of columns and their local indices 529 . y - a logically two-dimensional array of values 530 . addv - either INSERT_VALUES or ADD_VALUES, where 531 $ ADD_VALUES - adds values to any existing entries 532 $ INSERT_VALUES - replaces existing entries with new values 533 534 Notes: 535 Before calling MatSetValuesLocal(), the user must first set the 536 local-to-global mapping by calling MatSetLocalToGlobalMapping(). 537 538 Calls to MatSetValuesLocal() with the INSERT_VALUES and ADD_VALUES 539 options cannot be mixed without intervening calls to the assembly 540 routines. 541 542 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 543 MUST be called after all calls to MatSetValuesLocal() have been completed. 544 545 .keywords: matrix, set, values, local ordering 546 547 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), MatSetLocalToGlobalMapping() 548 @*/ 549 int MatSetValuesLocal(Mat mat,int nrow,int *irow,int ncol, int *icol,Scalar *y,InsertMode addv) 550 { 551 int ierr,irowm[2048],icolm[2048]; 552 553 PetscFunctionBegin; 554 PetscValidHeaderSpecific(mat,MAT_COOKIE); 555 PetscValidIntPointer(irow); 556 PetscValidIntPointer(icol); 557 PetscValidScalarPointer(y); 558 559 if (mat->insertmode == NOT_SET_VALUES) { 560 mat->insertmode = addv; 561 } 562 #if defined(USE_PETSC_BOPT_g) 563 else if (mat->insertmode != addv) { 564 SETERRQ(1,1,"Cannot mix add values and insert values"); 565 } 566 if (!mat->mapping) { 567 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMapping"); 568 } 569 if (nrow > 2048 || ncol > 2048) { 570 SETERRQ(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048"); 571 } 572 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 573 #endif 574 575 if (mat->assembled) { 576 mat->was_assembled = PETSC_TRUE; 577 mat->assembled = PETSC_FALSE; 578 } 579 PLogEventBegin(MAT_SetValues,mat,0,0,0); 580 ierr = ISLocalToGlobalMappingApply(mat->mapping,nrow,irow,irowm); CHKERRQ(ierr); 581 ierr = ISLocalToGlobalMappingApply(mat->mapping,ncol,icol,icolm);CHKERRQ(ierr); 582 ierr = (*mat->ops.setvalues)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr); 583 PLogEventEnd(MAT_SetValues,mat,0,0,0); 584 PetscFunctionReturn(0); 585 } 586 587 #undef __FUNC__ 588 #define __FUNC__ "MatSetValuesBlockedLocal" 589 /*@ 590 MatSetValuesBlockedLocal - Inserts or adds values into certain locations of a matrix, 591 using a local ordering of the nodes a block at a time. 592 593 Input Parameters: 594 . x - the matrix 595 . nrow, irow - number of rows and their local indices 596 . ncol, icol - number of columns and their local indices 597 . y - a logically two-dimensional array of values 598 . addv - either INSERT_VALUES or ADD_VALUES, where 599 $ ADD_VALUES - adds values to any existing entries 600 $ INSERT_VALUES - replaces existing entries with new values 601 602 Notes: 603 Before calling MatSetValuesBlockedLocal(), the user must first set the 604 local-to-global mapping by calling MatSetLocalToGlobalMappingBlocked(), 605 where the mapping MUST be set for matrix blocks, not for matrix elements. 606 607 Calls to MatSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES 608 options cannot be mixed without intervening calls to the assembly 609 routines. 610 611 These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd() 612 MUST be called after all calls to MatSetValuesBlockedLocal() have been completed. 613 614 .keywords: matrix, set, values, blocked, local 615 616 .seealso: MatAssemblyBegin(), MatAssemblyEnd(), MatSetValuesLocal(), MatSetLocalToGlobalMappingBlocked() 617 @*/ 618 int MatSetValuesBlockedLocal(Mat mat,int nrow,int *irow,int ncol,int *icol,Scalar *y,InsertMode addv) 619 { 620 int ierr,irowm[2048],icolm[2048]; 621 622 PetscFunctionBegin; 623 PetscValidHeaderSpecific(mat,MAT_COOKIE); 624 PetscValidIntPointer(irow); 625 PetscValidIntPointer(icol); 626 PetscValidScalarPointer(y); 627 if (mat->insertmode == NOT_SET_VALUES) { 628 mat->insertmode = addv; 629 } 630 #if defined(USE_PETSC_BOPT_g) 631 else if (mat->insertmode != addv) { 632 SETERRQ(1,1,"Cannot mix add values and insert values"); 633 } 634 if (!mat->bmapping) { 635 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Local to global never set with MatSetLocalToGlobalMappingBlocked"); 636 } 637 if (nrow > 2048 || ncol > 2048) { 638 SETERRQ(PETSC_ERR_SUP,0,"Number column/row indices must be <= 2048"); 639 } 640 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 641 #endif 642 643 if (mat->assembled) { 644 mat->was_assembled = PETSC_TRUE; 645 mat->assembled = PETSC_FALSE; 646 } 647 PLogEventBegin(MAT_SetValues,mat,0,0,0); 648 ierr = ISLocalToGlobalMappingApply(mat->bmapping,nrow,irow,irowm); CHKERRQ(ierr); 649 ierr = ISLocalToGlobalMappingApply(mat->bmapping,ncol,icol,icolm); CHKERRQ(ierr); 650 ierr = (*mat->ops.setvaluesblocked)(mat,nrow,irowm,ncol,icolm,y,addv);CHKERRQ(ierr); 651 PLogEventEnd(MAT_SetValues,mat,0,0,0); 652 PetscFunctionReturn(0); 653 } 654 655 /* --------------------------------------------------------*/ 656 #undef __FUNC__ 657 #define __FUNC__ "MatMult" 658 /*@ 659 MatMult - Computes the matrix-vector product, y = Ax. 660 661 Input Parameters: 662 . mat - the matrix 663 . x - the vector to be multilplied 664 665 Output Parameters: 666 . y - the result 667 668 Notes: 669 The vectors x and y cannot be the same. I.e., one cannot 670 call MatMult(A,y,y). 671 672 .keywords: matrix, multiply, matrix-vector product 673 674 .seealso: MatMultTrans(), MatMultAdd(), MatMultTransAdd() 675 @*/ 676 int MatMult(Mat mat,Vec x,Vec y) 677 { 678 int ierr; 679 680 PetscFunctionBegin; 681 PetscValidHeaderSpecific(mat,MAT_COOKIE); 682 PetscValidHeaderSpecific(x,VEC_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 683 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 684 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 685 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors"); 686 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 687 if (mat->M != y->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim"); 688 if (mat->m != y->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: local dim"); 689 690 PLogEventBegin(MAT_Mult,mat,x,y,0); 691 ierr = (*mat->ops.mult)(mat,x,y); CHKERRQ(ierr); 692 PLogEventEnd(MAT_Mult,mat,x,y,0); 693 694 PetscFunctionReturn(0); 695 } 696 697 #undef __FUNC__ 698 #define __FUNC__ "MatMultTrans" 699 /*@ 700 MatMultTrans - Computes matrix transpose times a vector. 701 702 Input Parameters: 703 . mat - the matrix 704 . x - the vector to be multilplied 705 706 Output Parameters: 707 . y - the result 708 709 Notes: 710 The vectors x and y cannot be the same. I.e., one cannot 711 call MatMultTrans(A,y,y). 712 713 .keywords: matrix, multiply, matrix-vector product, transpose 714 715 .seealso: MatMult(), MatMultAdd(), MatMultTransAdd() 716 @*/ 717 int MatMultTrans(Mat mat,Vec x,Vec y) 718 { 719 int ierr; 720 721 PetscFunctionBegin; 722 PetscValidHeaderSpecific(mat,MAT_COOKIE); 723 PetscValidHeaderSpecific(x,VEC_COOKIE); PetscValidHeaderSpecific(y,VEC_COOKIE); 724 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 725 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 726 if (x == y) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"x and y must be different vectors"); 727 if (mat->M != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 728 if (mat->N != y->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim"); 729 PLogEventBegin(MAT_MultTrans,mat,x,y,0); 730 ierr = (*mat->ops.multtrans)(mat,x,y); CHKERRQ(ierr); 731 PLogEventEnd(MAT_MultTrans,mat,x,y,0); 732 PetscFunctionReturn(0); 733 } 734 735 #undef __FUNC__ 736 #define __FUNC__ "MatMultAdd" 737 /*@ 738 MatMultAdd - Computes v3 = v2 + A * v1. 739 740 Input Parameters: 741 . mat - the matrix 742 . v1, v2 - the vectors 743 744 Output Parameters: 745 . v3 - the result 746 747 Notes: 748 The vectors v1 and v3 cannot be the same. I.e., one cannot 749 call MatMultAdd(A,v1,v2,v1). 750 751 .keywords: matrix, multiply, matrix-vector product, add 752 753 .seealso: MatMultTrans(), MatMult(), MatMultTransAdd() 754 @*/ 755 int MatMultAdd(Mat mat,Vec v1,Vec v2,Vec v3) 756 { 757 int ierr; 758 759 PetscFunctionBegin; 760 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v1,VEC_COOKIE); 761 PetscValidHeaderSpecific(v2,VEC_COOKIE); PetscValidHeaderSpecific(v3,VEC_COOKIE); 762 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 763 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 764 if (mat->N != v1->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim"); 765 if (mat->M != v2->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim"); 766 if (mat->M != v3->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim"); 767 if (mat->m != v3->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: local dim"); 768 if (mat->m != v2->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: local dim"); 769 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors"); 770 771 PLogEventBegin(MAT_MultAdd,mat,v1,v2,v3); 772 ierr = (*mat->ops.multadd)(mat,v1,v2,v3); CHKERRQ(ierr); 773 PLogEventEnd(MAT_MultAdd,mat,v1,v2,v3); 774 PetscFunctionReturn(0); 775 } 776 777 #undef __FUNC__ 778 #define __FUNC__ "MatMultTransAdd" 779 /*@ 780 MatMultTransAdd - Computes v3 = v2 + A' * v1. 781 782 Input Parameters: 783 . mat - the matrix 784 . v1, v2 - the vectors 785 786 Output Parameters: 787 . v3 - the result 788 789 Notes: 790 The vectors v1 and v3 cannot be the same. I.e., one cannot 791 call MatMultTransAdd(A,v1,v2,v1). 792 793 .keywords: matrix, multiply, matrix-vector product, transpose, add 794 795 .seealso: MatMultTrans(), MatMultAdd(), MatMult() 796 @*/ 797 int MatMultTransAdd(Mat mat,Vec v1,Vec v2,Vec v3) 798 { 799 int ierr; 800 801 PetscFunctionBegin; 802 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v1,VEC_COOKIE); 803 PetscValidHeaderSpecific(v2,VEC_COOKIE);PetscValidHeaderSpecific(v3,VEC_COOKIE); 804 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 805 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 806 if (!mat->ops.multtransadd) SETERRQ(PETSC_ERR_SUP,0,""); 807 if (v1 == v3) SETERRQ(PETSC_ERR_ARG_IDN,0,"v1 and v3 must be different vectors"); 808 if (mat->M != v1->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v1: global dim"); 809 if (mat->N != v2->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v2: global dim"); 810 if (mat->N != v3->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec v3: global dim"); 811 812 PLogEventBegin(MAT_MultTransAdd,mat,v1,v2,v3); 813 ierr = (*mat->ops.multtransadd)(mat,v1,v2,v3); CHKERRQ(ierr); 814 PLogEventEnd(MAT_MultTransAdd,mat,v1,v2,v3); 815 PetscFunctionReturn(0); 816 } 817 /* ------------------------------------------------------------*/ 818 #undef __FUNC__ 819 #define __FUNC__ "MatGetInfo" 820 /*@C 821 MatGetInfo - Returns information about matrix storage (number of 822 nonzeros, memory, etc.). 823 824 Input Parameters: 825 . mat - the matrix 826 827 Output Parameters: 828 . flag - flag indicating the type of parameters to be returned 829 $ flag = MAT_LOCAL: local matrix 830 $ flag = MAT_GLOBAL_MAX: maximum over all processors 831 $ flag = MAT_GLOBAL_SUM: sum over all processors 832 . info - matrix information context 833 834 Notes: 835 The MatInfo context contains a variety of matrix data, including 836 number of nonzeros allocated and used, number of mallocs during 837 matrix assembly, etc. Additional information for factored matrices 838 is provided (such as the fill ratio, number of mallocs during 839 factorization, etc.). Much of this info is printed to STDOUT 840 when using the runtime options 841 $ -log_info -mat_view_info 842 843 Example for C/C++ Users: 844 See the file $(PETSC_DIR)/include/mat.h for a complete list of 845 data within the MatInfo context. For example, 846 $ 847 $ MatInfo *info; 848 $ Mat A; 849 $ double mal, nz_a, nz_u; 850 $ 851 $ MatGetInfo(A,MAT_LOCAL,&info); 852 $ mal = info->mallocs; 853 $ nz_a = info->nz_allocated; 854 $ 855 856 Example for Fortran Users: 857 Fortran users should declare info as a double precision 858 array of dimension MAT_INFO_SIZE, and then extract the parameters 859 of interest. See the file $(PETSC_DIR)/include/FINCLUDE/mat.h 860 a complete list of parameter names. 861 $ 862 $ double precision info(MAT_INFO_SIZE) 863 $ double precision mal, nz_a 864 $ Mat A 865 $ integer ierr 866 $ 867 $ call MatGetInfo(A,MAT_LOCAL,info,ierr) 868 $ mal = info(MAT_INFO_MALLOCS) 869 $ nz_a = info(MAT_INFO_NZ_ALLOCATED) 870 $ 871 872 .keywords: matrix, get, info, storage, nonzeros, memory, fill 873 @*/ 874 int MatGetInfo(Mat mat,MatInfoType flag,MatInfo *info) 875 { 876 int ierr; 877 878 PetscFunctionBegin; 879 PetscValidHeaderSpecific(mat,MAT_COOKIE); 880 PetscValidPointer(info); 881 if (!mat->ops.getinfo) SETERRQ(PETSC_ERR_SUP,0,""); 882 ierr = (*mat->ops.getinfo)(mat,flag,info);CHKERRQ(ierr); 883 PetscFunctionReturn(0); 884 } 885 886 /* ----------------------------------------------------------*/ 887 #undef __FUNC__ 888 #define __FUNC__ "MatILUDTFactor" 889 /*@ 890 MatILUDTFactor - Performs a drop tolerance ILU factorization. 891 892 Input Parameters: 893 . mat - the matrix 894 . dt - the drop tolerance 895 . maxnz - the maximum number of nonzeros per row allowed? 896 . row - row permutation 897 . col - column permutation 898 899 Output Parameters: 900 . fact - the factored matrix 901 902 .keywords: matrix, factor, LU, in-place 903 904 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 905 @*/ 906 int MatILUDTFactor(Mat mat,double dt,int maxnz,IS row,IS col,Mat *fact) 907 { 908 int ierr; 909 910 PetscFunctionBegin; 911 PetscValidHeaderSpecific(mat,MAT_COOKIE); 912 PetscValidPointer(fact); 913 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 914 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 915 if (!mat->ops.iludtfactor) SETERRQ(PETSC_ERR_SUP,0,""); 916 917 PLogEventBegin(MAT_ILUFactor,mat,row,col,0); 918 ierr = (*mat->ops.iludtfactor)(mat,dt,maxnz,row,col,fact); CHKERRQ(ierr); 919 PLogEventEnd(MAT_ILUFactor,mat,row,col,0); 920 921 PetscFunctionReturn(0); 922 } 923 924 #undef __FUNC__ 925 #define __FUNC__ "MatLUFactor" 926 /*@ 927 MatLUFactor - Performs in-place LU factorization of matrix. 928 929 Input Parameters: 930 . mat - the matrix 931 . row - row permutation 932 . col - column permutation 933 . f - expected fill as ratio of original fill. 934 935 .keywords: matrix, factor, LU, in-place 936 937 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 938 @*/ 939 int MatLUFactor(Mat mat,IS row,IS col,double f) 940 { 941 int ierr; 942 943 PetscFunctionBegin; 944 PetscValidHeaderSpecific(mat,MAT_COOKIE); 945 if (mat->M != mat->N) SETERRQ(1,0,"matrix must be square"); 946 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 947 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 948 if (!mat->ops.lufactor) SETERRQ(PETSC_ERR_SUP,0,""); 949 950 PLogEventBegin(MAT_LUFactor,mat,row,col,0); 951 ierr = (*mat->ops.lufactor)(mat,row,col,f); CHKERRQ(ierr); 952 PLogEventEnd(MAT_LUFactor,mat,row,col,0); 953 PetscFunctionReturn(0); 954 } 955 956 #undef __FUNC__ 957 #define __FUNC__ "MatILUFactor" 958 /*@ 959 MatILUFactor - Performs in-place ILU factorization of matrix. 960 961 Input Parameters: 962 . mat - the matrix 963 . row - row permutation 964 . col - column permutation 965 . f - expected fill as ratio of original fill. 966 . level - number of levels of fill. 967 968 Note: probably really only in-place when level is zero. 969 .keywords: matrix, factor, ILU, in-place 970 971 .seealso: MatILUFactorSymbolic(), MatLUFactorNumeric(), MatCholeskyFactor() 972 @*/ 973 int MatILUFactor(Mat mat,IS row,IS col,double f,int level) 974 { 975 int ierr; 976 977 PetscFunctionBegin; 978 PetscValidHeaderSpecific(mat,MAT_COOKIE); 979 if (mat->M != mat->N) SETERRQ(1,0,"matrix must be square"); 980 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 981 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 982 if (!mat->ops.ilufactor) SETERRQ(PETSC_ERR_SUP,0,""); 983 984 PLogEventBegin(MAT_ILUFactor,mat,row,col,0); 985 ierr = (*mat->ops.ilufactor)(mat,row,col,f,level); CHKERRQ(ierr); 986 PLogEventEnd(MAT_ILUFactor,mat,row,col,0); 987 PetscFunctionReturn(0); 988 } 989 990 #undef __FUNC__ 991 #define __FUNC__ "MatLUFactorSymbolic" 992 /*@ 993 MatLUFactorSymbolic - Performs symbolic LU factorization of matrix. 994 Call this routine before calling MatLUFactorNumeric(). 995 996 Input Parameters: 997 . mat - the matrix 998 . row, col - row and column permutations 999 . f - expected fill as ratio of the original number of nonzeros, 1000 for example 3.0; choosing this parameter well can result in 1001 more efficient use of time and space. Run with the option -log_info 1002 to determine an optimal value to use 1003 1004 Output Parameter: 1005 . fact - new matrix that has been symbolically factored 1006 1007 Notes: 1008 See the file $(PETSC_DIR)/Performace for additional information about 1009 choosing the fill factor for better efficiency. 1010 1011 .keywords: matrix, factor, LU, symbolic, fill 1012 1013 .seealso: MatLUFactor(), MatLUFactorNumeric(), MatCholeskyFactor() 1014 @*/ 1015 int MatLUFactorSymbolic(Mat mat,IS row,IS col,double f,Mat *fact) 1016 { 1017 int ierr; 1018 1019 PetscFunctionBegin; 1020 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1021 if (mat->M != mat->N) SETERRQ(1,0,"matrix must be square"); 1022 PetscValidPointer(fact); 1023 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1024 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1025 if (!mat->ops.lufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,""); 1026 1027 PLogEventBegin(MAT_LUFactorSymbolic,mat,row,col,0); 1028 ierr = (*mat->ops.lufactorsymbolic)(mat,row,col,f,fact); CHKERRQ(ierr); 1029 PLogEventEnd(MAT_LUFactorSymbolic,mat,row,col,0); 1030 PetscFunctionReturn(0); 1031 } 1032 1033 #undef __FUNC__ 1034 #define __FUNC__ "MatLUFactorNumeric" 1035 /*@ 1036 MatLUFactorNumeric - Performs numeric LU factorization of a matrix. 1037 Call this routine after first calling MatLUFactorSymbolic(). 1038 1039 Input Parameters: 1040 . mat - the matrix 1041 . row, col - row and column permutations 1042 1043 Output Parameters: 1044 . fact - symbolically factored matrix that must have been generated 1045 by MatLUFactorSymbolic() 1046 1047 Notes: 1048 See MatLUFactor() for in-place factorization. See 1049 MatCholeskyFactorNumeric() for the symmetric, positive definite case. 1050 1051 .keywords: matrix, factor, LU, numeric 1052 1053 .seealso: MatLUFactorSymbolic(), MatLUFactor(), MatCholeskyFactor() 1054 @*/ 1055 int MatLUFactorNumeric(Mat mat,Mat *fact) 1056 { 1057 int ierr,flg; 1058 1059 PetscFunctionBegin; 1060 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1061 PetscValidPointer(fact); 1062 PetscValidHeaderSpecific(*fact,MAT_COOKIE); 1063 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1064 if (mat->M != (*fact)->M || mat->N != (*fact)->N) 1065 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dim"); 1066 if (!(*fact)->ops.lufactornumeric) SETERRQ(PETSC_ERR_SUP,0,""); 1067 1068 PLogEventBegin(MAT_LUFactorNumeric,mat,*fact,0,0); 1069 ierr = (*(*fact)->ops.lufactornumeric)(mat,fact); CHKERRQ(ierr); 1070 PLogEventEnd(MAT_LUFactorNumeric,mat,*fact,0,0); 1071 ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg); CHKERRQ(ierr); 1072 if (flg) { 1073 ierr = OptionsHasName(0,"-mat_view_contour",&flg); CHKERRQ(ierr); 1074 if (flg) { 1075 ViewerPushFormat(VIEWER_DRAWX_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr); 1076 } 1077 ierr = MatView(*fact,VIEWER_DRAWX_(mat->comm)); CHKERRQ(ierr); 1078 ierr = ViewerFlush(VIEWER_DRAWX_(mat->comm)); CHKERRQ(ierr); 1079 if (flg) { 1080 ViewerPopFormat(VIEWER_DRAWX_(mat->comm));CHKERRQ(ierr); 1081 } 1082 } 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNC__ 1087 #define __FUNC__ "MatCholeskyFactor" 1088 /*@ 1089 MatCholeskyFactor - Performs in-place Cholesky factorization of a 1090 symmetric matrix. 1091 1092 Input Parameters: 1093 . mat - the matrix 1094 . perm - row and column permutations 1095 . f - expected fill as ratio of original fill 1096 1097 Notes: 1098 See MatLUFactor() for the nonsymmetric case. See also 1099 MatCholeskyFactorSymbolic(), and MatCholeskyFactorNumeric(). 1100 1101 .keywords: matrix, factor, in-place, Cholesky 1102 1103 .seealso: MatLUFactor(), MatCholeskyFactorSymbolic(), MatCholeskyFactorNumeric() 1104 @*/ 1105 int MatCholeskyFactor(Mat mat,IS perm,double f) 1106 { 1107 int ierr; 1108 1109 PetscFunctionBegin; 1110 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1111 if (mat->M != mat->N) SETERRQ(1,0,"matrix must be square"); 1112 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1113 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1114 if (!mat->ops.choleskyfactor) SETERRQ(PETSC_ERR_SUP,0,""); 1115 1116 PLogEventBegin(MAT_CholeskyFactor,mat,perm,0,0); 1117 ierr = (*mat->ops.choleskyfactor)(mat,perm,f); CHKERRQ(ierr); 1118 PLogEventEnd(MAT_CholeskyFactor,mat,perm,0,0); 1119 PetscFunctionReturn(0); 1120 } 1121 1122 #undef __FUNC__ 1123 #define __FUNC__ "MatCholeskyFactorSymbolic" 1124 /*@ 1125 MatCholeskyFactorSymbolic - Performs symbolic Cholesky factorization 1126 of a symmetric matrix. 1127 1128 Input Parameters: 1129 . mat - the matrix 1130 . perm - row and column permutations 1131 . f - expected fill as ratio of original 1132 1133 Output Parameter: 1134 . fact - the factored matrix 1135 1136 Notes: 1137 See MatLUFactorSymbolic() for the nonsymmetric case. See also 1138 MatCholeskyFactor() and MatCholeskyFactorNumeric(). 1139 1140 .keywords: matrix, factor, factorization, symbolic, Cholesky 1141 1142 .seealso: MatLUFactorSymbolic(), MatCholeskyFactor(), MatCholeskyFactorNumeric() 1143 @*/ 1144 int MatCholeskyFactorSymbolic(Mat mat,IS perm,double f,Mat *fact) 1145 { 1146 int ierr; 1147 1148 PetscFunctionBegin; 1149 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1150 PetscValidPointer(fact); 1151 if (mat->M != mat->N) SETERRQ(1,0,"matrix must be square"); 1152 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1153 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1154 if (!mat->ops.choleskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,""); 1155 1156 PLogEventBegin(MAT_CholeskyFactorSymbolic,mat,perm,0,0); 1157 ierr = (*mat->ops.choleskyfactorsymbolic)(mat,perm,f,fact); CHKERRQ(ierr); 1158 PLogEventEnd(MAT_CholeskyFactorSymbolic,mat,perm,0,0); 1159 PetscFunctionReturn(0); 1160 } 1161 1162 #undef __FUNC__ 1163 #define __FUNC__ "MatCholeskyFactorNumeric" 1164 /*@ 1165 MatCholeskyFactorNumeric - Performs numeric Cholesky factorization 1166 of a symmetric matrix. Call this routine after first calling 1167 MatCholeskyFactorSymbolic(). 1168 1169 Input Parameter: 1170 . mat - the initial matrix 1171 1172 Output Parameter: 1173 . fact - the factored matrix 1174 1175 .keywords: matrix, factor, numeric, Cholesky 1176 1177 .seealso: MatCholeskyFactorSymbolic(), MatCholeskyFactor(), MatLUFactorNumeric() 1178 @*/ 1179 int MatCholeskyFactorNumeric(Mat mat,Mat *fact) 1180 { 1181 int ierr; 1182 1183 PetscFunctionBegin; 1184 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1185 PetscValidPointer(fact); 1186 if (!mat->ops.choleskyfactornumeric) SETERRQ(PETSC_ERR_SUP,0,""); 1187 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1188 if (mat->M != (*fact)->M || mat->N != (*fact)->N) 1189 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Mat *fact: global dim"); 1190 1191 PLogEventBegin(MAT_CholeskyFactorNumeric,mat,*fact,0,0); 1192 ierr = (*mat->ops.choleskyfactornumeric)(mat,fact); CHKERRQ(ierr); 1193 PLogEventEnd(MAT_CholeskyFactorNumeric,mat,*fact,0,0); 1194 PetscFunctionReturn(0); 1195 } 1196 1197 /* ----------------------------------------------------------------*/ 1198 #undef __FUNC__ 1199 #define __FUNC__ "MatSolve" 1200 /*@ 1201 MatSolve - Solves A x = b, given a factored matrix. 1202 1203 Input Parameters: 1204 . mat - the factored matrix 1205 . b - the right-hand-side vector 1206 1207 Output Parameter: 1208 . x - the result vector 1209 1210 Notes: 1211 The vectors b and x cannot be the same. I.e., one cannot 1212 call MatSolve(A,x,x). 1213 1214 .keywords: matrix, linear system, solve, LU, Cholesky, triangular solve 1215 1216 .seealso: MatSolveAdd(), MatSolveTrans(), MatSolveTransAdd() 1217 @*/ 1218 int MatSolve(Mat mat,Vec b,Vec x) 1219 { 1220 int ierr; 1221 1222 PetscFunctionBegin; 1223 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1224 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1225 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1226 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1227 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1228 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1229 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1230 if (mat->M == 0 && mat->N == 0) PetscFunctionReturn(0); 1231 1232 if (!mat->ops.solve) SETERRQ(PETSC_ERR_SUP,0,""); 1233 PLogEventBegin(MAT_Solve,mat,b,x,0); 1234 ierr = (*mat->ops.solve)(mat,b,x); CHKERRQ(ierr); 1235 PLogEventEnd(MAT_Solve,mat,b,x,0); 1236 PetscFunctionReturn(0); 1237 } 1238 1239 #undef __FUNC__ 1240 #define __FUNC__ "MatForwardSolve" 1241 /* @ 1242 MatForwardSolve - Solves L x = b, given a factored matrix, A = LU. 1243 1244 Input Parameters: 1245 . mat - the factored matrix 1246 . b - the right-hand-side vector 1247 1248 Output Parameter: 1249 . x - the result vector 1250 1251 Notes: 1252 MatSolve() should be used for most applications, as it performs 1253 a forward solve followed by a backward solve. 1254 1255 The vectors b and x cannot be the same. I.e., one cannot 1256 call MatForwardSolve(A,x,x). 1257 1258 .keywords: matrix, forward, LU, Cholesky, triangular solve 1259 1260 .seealso: MatSolve(), MatBackwardSolve() 1261 @ */ 1262 int MatForwardSolve(Mat mat,Vec b,Vec x) 1263 { 1264 int ierr; 1265 1266 PetscFunctionBegin; 1267 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1268 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1269 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1270 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1271 if (!mat->ops.forwardsolve) SETERRQ(PETSC_ERR_SUP,0,""); 1272 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1273 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1274 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1275 1276 PLogEventBegin(MAT_ForwardSolve,mat,b,x,0); 1277 ierr = (*mat->ops.forwardsolve)(mat,b,x); CHKERRQ(ierr); 1278 PLogEventEnd(MAT_ForwardSolve,mat,b,x,0); 1279 PetscFunctionReturn(0); 1280 } 1281 1282 #undef __FUNC__ 1283 #define __FUNC__ "MatBackwardSolve" 1284 /* @ 1285 MatBackwardSolve - Solves U x = b, given a factored matrix, A = LU. 1286 1287 Input Parameters: 1288 . mat - the factored matrix 1289 . b - the right-hand-side vector 1290 1291 Output Parameter: 1292 . x - the result vector 1293 1294 Notes: 1295 MatSolve() should be used for most applications, as it performs 1296 a forward solve followed by a backward solve. 1297 1298 The vectors b and x cannot be the same. I.e., one cannot 1299 call MatBackwardSolve(A,x,x). 1300 1301 .keywords: matrix, backward, LU, Cholesky, triangular solve 1302 1303 .seealso: MatSolve(), MatForwardSolve() 1304 @ */ 1305 int MatBackwardSolve(Mat mat,Vec b,Vec x) 1306 { 1307 int ierr; 1308 1309 PetscFunctionBegin; 1310 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1311 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1312 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1313 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1314 if (!mat->ops.backwardsolve) SETERRQ(PETSC_ERR_SUP,0,""); 1315 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1316 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1317 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1318 1319 PLogEventBegin(MAT_BackwardSolve,mat,b,x,0); 1320 ierr = (*mat->ops.backwardsolve)(mat,b,x); CHKERRQ(ierr); 1321 PLogEventEnd(MAT_BackwardSolve,mat,b,x,0); 1322 PetscFunctionReturn(0); 1323 } 1324 1325 #undef __FUNC__ 1326 #define __FUNC__ "MatSolveAdd" 1327 /*@ 1328 MatSolveAdd - Computes x = y + inv(A)*b, given a factored matrix. 1329 1330 Input Parameters: 1331 . mat - the factored matrix 1332 . b - the right-hand-side vector 1333 . y - the vector to be added to 1334 1335 Output Parameter: 1336 . x - the result vector 1337 1338 Notes: 1339 The vectors b and x cannot be the same. I.e., one cannot 1340 call MatSolveAdd(A,x,y,x). 1341 1342 .keywords: matrix, linear system, solve, LU, Cholesky, add 1343 1344 .seealso: MatSolve(), MatSolveTrans(), MatSolveTransAdd() 1345 @*/ 1346 int MatSolveAdd(Mat mat,Vec b,Vec y,Vec x) 1347 { 1348 Scalar one = 1.0; 1349 Vec tmp; 1350 int ierr; 1351 1352 PetscFunctionBegin; 1353 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1354 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1355 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1356 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1357 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1358 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1359 if (mat->M != y->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim"); 1360 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1361 if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim"); 1362 1363 PLogEventBegin(MAT_SolveAdd,mat,b,x,y); 1364 if (mat->ops.solveadd) { 1365 ierr = (*mat->ops.solveadd)(mat,b,y,x); CHKERRQ(ierr); 1366 } 1367 else { 1368 /* do the solve then the add manually */ 1369 if (x != y) { 1370 ierr = MatSolve(mat,b,x); CHKERRQ(ierr); 1371 ierr = VecAXPY(&one,y,x); CHKERRQ(ierr); 1372 } else { 1373 ierr = VecDuplicate(x,&tmp); CHKERRQ(ierr); 1374 PLogObjectParent(mat,tmp); 1375 ierr = VecCopy(x,tmp); CHKERRQ(ierr); 1376 ierr = MatSolve(mat,b,x); CHKERRQ(ierr); 1377 ierr = VecAXPY(&one,tmp,x); CHKERRQ(ierr); 1378 ierr = VecDestroy(tmp); CHKERRQ(ierr); 1379 } 1380 } 1381 PLogEventEnd(MAT_SolveAdd,mat,b,x,y); 1382 PetscFunctionReturn(0); 1383 } 1384 1385 #undef __FUNC__ 1386 #define __FUNC__ "MatSolveTrans" 1387 /*@ 1388 MatSolveTrans - Solves A' x = b, given a factored matrix. 1389 1390 Input Parameters: 1391 . mat - the factored matrix 1392 . b - the right-hand-side vector 1393 1394 Output Parameter: 1395 . x - the result vector 1396 1397 Notes: 1398 The vectors b and x cannot be the same. I.e., one cannot 1399 call MatSolveTrans(A,x,x). 1400 1401 .keywords: matrix, linear system, solve, LU, Cholesky, transpose 1402 1403 .seealso: MatSolve(), MatSolveAdd(), MatSolveTransAdd() 1404 @*/ 1405 int MatSolveTrans(Mat mat,Vec b,Vec x) 1406 { 1407 int ierr; 1408 1409 PetscFunctionBegin; 1410 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1411 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1412 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1413 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1414 if (!mat->ops.solvetrans) SETERRQ(PETSC_ERR_SUP,0,""); 1415 if (mat->M != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1416 if (mat->N != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1417 1418 PLogEventBegin(MAT_SolveTrans,mat,b,x,0); 1419 ierr = (*mat->ops.solvetrans)(mat,b,x); CHKERRQ(ierr); 1420 PLogEventEnd(MAT_SolveTrans,mat,b,x,0); 1421 PetscFunctionReturn(0); 1422 } 1423 1424 #undef __FUNC__ 1425 #define __FUNC__ "MatSolveTransAdd" 1426 /*@ 1427 MatSolveTransAdd - Computes x = y + inv(trans(A)) b, given a 1428 factored matrix. 1429 1430 Input Parameters: 1431 . mat - the factored matrix 1432 . b - the right-hand-side vector 1433 . y - the vector to be added to 1434 1435 Output Parameter: 1436 . x - the result vector 1437 1438 Notes: 1439 The vectors b and x cannot be the same. I.e., one cannot 1440 call MatSolveTransAdd(A,x,y,x). 1441 1442 .keywords: matrix, linear system, solve, LU, Cholesky, transpose, add 1443 1444 .seealso: MatSolve(), MatSolveAdd(), MatSolveTrans() 1445 @*/ 1446 int MatSolveTransAdd(Mat mat,Vec b,Vec y,Vec x) 1447 { 1448 Scalar one = 1.0; 1449 int ierr; 1450 Vec tmp; 1451 1452 PetscFunctionBegin; 1453 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(y,VEC_COOKIE); 1454 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1455 if (x == b) SETERRQ(PETSC_ERR_ARG_IDN,0,"x and b must be different vectors"); 1456 if (!mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Unfactored matrix"); 1457 if (mat->M != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1458 if (mat->N != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1459 if (mat->N != y->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec y: global dim"); 1460 if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Vec x,Vec y: local dim"); 1461 1462 PLogEventBegin(MAT_SolveTransAdd,mat,b,x,y); 1463 if (mat->ops.solvetransadd) { 1464 ierr = (*mat->ops.solvetransadd)(mat,b,y,x); CHKERRQ(ierr); 1465 } 1466 else { 1467 /* do the solve then the add manually */ 1468 if (x != y) { 1469 ierr = MatSolveTrans(mat,b,x); CHKERRQ(ierr); 1470 ierr = VecAXPY(&one,y,x); CHKERRQ(ierr); 1471 } 1472 else { 1473 ierr = VecDuplicate(x,&tmp); CHKERRQ(ierr); 1474 PLogObjectParent(mat,tmp); 1475 ierr = VecCopy(x,tmp); CHKERRQ(ierr); 1476 ierr = MatSolveTrans(mat,b,x); CHKERRQ(ierr); 1477 ierr = VecAXPY(&one,tmp,x); CHKERRQ(ierr); 1478 ierr = VecDestroy(tmp); CHKERRQ(ierr); 1479 } 1480 } 1481 PLogEventEnd(MAT_SolveTransAdd,mat,b,x,y); 1482 PetscFunctionReturn(0); 1483 } 1484 /* ----------------------------------------------------------------*/ 1485 1486 #undef __FUNC__ 1487 #define __FUNC__ "MatRelax" 1488 /*@ 1489 MatRelax - Computes one relaxation sweep. 1490 1491 Input Parameters: 1492 . mat - the matrix 1493 . b - the right hand side 1494 . omega - the relaxation factor 1495 . flag - flag indicating the type of SOR, one of 1496 $ SOR_FORWARD_SWEEP 1497 $ SOR_BACKWARD_SWEEP 1498 $ SOR_SYMMETRIC_SWEEP (SSOR method) 1499 $ SOR_LOCAL_FORWARD_SWEEP 1500 $ SOR_LOCAL_BACKWARD_SWEEP 1501 $ SOR_LOCAL_SYMMETRIC_SWEEP (local SSOR) 1502 $ SOR_APPLY_UPPER, SOR_APPLY_LOWER - applies 1503 $ upper/lower triangular part of matrix to 1504 $ vector (with omega) 1505 $ SOR_ZERO_INITIAL_GUESS - zero initial guess 1506 . shift - diagonal shift 1507 . its - the number of iterations 1508 1509 Output Parameters: 1510 . x - the solution (can contain an initial guess) 1511 1512 Notes: 1513 SOR_LOCAL_FORWARD_SWEEP, SOR_LOCAL_BACKWARD_SWEEP, and 1514 SOR_LOCAL_SYMMETRIC_SWEEP perform seperate independent smoothings 1515 on each processor. 1516 1517 Application programmers will not generally use MatRelax() directly, 1518 but instead will employ the SLES/PC interface. 1519 1520 Notes for Advanced Users: 1521 The flags are implemented as bitwise inclusive or operations. 1522 For example, use (SOR_ZERO_INITIAL_GUESS | SOR_SYMMETRIC_SWEEP) 1523 to specify a zero initial guess for SSOR. 1524 1525 .keywords: matrix, relax, relaxation, sweep 1526 @*/ 1527 int MatRelax(Mat mat,Vec b,double omega,MatSORType flag,double shift, 1528 int its,Vec x) 1529 { 1530 int ierr; 1531 1532 PetscFunctionBegin; 1533 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1534 PetscValidHeaderSpecific(b,VEC_COOKIE); PetscValidHeaderSpecific(x,VEC_COOKIE); 1535 if (!mat->ops.relax) SETERRQ(PETSC_ERR_SUP,0,""); 1536 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1537 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1538 if (mat->N != x->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec x: global dim"); 1539 if (mat->M != b->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: global dim"); 1540 if (mat->m != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat mat,Vec b: local dim"); 1541 1542 PLogEventBegin(MAT_Relax,mat,b,x,0); 1543 ierr =(*mat->ops.relax)(mat,b,omega,flag,shift,its,x); CHKERRQ(ierr); 1544 PLogEventEnd(MAT_Relax,mat,b,x,0); 1545 PetscFunctionReturn(0); 1546 } 1547 1548 #undef __FUNC__ 1549 #define __FUNC__ "MatCopy_Basic" 1550 /* 1551 Default matrix copy routine. 1552 */ 1553 int MatCopy_Basic(Mat A,Mat B) 1554 { 1555 int ierr,i,rstart,rend,nz,*cwork; 1556 Scalar *vwork; 1557 1558 PetscFunctionBegin; 1559 ierr = MatZeroEntries(B); CHKERRQ(ierr); 1560 ierr = MatGetOwnershipRange(A,&rstart,&rend); CHKERRQ(ierr); 1561 for (i=rstart; i<rend; i++) { 1562 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1563 ierr = MatSetValues(B,1,&i,nz,cwork,vwork,INSERT_VALUES); CHKERRQ(ierr); 1564 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1565 } 1566 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1567 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1568 PetscFunctionReturn(0); 1569 } 1570 1571 #undef __FUNC__ 1572 #define __FUNC__ "MatCopy" 1573 /*@C 1574 MatCopy - Copys a matrix to another matrix. 1575 1576 Input Parameters: 1577 . A - the matrix 1578 1579 Output Parameter: 1580 . B - where the copy is put 1581 1582 Notes: 1583 MatCopy() copies the matrix entries of a matrix to another existing 1584 matrix (after first zeroing the second matrix). A related routine is 1585 MatConvert(), which first creates a new matrix and then copies the data. 1586 1587 .keywords: matrix, copy, convert 1588 1589 .seealso: MatConvert() 1590 @*/ 1591 int MatCopy(Mat A,Mat B) 1592 { 1593 int ierr; 1594 1595 PetscFunctionBegin; 1596 PetscValidHeaderSpecific(A,MAT_COOKIE); PetscValidHeaderSpecific(B,MAT_COOKIE); 1597 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1598 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1599 if (A->M != B->M || A->N != B->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim"); 1600 1601 PLogEventBegin(MAT_Copy,A,B,0,0); 1602 if (A->ops.copy) { 1603 ierr = (*A->ops.copy)(A,B); CHKERRQ(ierr); 1604 } 1605 else { /* generic conversion */ 1606 ierr = MatCopy_Basic(A,B); CHKERRQ(ierr); 1607 } 1608 PLogEventEnd(MAT_Copy,A,B,0,0); 1609 PetscFunctionReturn(0); 1610 } 1611 1612 static int MatConvertersSet = 0; 1613 static int (*MatConverters[MAX_MATRIX_TYPES][MAX_MATRIX_TYPES])(Mat,MatType,Mat*) = 1614 {{0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1615 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1616 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1617 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1618 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}, 1619 {0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0}}; 1620 1621 #undef __FUNC__ 1622 #define __FUNC__ "MatConvertRegister" 1623 /*@C 1624 MatConvertRegister - Allows one to register a routine that converts between 1625 two matrix types. 1626 1627 Input Parameters: 1628 . intype - the type of matrix (defined in include/mat.h), for example, MATSEQAIJ. 1629 . outtype - new matrix type, or MATSAME 1630 1631 .seealso: MatConvertRegisterAll() 1632 1633 @*/ 1634 int MatConvertRegister(MatType intype,MatType outtype,int (*converter)(Mat,MatType,Mat*)) 1635 { 1636 PetscFunctionBegin; 1637 MatConverters[intype][outtype] = converter; 1638 MatConvertersSet = 1; 1639 PetscFunctionReturn(0); 1640 } 1641 1642 #undef __FUNC__ 1643 #define __FUNC__ "MatConvert" 1644 /*@C 1645 MatConvert - Converts a matrix to another matrix, either of the same 1646 or different type. 1647 1648 Input Parameters: 1649 . mat - the matrix 1650 . newtype - new matrix type. Use MATSAME to create a new matrix of the 1651 same type as the original matrix. 1652 1653 Output Parameter: 1654 . M - pointer to place new matrix 1655 1656 Notes: 1657 MatConvert() first creates a new matrix and then copies the data from 1658 the first matrix. A related routine is MatCopy(), which copies the matrix 1659 entries of one matrix to another already existing matrix context. 1660 1661 .keywords: matrix, copy, convert 1662 1663 .seealso: MatCopy(), MatDuplicate() 1664 @*/ 1665 int MatConvert(Mat mat,MatType newtype,Mat *M) 1666 { 1667 int ierr; 1668 1669 PetscFunctionBegin; 1670 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1671 PetscValidPointer(M); 1672 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1673 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1674 1675 if (newtype > MAX_MATRIX_TYPES || newtype < -1) { 1676 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Not a valid matrix type"); 1677 } 1678 *M = 0; 1679 1680 if (!MatConvertersSet) { 1681 ierr = MatLoadRegisterAll(); CHKERRQ(ierr); 1682 } 1683 1684 PLogEventBegin(MAT_Convert,mat,0,0,0); 1685 if ((newtype == mat->type || newtype == MATSAME) && mat->ops.convertsametype) { 1686 ierr = (*mat->ops.convertsametype)(mat,M,COPY_VALUES); CHKERRQ(ierr); 1687 } else { 1688 if (!MatConvertersSet) { 1689 ierr = MatConvertRegisterAll(); CHKERRQ(ierr); 1690 } 1691 if (!MatConverters[mat->type][newtype]) { 1692 SETERRQ(PETSC_ERR_ARG_WRONG,1,"Invalid matrix type, or matrix converter not registered"); 1693 } 1694 ierr = (*MatConverters[mat->type][newtype])(mat,newtype,M); CHKERRQ(ierr); 1695 } 1696 PLogEventEnd(MAT_Convert,mat,0,0,0); 1697 PetscFunctionReturn(0); 1698 } 1699 1700 #undef __FUNC__ 1701 #define __FUNC__ "MatDuplicate" 1702 /*@C 1703 MatDuplicate - Duplicates a matrix including the non-zero structure, but 1704 does not copy over the values. 1705 1706 Input Parameters: 1707 . mat - the matrix 1708 1709 Output Parameter: 1710 . M - pointer to place new matrix 1711 1712 .keywords: matrix, copy, convert, duplicate 1713 1714 .seealso: MatCopy(), MatDuplicate(), MatConvert() 1715 @*/ 1716 int MatDuplicate(Mat mat,Mat *M) 1717 { 1718 int ierr; 1719 1720 PetscFunctionBegin; 1721 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1722 PetscValidPointer(M); 1723 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1724 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1725 1726 *M = 0; 1727 PLogEventBegin(MAT_Convert,mat,0,0,0); 1728 if (!mat->ops.convertsametype) { 1729 SETERRQ(1,1,"Not written for this matrix type"); 1730 } 1731 ierr = (*mat->ops.convertsametype)(mat,M,DO_NOT_COPY_VALUES); CHKERRQ(ierr); 1732 PLogEventEnd(MAT_Convert,mat,0,0,0); 1733 PetscFunctionReturn(0); 1734 } 1735 1736 #undef __FUNC__ 1737 #define __FUNC__ "MatGetDiagonal" 1738 /*@ 1739 MatGetDiagonal - Gets the diagonal of a matrix. 1740 1741 Input Parameters: 1742 . mat - the matrix 1743 . v - the vector for storing the diagonal 1744 1745 Output Parameter: 1746 . v - the diagonal of the matrix 1747 1748 Notes: For the SeqAIJ matrix format, this routine may also be called 1749 on a LU factored matrix; in that case it routines the reciprocal of 1750 the diagonal entries in U. It returns the entries permuted by the 1751 row and column permutation used during the symbolic factorization. 1752 1753 .keywords: matrix, get, diagonal 1754 @*/ 1755 int MatGetDiagonal(Mat mat,Vec v) 1756 { 1757 int ierr; 1758 1759 PetscFunctionBegin; 1760 PetscValidHeaderSpecific(mat,MAT_COOKIE);PetscValidHeaderSpecific(v,VEC_COOKIE); 1761 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1762 /* 1763 The error checking for a factored matrix is handled inside 1764 each matrix type, since MatGetDiagonal() is supported by 1765 factored AIJ matrices 1766 */ 1767 /* if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); */ 1768 if (!mat->ops.getdiagonal) SETERRQ(PETSC_ERR_SUP,0,""); 1769 ierr = (*mat->ops.getdiagonal)(mat,v);CHKERRQ(ierr); 1770 PetscFunctionReturn(0); 1771 } 1772 1773 #undef __FUNC__ 1774 #define __FUNC__ "MatTranspose" 1775 /*@C 1776 MatTranspose - Computes an in-place or out-of-place transpose of a matrix. 1777 1778 Input Parameter: 1779 . mat - the matrix to transpose 1780 1781 Output Parameters: 1782 . B - the transpose (or pass in PETSC_NULL for an in-place transpose) 1783 1784 .keywords: matrix, transpose 1785 1786 .seealso: MatMultTrans(), MatMultTransAdd() 1787 @*/ 1788 int MatTranspose(Mat mat,Mat *B) 1789 { 1790 int ierr; 1791 1792 PetscFunctionBegin; 1793 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1794 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1795 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1796 if (!mat->ops.transpose) SETERRQ(PETSC_ERR_SUP,0,""); 1797 ierr = (*mat->ops.transpose)(mat,B);CHKERRQ(ierr); 1798 PetscFunctionReturn(0); 1799 } 1800 1801 #undef __FUNC__ 1802 #define __FUNC__ "MatPermute" 1803 /*@C 1804 MatPermute - Creates a new matrix with rows and columns permuted from the 1805 original. 1806 1807 Input Parameter: 1808 . mat - the matrix to permute 1809 . row - row permutation 1810 . col - column permutation 1811 1812 Output Parameters: 1813 . B - the permuted matrix 1814 1815 .keywords: matrix, transpose 1816 1817 .seealso: MatGetReordering() 1818 @*/ 1819 int MatPermute(Mat mat,IS row,IS col,Mat *B) 1820 { 1821 int ierr; 1822 1823 PetscFunctionBegin; 1824 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1825 PetscValidHeaderSpecific(row,IS_COOKIE); 1826 PetscValidHeaderSpecific(col,IS_COOKIE); 1827 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1828 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1829 if (!mat->ops.permute) SETERRQ(PETSC_ERR_SUP,0,""); 1830 ierr = (*mat->ops.permute)(mat,row,col,B);CHKERRQ(ierr); 1831 PetscFunctionReturn(0); 1832 } 1833 1834 #undef __FUNC__ 1835 #define __FUNC__ "MatEqual" 1836 /*@ 1837 MatEqual - Compares two matrices. 1838 1839 Input Parameters: 1840 . A - the first matrix 1841 . B - the second matrix 1842 1843 Output Parameter: 1844 . flg : PETSC_TRUE if the matrices are equal; 1845 PETSC_FALSE otherwise. 1846 1847 .keywords: matrix, equal, equivalent 1848 @*/ 1849 int MatEqual(Mat A,Mat B,PetscTruth *flg) 1850 { 1851 int ierr; 1852 1853 PetscFunctionBegin; 1854 PetscValidHeaderSpecific(A,MAT_COOKIE); PetscValidHeaderSpecific(B,MAT_COOKIE); 1855 PetscValidIntPointer(flg); 1856 if (!A->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1857 if (!B->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1858 if (A->M != B->M || A->N != B->N) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Mat A,Mat B: global dim"); 1859 if (!A->ops.equal) SETERRQ(PETSC_ERR_SUP,0,""); 1860 ierr = (*A->ops.equal)(A,B,flg);CHKERRQ(ierr); 1861 PetscFunctionReturn(0); 1862 } 1863 1864 #undef __FUNC__ 1865 #define __FUNC__ "MatDiagonalScale" 1866 /*@ 1867 MatDiagonalScale - Scales a matrix on the left and right by diagonal 1868 matrices that are stored as vectors. Either of the two scaling 1869 matrices can be PETSC_NULL. 1870 1871 Input Parameters: 1872 . mat - the matrix to be scaled 1873 . l - the left scaling vector (or PETSC_NULL) 1874 . r - the right scaling vector (or PETSC_NULL) 1875 1876 Notes: 1877 MatDiagonalScale() computes A <- LAR, where 1878 $ L = a diagonal matrix 1879 $ R = a diagonal matrix 1880 1881 .keywords: matrix, diagonal, scale 1882 1883 .seealso: MatDiagonalScale() 1884 @*/ 1885 int MatDiagonalScale(Mat mat,Vec l,Vec r) 1886 { 1887 int ierr; 1888 1889 PetscFunctionBegin; 1890 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1891 if (!mat->ops.diagonalscale) SETERRQ(PETSC_ERR_SUP,0,""); 1892 if (l) PetscValidHeaderSpecific(l,VEC_COOKIE); 1893 if (r) PetscValidHeaderSpecific(r,VEC_COOKIE); 1894 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1895 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1896 1897 PLogEventBegin(MAT_Scale,mat,0,0,0); 1898 ierr = (*mat->ops.diagonalscale)(mat,l,r); CHKERRQ(ierr); 1899 PLogEventEnd(MAT_Scale,mat,0,0,0); 1900 PetscFunctionReturn(0); 1901 } 1902 1903 #undef __FUNC__ 1904 #define __FUNC__ "MatScale" 1905 /*@ 1906 MatScale - Scales all elements of a matrix by a given number. 1907 1908 Input Parameters: 1909 . mat - the matrix to be scaled 1910 . a - the scaling value 1911 1912 Output Parameter: 1913 . mat - the scaled matrix 1914 1915 .keywords: matrix, scale 1916 1917 .seealso: MatDiagonalScale() 1918 @*/ 1919 int MatScale(Scalar *a,Mat mat) 1920 { 1921 int ierr; 1922 1923 PetscFunctionBegin; 1924 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1925 PetscValidScalarPointer(a); 1926 if (!mat->ops.scale) SETERRQ(PETSC_ERR_SUP,0,""); 1927 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1928 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1929 1930 PLogEventBegin(MAT_Scale,mat,0,0,0); 1931 ierr = (*mat->ops.scale)(a,mat); CHKERRQ(ierr); 1932 PLogEventEnd(MAT_Scale,mat,0,0,0); 1933 PetscFunctionReturn(0); 1934 } 1935 1936 #undef __FUNC__ 1937 #define __FUNC__ "MatNorm" 1938 /*@ 1939 MatNorm - Calculates various norms of a matrix. 1940 1941 Input Parameters: 1942 . mat - the matrix 1943 . type - the type of norm, NORM_1, NORM_2, NORM_FROBENIUS, NORM_INFINITY 1944 1945 Output Parameters: 1946 . norm - the resulting norm 1947 1948 .keywords: matrix, norm, Frobenius 1949 @*/ 1950 int MatNorm(Mat mat,NormType type,double *norm) 1951 { 1952 int ierr; 1953 1954 PetscFunctionBegin; 1955 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1956 PetscValidScalarPointer(norm); 1957 1958 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1959 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1960 if (!mat->ops.norm) SETERRQ(PETSC_ERR_SUP,0,"Not for this matrix type"); 1961 ierr = (*mat->ops.norm)(mat,type,norm);CHKERRQ(ierr); 1962 PetscFunctionReturn(0); 1963 } 1964 1965 /* 1966 This variable is used to prevent counting of MatAssemblyBegin() that 1967 are called from within a MatAssemblyEnd(). 1968 */ 1969 static int MatAssemblyEnd_InUse = 0; 1970 #undef __FUNC__ 1971 #define __FUNC__ "MatAssemblyBegin" 1972 /*@ 1973 MatAssemblyBegin - Begins assembling the matrix. This routine should 1974 be called after completing all calls to MatSetValues(). 1975 1976 Input Parameters: 1977 . mat - the matrix 1978 . type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 1979 1980 Notes: 1981 MatSetValues() generally caches the values. The matrix is ready to 1982 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 1983 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 1984 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 1985 using the matrix. 1986 1987 .keywords: matrix, assembly, assemble, begin 1988 1989 .seealso: MatAssemblyEnd(), MatSetValues() 1990 @*/ 1991 int MatAssemblyBegin(Mat mat,MatAssemblyType type) 1992 { 1993 int ierr; 1994 1995 PetscFunctionBegin; 1996 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1997 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1998 if (mat->assembled) { 1999 mat->was_assembled = PETSC_TRUE; 2000 mat->assembled = PETSC_FALSE; 2001 } 2002 if (!MatAssemblyEnd_InUse) { 2003 PLogEventBegin(MAT_AssemblyBegin,mat,0,0,0); 2004 if (mat->ops.assemblybegin){ierr = (*mat->ops.assemblybegin)(mat,type);CHKERRQ(ierr);} 2005 PLogEventEnd(MAT_AssemblyBegin,mat,0,0,0); 2006 } else { 2007 if (mat->ops.assemblybegin){ierr = (*mat->ops.assemblybegin)(mat,type);CHKERRQ(ierr);} 2008 } 2009 PetscFunctionReturn(0); 2010 } 2011 2012 2013 #undef __FUNC__ 2014 #define __FUNC__ "MatView_Private" 2015 /* 2016 Processes command line options to determine if/how a matrix 2017 is to be viewed. Called by MatAssemblyEnd() and MatLoad(). 2018 */ 2019 int MatView_Private(Mat mat) 2020 { 2021 int ierr,flg; 2022 2023 PetscFunctionBegin; 2024 ierr = OptionsHasName(PETSC_NULL,"-mat_view_info",&flg); CHKERRQ(ierr); 2025 if (flg) { 2026 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO,0);CHKERRQ(ierr); 2027 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2028 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2029 } 2030 ierr = OptionsHasName(PETSC_NULL,"-mat_view_info_detailed",&flg);CHKERRQ(ierr); 2031 if (flg) { 2032 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_INFO_LONG,0);CHKERRQ(ierr); 2033 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2034 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2035 } 2036 ierr = OptionsHasName(PETSC_NULL,"-mat_view",&flg); CHKERRQ(ierr); 2037 if (flg) { 2038 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2039 } 2040 ierr = OptionsHasName(PETSC_NULL,"-mat_view_matlab",&flg); CHKERRQ(ierr); 2041 if (flg) { 2042 ierr = ViewerPushFormat(VIEWER_STDOUT_(mat->comm),VIEWER_FORMAT_ASCII_MATLAB,"M");CHKERRQ(ierr); 2043 ierr = MatView(mat,VIEWER_STDOUT_(mat->comm)); CHKERRQ(ierr); 2044 ierr = ViewerPopFormat(VIEWER_STDOUT_(mat->comm));CHKERRQ(ierr); 2045 } 2046 ierr = OptionsHasName(PETSC_NULL,"-mat_view_draw",&flg); CHKERRQ(ierr); 2047 if (flg) { 2048 ierr = OptionsHasName(0,"-mat_view_contour",&flg); CHKERRQ(ierr); 2049 if (flg) { 2050 ViewerPushFormat(VIEWER_DRAWX_(mat->comm),VIEWER_FORMAT_DRAW_CONTOUR,0);CHKERRQ(ierr); 2051 } 2052 ierr = MatView(mat,VIEWER_DRAWX_(mat->comm)); CHKERRQ(ierr); 2053 ierr = ViewerFlush(VIEWER_DRAWX_(mat->comm)); CHKERRQ(ierr); 2054 if (flg) { 2055 ViewerPopFormat(VIEWER_DRAWX_(mat->comm));CHKERRQ(ierr); 2056 } 2057 } 2058 PetscFunctionReturn(0); 2059 } 2060 2061 #undef __FUNC__ 2062 #define __FUNC__ "MatAssemblyEnd" 2063 /*@ 2064 MatAssemblyEnd - Completes assembling the matrix. This routine should 2065 be called after MatAssemblyBegin(). 2066 2067 Input Parameters: 2068 . mat - the matrix 2069 . type - type of assembly, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY 2070 2071 Options Database Keys: 2072 $ -mat_view_info : Prints info on matrix at 2073 $ conclusion of MatEndAssembly() 2074 $ -mat_view_info_detailed: Prints more detailed info. 2075 $ -mat_view : Prints matrix in ASCII format. 2076 $ -mat_view_matlab : Prints matrix in Matlab format. 2077 $ -mat_view_draw : Draws nonzero structure of matrix, 2078 $ using MatView() and DrawOpenX(). 2079 $ -display <name> : Set display name (default is host) 2080 $ -draw_pause <sec> : Set number of seconds to pause after display 2081 2082 Notes: 2083 MatSetValues() generally caches the values. The matrix is ready to 2084 use only after MatAssemblyBegin() and MatAssemblyEnd() have been called. 2085 Use MAT_FLUSH_ASSEMBLY when switching between ADD_VALUES and INSERT_VALUES 2086 in MatSetValues(); use MAT_FINAL_ASSEMBLY for the final assembly before 2087 using the matrix. 2088 2089 .keywords: matrix, assembly, assemble, end 2090 2091 .seealso: MatAssemblyBegin(), MatSetValues(), DrawOpenX(), MatView() 2092 @*/ 2093 int MatAssemblyEnd(Mat mat,MatAssemblyType type) 2094 { 2095 int ierr; 2096 static int inassm = 0; 2097 2098 PetscFunctionBegin; 2099 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2100 2101 inassm++; 2102 MatAssemblyEnd_InUse++; 2103 PLogEventBegin(MAT_AssemblyEnd,mat,0,0,0); 2104 if (mat->ops.assemblyend) { 2105 ierr = (*mat->ops.assemblyend)(mat,type); CHKERRQ(ierr); 2106 } 2107 2108 /* Flush assembly is not a true assembly */ 2109 if (type != MAT_FLUSH_ASSEMBLY) { 2110 mat->assembled = PETSC_TRUE; mat->num_ass++; 2111 } 2112 mat->insertmode = NOT_SET_VALUES; 2113 PLogEventEnd(MAT_AssemblyEnd,mat,0,0,0); 2114 MatAssemblyEnd_InUse--; 2115 2116 if (inassm == 1 && type != MAT_FLUSH_ASSEMBLY) { 2117 ierr = MatView_Private(mat); CHKERRQ(ierr); 2118 } 2119 inassm--; 2120 PetscFunctionReturn(0); 2121 } 2122 2123 2124 #undef __FUNC__ 2125 #define __FUNC__ "MatCompress" 2126 /*@ 2127 MatCompress - Tries to store the matrix in as little space as 2128 possible. May fail if memory is already fully used, since it 2129 tries to allocate new space. 2130 2131 Input Parameters: 2132 . mat - the matrix 2133 2134 .keywords: matrix, compress 2135 @*/ 2136 int MatCompress(Mat mat) 2137 { 2138 int ierr; 2139 2140 PetscFunctionBegin; 2141 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2142 if (mat->ops.compress) {ierr = (*mat->ops.compress)(mat);CHKERRQ(ierr);} 2143 PetscFunctionReturn(0); 2144 } 2145 2146 #undef __FUNC__ 2147 #define __FUNC__ "MatSetOption" 2148 /*@ 2149 MatSetOption - Sets a parameter option for a matrix. Some options 2150 may be specific to certain storage formats. Some options 2151 determine how values will be inserted (or added). Sorted, 2152 row-oriented input will generally assemble the fastest. The default 2153 is row-oriented, nonsorted input. 2154 2155 Input Parameters: 2156 . mat - the matrix 2157 . option - the option, one of those listed below (and possibly others), 2158 e.g., MAT_ROWS_SORTED, MAT_NEW_NONZERO_LOCATION_ERROR 2159 2160 Options Describing Matrix Structure: 2161 $ MAT_SYMMETRIC - symmetric in terms of both structure and value 2162 $ MAT_STRUCTURALLY_SYMMETRIC - symmetric nonzero structure 2163 2164 Options For Use with MatSetValues(): 2165 Insert a logically dense subblock, which can be 2166 2167 $ MAT_ROW_ORIENTED - row-oriented 2168 $ MAT_COLUMN_ORIENTED - column-oriented 2169 $ MAT_ROWS_SORTED - sorted by row 2170 $ MAT_ROWS_UNSORTED - not sorted by row 2171 $ MAT_COLUMNS_SORTED - sorted by column 2172 $ MAT_COLUMNS_UNSORTED - not sorted by column 2173 2174 When (re)assembling a matrix, we can restrict the input for 2175 efficiency/debugging purposes. 2176 2177 $ MAT_NO_NEW_NONZERO_LOCATIONS - additional insertions will not be 2178 allowed if they generate a new nonzero 2179 $ MAT_YES_NEW_NONZERO_LOCATIONS - additional insertions will be allowed 2180 $ MAT_NO_NEW_DIAGONALS - additional insertions will not be allowed if 2181 they generate a nonzero in a new diagonal (for block diagonal format only) 2182 $ MAT_YES_NEW_DIAGONALS - new diagonals will be allowed (for block diagonal format only) 2183 $ MAT_IGNORE_OFF_PROC_ENTRIES - drop off-processor entries 2184 $ MAT_NEW_NONZERO_LOCATION_ERROR - generate error for new matrix entry 2185 2186 Notes: 2187 Some options are relevant only for particular matrix types and 2188 are thus ignored by others. Other options are not supported by 2189 certain matrix types and will generate an error message if set. 2190 2191 If using a Fortran 77 module to compute a matrix, one may need to 2192 use the column-oriented option (or convert to the row-oriented 2193 format). 2194 2195 MAT_NO_NEW_NONZERO_LOCATIONS indicates that any add or insertion 2196 that would generate a new entry in the nonzero structure is instead 2197 ignored. Thus, if memory has not alredy been allocated for this particular 2198 data, then the insertion is ignored. For dense matrices, in which 2199 the entire array is allocated, no entries are ever ignored. 2200 2201 MAT_NEW_NONZERO_LOCATION_ERROR indicates that any add or insertion 2202 that would generate a new entry in the nonzero structure instead produces 2203 an error. (Currently supported for AIJ and BAIJ formats only.) 2204 This is a useful flag when using SAME_NONZERO_PATTERN in calling 2205 SLESSetOperators() to ensure that the nonzero pattern truely does 2206 remain unchanged. 2207 2208 MAT_NEW_NONZERO_ALLOCATION_ERROR indicates that any add or insertion 2209 that would generate a new entry that has not been preallocated will 2210 instead produce an error. (Currently supported for AIJ and BAIJ formats 2211 only.) This is a useful flag when debugging matrix memory preallocation. 2212 2213 MAT_IGNORE_OFF_PROC_ENTRIES indicates entries destined for 2214 other processors should be dropped, rather than stashed. 2215 This is useful if you know that the "owning" processor is also 2216 always generating the correct matrix entries, so that PETSc need 2217 not transfer duplicate entries generated on another processor. 2218 2219 .keywords: matrix, option, row-oriented, column-oriented, sorted, nonzero 2220 @*/ 2221 int MatSetOption(Mat mat,MatOption op) 2222 { 2223 int ierr; 2224 2225 PetscFunctionBegin; 2226 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2227 if (mat->ops.setoption) {ierr = (*mat->ops.setoption)(mat,op);CHKERRQ(ierr);} 2228 PetscFunctionReturn(0); 2229 } 2230 2231 #undef __FUNC__ 2232 #define __FUNC__ "MatZeroEntries" 2233 /*@ 2234 MatZeroEntries - Zeros all entries of a matrix. For sparse matrices 2235 this routine retains the old nonzero structure. 2236 2237 Input Parameters: 2238 . mat - the matrix 2239 2240 .keywords: matrix, zero, entries 2241 2242 .seealso: MatZeroRows() 2243 @*/ 2244 int MatZeroEntries(Mat mat) 2245 { 2246 int ierr; 2247 2248 PetscFunctionBegin; 2249 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2250 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2251 if (!mat->ops.zeroentries) SETERRQ(PETSC_ERR_SUP,0,""); 2252 2253 PLogEventBegin(MAT_ZeroEntries,mat,0,0,0); 2254 ierr = (*mat->ops.zeroentries)(mat); CHKERRQ(ierr); 2255 PLogEventEnd(MAT_ZeroEntries,mat,0,0,0); 2256 PetscFunctionReturn(0); 2257 } 2258 2259 #undef __FUNC__ 2260 #define __FUNC__ "MatZeroRows" 2261 /*@ 2262 MatZeroRows - Zeros all entries (except possibly the main diagonal) 2263 of a set of rows of a matrix. 2264 2265 Input Parameters: 2266 . mat - the matrix 2267 . is - index set of rows to remove 2268 . diag - pointer to value put in all diagonals of eliminated rows. 2269 Note that diag is not a pointer to an array, but merely a 2270 pointer to a single value. 2271 2272 Notes: 2273 For the AIJ matrix formats this removes the old nonzero structure, 2274 but does not release memory. For the dense and block diagonal 2275 formats this does not alter the nonzero structure. 2276 2277 The user can set a value in the diagonal entry (or for the AIJ and 2278 row formats can optionally remove the main diagonal entry from the 2279 nonzero structure as well, by passing a null pointer as the final 2280 argument). 2281 2282 For the parallel case, all processes that share the matrix (i.e., 2283 those in the communicator used for matrix creation) MUST call this 2284 routine, regardless of whether any rows being zeroed are owned by 2285 them. 2286 2287 .keywords: matrix, zero, rows, boundary conditions 2288 2289 .seealso: MatZeroEntries(), 2290 @*/ 2291 int MatZeroRows(Mat mat,IS is, Scalar *diag) 2292 { 2293 int ierr; 2294 2295 PetscFunctionBegin; 2296 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2297 PetscValidHeaderSpecific(is,IS_COOKIE); 2298 if (diag) PetscValidScalarPointer(diag); 2299 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2300 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2301 if (!mat->ops.zerorows) SETERRQ(PETSC_ERR_SUP,0,""); 2302 2303 ierr = (*mat->ops.zerorows)(mat,is,diag); CHKERRQ(ierr); 2304 ierr = MatView_Private(mat); CHKERRQ(ierr); 2305 PetscFunctionReturn(0); 2306 } 2307 2308 #undef __FUNC__ 2309 #define __FUNC__ "MatZeroRowsLocal" 2310 /*@ 2311 MatZeroRowsLocal - Zeros all entries (except possibly the main diagonal) 2312 of a set of rows of a matrix; using local numbering of rows. 2313 2314 Input Parameters: 2315 . mat - the matrix 2316 . is - index set of rows to remove 2317 . diag - pointer to value put in all diagonals of eliminated rows. 2318 Note that diag is not a pointer to an array, but merely a 2319 pointer to a single value. 2320 2321 Notes: 2322 For the AIJ matrix formats this removes the old nonzero structure, 2323 but does not release memory. For the dense and block diagonal 2324 formats this does not alter the nonzero structure. 2325 2326 The user can set a value in the diagonal entry (or for the AIJ and 2327 row formats can optionally remove the main diagonal entry from the 2328 nonzero structure as well, by passing a null pointer as the final 2329 argument). 2330 2331 .keywords: matrix, zero, rows, boundary conditions 2332 2333 .seealso: MatZeroEntries(), 2334 @*/ 2335 int MatZeroRowsLocal(Mat mat,IS is, Scalar *diag) 2336 { 2337 int ierr; 2338 IS newis; 2339 2340 PetscFunctionBegin; 2341 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2342 PetscValidHeaderSpecific(is,IS_COOKIE); 2343 if (diag) PetscValidScalarPointer(diag); 2344 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2345 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2346 if (!mat->ops.zerorows) SETERRQ(PETSC_ERR_SUP,0,""); 2347 2348 ierr = ISLocalToGlobalMappingApplyIS(mat->mapping,is,&newis); CHKERRQ(ierr); 2349 ierr = (*mat->ops.zerorows)(mat,newis,diag); CHKERRQ(ierr); 2350 ierr = ISDestroy(newis); 2351 PetscFunctionReturn(0); 2352 } 2353 2354 #undef __FUNC__ 2355 #define __FUNC__ "MatGetSize" 2356 /*@ 2357 MatGetSize - Returns the numbers of rows and columns in a matrix. 2358 2359 Input Parameter: 2360 . mat - the matrix 2361 2362 Output Parameters: 2363 . m - the number of global rows 2364 . n - the number of global columns 2365 2366 .keywords: matrix, dimension, size, rows, columns, global, get 2367 2368 .seealso: MatGetLocalSize() 2369 @*/ 2370 int MatGetSize(Mat mat,int *m,int* n) 2371 { 2372 int ierr; 2373 2374 PetscFunctionBegin; 2375 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2376 ierr = (*mat->ops.getsize)(mat,m,n);CHKERRQ(ierr); 2377 PetscFunctionReturn(0); 2378 } 2379 2380 #undef __FUNC__ 2381 #define __FUNC__ "MatGetLocalSize" 2382 /*@ 2383 MatGetLocalSize - Returns the number of rows and columns in a matrix 2384 stored locally. This information may be implementation dependent, so 2385 use with care. 2386 2387 Input Parameters: 2388 . mat - the matrix 2389 2390 Output Parameters: 2391 . m - the number of local rows 2392 . n - the number of local columns 2393 2394 .keywords: matrix, dimension, size, local, rows, columns, get 2395 2396 .seealso: MatGetSize() 2397 @*/ 2398 int MatGetLocalSize(Mat mat,int *m,int* n) 2399 { 2400 int ierr; 2401 2402 PetscFunctionBegin; 2403 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2404 ierr = (*mat->ops.getlocalsize)(mat,m,n);CHKERRQ(ierr); 2405 PetscFunctionReturn(0); 2406 } 2407 2408 #undef __FUNC__ 2409 #define __FUNC__ "MatGetOwnershipRange" 2410 /*@ 2411 MatGetOwnershipRange - Returns the range of matrix rows owned by 2412 this processor, assuming that the matrix is laid out with the first 2413 n1 rows on the first processor, the next n2 rows on the second, etc. 2414 For certain parallel layouts this range may not be well defined. 2415 2416 Input Parameters: 2417 . mat - the matrix 2418 2419 Output Parameters: 2420 . m - the global index of the first local row 2421 . n - one more than the global index of the last local row 2422 2423 .keywords: matrix, get, range, ownership 2424 @*/ 2425 int MatGetOwnershipRange(Mat mat,int *m,int* n) 2426 { 2427 int ierr; 2428 2429 PetscFunctionBegin; 2430 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2431 PetscValidIntPointer(m); 2432 PetscValidIntPointer(n); 2433 if (!mat->ops.getownershiprange) SETERRQ(PETSC_ERR_SUP,0,""); 2434 ierr = (*mat->ops.getownershiprange)(mat,m,n);CHKERRQ(ierr); 2435 PetscFunctionReturn(0); 2436 } 2437 2438 #undef __FUNC__ 2439 #define __FUNC__ "MatILUFactorSymbolic" 2440 /*@ 2441 MatILUFactorSymbolic - Performs symbolic ILU factorization of a matrix. 2442 Uses levels of fill only, not drop tolerance. Use MatLUFactorNumeric() 2443 to complete the factorization. 2444 2445 Input Parameters: 2446 . mat - the matrix 2447 . row - row permutation 2448 . column - column permutation 2449 . fill - number of levels of fill 2450 . f - expected fill as ratio of the original number of nonzeros, 2451 for example 3.0; choosing this parameter well can result in 2452 more efficient use of time and space. Run your code with -log_info 2453 to determine an optimal value to use. 2454 2455 Output Parameters: 2456 . fact - new matrix that has been symbolically factored 2457 2458 Notes: 2459 See the file $(PETSC_DIR)/Performace for additional information about 2460 choosing the fill factor for better efficiency. 2461 2462 .keywords: matrix, factor, incomplete, ILU, symbolic, fill 2463 2464 .seealso: MatLUFactorSymbolic(), MatLUFactorNumeric() 2465 @*/ 2466 int MatILUFactorSymbolic(Mat mat,IS row,IS col,double f,int fill,Mat *fact) 2467 { 2468 int ierr; 2469 2470 PetscFunctionBegin; 2471 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2472 PetscValidPointer(fact); 2473 if (fill < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Levels of fill negative"); 2474 if (!mat->ops.ilufactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support parallel ILU"); 2475 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2476 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2477 2478 PLogEventBegin(MAT_ILUFactorSymbolic,mat,row,col,0); 2479 ierr = (*mat->ops.ilufactorsymbolic)(mat,row,col,f,fill,fact); CHKERRQ(ierr); 2480 PLogEventEnd(MAT_ILUFactorSymbolic,mat,row,col,0); 2481 PetscFunctionReturn(0); 2482 } 2483 2484 #undef __FUNC__ 2485 #define __FUNC__ "MatIncompleteCholeskyFactorSymbolic" 2486 /*@ 2487 MatIncompleteCholeskyFactorSymbolic - Performs symbolic incomplete 2488 Cholesky factorization for a symmetric matrix. Use 2489 MatCholeskyFactorNumeric() to complete the factorization. 2490 2491 Input Parameters: 2492 . mat - the matrix 2493 . perm - row and column permutation 2494 . fill - levels of fill 2495 . f - expected fill as ratio of original fill 2496 2497 Output Parameter: 2498 . fact - the factored matrix 2499 2500 Note: Currently only no-fill factorization is supported. 2501 2502 .keywords: matrix, factor, incomplete, ICC, Cholesky, symbolic, fill 2503 2504 .seealso: MatCholeskyFactorNumeric(), MatCholeskyFactor() 2505 @*/ 2506 int MatIncompleteCholeskyFactorSymbolic(Mat mat,IS perm,double f,int fill, 2507 Mat *fact) 2508 { 2509 int ierr; 2510 2511 PetscFunctionBegin; 2512 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2513 PetscValidPointer(fact); 2514 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2515 if (fill < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Fill negative"); 2516 if (!mat->ops.incompletecholeskyfactorsymbolic) SETERRQ(PETSC_ERR_SUP,0,"Currently only MatCreateMPIRowbs() matrices support ICC in parallel"); 2517 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2518 2519 PLogEventBegin(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0); 2520 ierr = (*mat->ops.incompletecholeskyfactorsymbolic)(mat,perm,f,fill,fact);CHKERRQ(ierr); 2521 PLogEventEnd(MAT_IncompleteCholeskyFactorSymbolic,mat,perm,0,0); 2522 PetscFunctionReturn(0); 2523 } 2524 2525 #undef __FUNC__ 2526 #define __FUNC__ "MatGetArray" 2527 /*@C 2528 MatGetArray - Returns a pointer to the element values in the matrix. 2529 This routine is implementation dependent, and may not even work for 2530 certain matrix types. You MUST call MatRestoreArray() when you no 2531 longer need to access the array. 2532 2533 Input Parameter: 2534 . mat - the matrix 2535 2536 Output Parameter: 2537 . v - the location of the values 2538 2539 Fortran Note: 2540 The Fortran interface is slightly different from that given below. 2541 See the Fortran chapter of the users manual and 2542 petsc/src/mat/examples for details. 2543 2544 .keywords: matrix, array, elements, values 2545 2546 .seealso: MatRestoreArray() 2547 @*/ 2548 int MatGetArray(Mat mat,Scalar **v) 2549 { 2550 int ierr; 2551 2552 PetscFunctionBegin; 2553 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2554 PetscValidPointer(v); 2555 if (!mat->ops.getarray) SETERRQ(PETSC_ERR_SUP,0,""); 2556 ierr = (*mat->ops.getarray)(mat,v);CHKERRQ(ierr); 2557 PetscFunctionReturn(0); 2558 } 2559 2560 #undef __FUNC__ 2561 #define __FUNC__ "MatRestoreArray" 2562 /*@C 2563 MatRestoreArray - Restores the matrix after MatGetArray has been called. 2564 2565 Input Parameter: 2566 . mat - the matrix 2567 . v - the location of the values 2568 2569 Fortran Note: 2570 The Fortran interface is slightly different from that given below. 2571 See the users manual and petsc/src/mat/examples for details. 2572 2573 .keywords: matrix, array, elements, values, restore 2574 2575 .seealso: MatGetArray() 2576 @*/ 2577 int MatRestoreArray(Mat mat,Scalar **v) 2578 { 2579 int ierr; 2580 2581 PetscFunctionBegin; 2582 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2583 PetscValidPointer(v); 2584 if (!mat->ops.restorearray) SETERRQ(PETSC_ERR_SUP,0,""); 2585 ierr = (*mat->ops.restorearray)(mat,v);CHKERRQ(ierr); 2586 PetscFunctionReturn(0); 2587 } 2588 2589 #undef __FUNC__ 2590 #define __FUNC__ "MatGetSubMatrices" 2591 /*@C 2592 MatGetSubMatrices - Extracts several submatrices from a matrix. If submat 2593 points to an array of valid matrices, it may be reused. 2594 2595 Input Parameters: 2596 . mat - the matrix 2597 . n - the number of submatrixes to be extracted 2598 . irow, icol - index sets of rows and columns to extract 2599 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 2600 2601 Output Parameter: 2602 . submat - the array of submatrices 2603 2604 Notes: 2605 MatGetSubMatrices() can extract only sequential submatrices 2606 (from both sequential and parallel matrices). Use MatGetSubMatrix() 2607 to extract a parallel submatrix. 2608 2609 When extracting submatrices from a parallel matrix, each processor can 2610 form a different submatrix by setting the rows and columns of its 2611 individual index sets according to the local submatrix desired. 2612 2613 When finished using the submatrices, the user should destroy 2614 them with MatDestroySubMatrices(). 2615 2616 .keywords: matrix, get, submatrix, submatrices 2617 2618 .seealso: MatDestroyMatrices(), MatGetSubMatrix() 2619 @*/ 2620 int MatGetSubMatrices(Mat mat,int n,IS *irow,IS *icol,MatGetSubMatrixCall scall, 2621 Mat **submat) 2622 { 2623 int ierr; 2624 2625 PetscFunctionBegin; 2626 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2627 if (!mat->ops.getsubmatrices) SETERRQ(PETSC_ERR_SUP,0,""); 2628 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2629 2630 PLogEventBegin(MAT_GetSubMatrices,mat,0,0,0); 2631 ierr = (*mat->ops.getsubmatrices)(mat,n,irow,icol,scall,submat); CHKERRQ(ierr); 2632 PLogEventEnd(MAT_GetSubMatrices,mat,0,0,0); 2633 2634 PetscFunctionReturn(0); 2635 } 2636 2637 #undef __FUNC__ 2638 #define __FUNC__ "MatDestroyMatrices" 2639 /*@C 2640 MatDestroyMatrices - Destroys a set of matrices obtained with MatGetSubMatrices(). 2641 2642 Input Parameters: 2643 . n - the number of local matrices 2644 . mat - the matrices 2645 2646 .keywords: matrix, destroy, submatrix, submatrices 2647 2648 .seealso: MatGetSubMatrices() 2649 @*/ 2650 int MatDestroyMatrices(int n,Mat **mat) 2651 { 2652 int ierr,i; 2653 2654 PetscFunctionBegin; 2655 if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,1,"Trying to destroy negative number of matrices"); 2656 PetscValidPointer(mat); 2657 for ( i=0; i<n; i++ ) { 2658 ierr = MatDestroy((*mat)[i]); CHKERRQ(ierr); 2659 } 2660 if (n) PetscFree(*mat); 2661 PetscFunctionReturn(0); 2662 } 2663 2664 #undef __FUNC__ 2665 #define __FUNC__ "MatIncreaseOverlap" 2666 /*@ 2667 MatIncreaseOverlap - Given a set of submatrices indicated by index sets, 2668 replaces the index sets by larger ones that represent submatrices with 2669 additional overlap. 2670 2671 Input Parameters: 2672 . mat - the matrix 2673 . n - the number of index sets 2674 . is - the array of pointers to index sets 2675 . ov - the additional overlap requested 2676 2677 .keywords: matrix, overlap, Schwarz 2678 2679 .seealso: MatGetSubMatrices() 2680 @*/ 2681 int MatIncreaseOverlap(Mat mat,int n, IS *is,int ov) 2682 { 2683 int ierr; 2684 2685 PetscFunctionBegin; 2686 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2687 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2688 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2689 2690 if (ov == 0) PetscFunctionReturn(0); 2691 if (!mat->ops.increaseoverlap) SETERRQ(PETSC_ERR_SUP,0,""); 2692 PLogEventBegin(MAT_IncreaseOverlap,mat,0,0,0); 2693 ierr = (*mat->ops.increaseoverlap)(mat,n,is,ov); CHKERRQ(ierr); 2694 PLogEventEnd(MAT_IncreaseOverlap,mat,0,0,0); 2695 PetscFunctionReturn(0); 2696 } 2697 2698 #undef __FUNC__ 2699 #define __FUNC__ "MatPrintHelp" 2700 /*@ 2701 MatPrintHelp - Prints all the options for the matrix. 2702 2703 Input Parameter: 2704 . mat - the matrix 2705 2706 Options Database Keys: 2707 $ -help, -h 2708 2709 .keywords: mat, help 2710 2711 .seealso: MatCreate(), MatCreateXXX() 2712 @*/ 2713 int MatPrintHelp(Mat mat) 2714 { 2715 static int called = 0; 2716 int ierr; 2717 MPI_Comm comm; 2718 2719 PetscFunctionBegin; 2720 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2721 2722 comm = mat->comm; 2723 if (!called) { 2724 PetscPrintf(comm,"General matrix options:\n"); 2725 PetscPrintf(comm," -mat_view_info: view basic matrix info during MatAssemblyEnd()\n"); 2726 PetscPrintf(comm," -mat_view_info_detailed: view detailed matrix info during MatAssemblyEnd()\n"); 2727 PetscPrintf(comm," -mat_view_draw: draw nonzero matrix structure during MatAssemblyEnd()\n"); 2728 PetscPrintf(comm," -draw_pause <sec>: set seconds of display pause\n"); 2729 PetscPrintf(comm," -display <name>: set alternate display\n"); 2730 called = 1; 2731 } 2732 if (mat->ops.printhelp) { 2733 ierr = (*mat->ops.printhelp)(mat); CHKERRQ(ierr); 2734 } 2735 PetscFunctionReturn(0); 2736 } 2737 2738 #undef __FUNC__ 2739 #define __FUNC__ "MatGetBlockSize" 2740 /*@ 2741 MatGetBlockSize - Returns the matrix block size; useful especially for the 2742 block row and block diagonal formats. 2743 2744 Input Parameter: 2745 . mat - the matrix 2746 2747 Output Parameter: 2748 . bs - block size 2749 2750 Notes: 2751 $ block diagonal formats: MATSEQBDIAG, MATMPIBDIAG 2752 $ block row formats: MATSEQBAIJ, MATMPIBAIJ 2753 2754 .keywords: matrix, get, block, size 2755 2756 .seealso: MatCreateSeqBAIJ(), MatCreateMPIBAIJ(), MatCreateSeqBDiag(), MatCreateMPIBDiag() 2757 @*/ 2758 int MatGetBlockSize(Mat mat,int *bs) 2759 { 2760 int ierr; 2761 2762 PetscFunctionBegin; 2763 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2764 PetscValidIntPointer(bs); 2765 if (!mat->ops.getblocksize) SETERRQ(PETSC_ERR_SUP,0,""); 2766 ierr = (*mat->ops.getblocksize)(mat,bs);CHKERRQ(ierr); 2767 PetscFunctionReturn(0); 2768 } 2769 2770 #undef __FUNC__ 2771 #define __FUNC__ "MatGetRowIJ" 2772 /*@C 2773 MatGetRowIJ - Returns the compress row storage i and j indices for sequential matrices. 2774 EXPERTS ONLY. 2775 2776 Input Parameters: 2777 . mat - the matrix 2778 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2779 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2780 symmetrized 2781 2782 Output Parameters: 2783 . n - number of rows and columns in the (possibly compressed) matrix 2784 . ia - the row indices 2785 . ja - the column indices 2786 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2787 @*/ 2788 int MatGetRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2789 { 2790 int ierr; 2791 2792 PetscFunctionBegin; 2793 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2794 if (ia) PetscValidIntPointer(ia); 2795 if (ja) PetscValidIntPointer(ja); 2796 PetscValidIntPointer(done); 2797 if (!mat->ops.getrowij) *done = PETSC_FALSE; 2798 else { 2799 *done = PETSC_TRUE; 2800 ierr = (*mat->ops.getrowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2801 } 2802 PetscFunctionReturn(0); 2803 } 2804 2805 #undef __FUNC__ 2806 #define __FUNC__ "MatGetColumnIJ" 2807 /*@C 2808 MatGetColumnIJ - Returns the compress Column storage i and j indices for sequential matrices. 2809 EXPERTS ONLY. 2810 2811 Input Parameters: 2812 . mat - the matrix 2813 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2814 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2815 symmetrized 2816 2817 Output Parameters: 2818 . n - number of Columns and columns in the (possibly compressed) matrix 2819 . ia - the Column indices 2820 . ja - the column indices 2821 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2822 @*/ 2823 int MatGetColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2824 { 2825 int ierr; 2826 2827 PetscFunctionBegin; 2828 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2829 if (ia) PetscValidIntPointer(ia); 2830 if (ja) PetscValidIntPointer(ja); 2831 PetscValidIntPointer(done); 2832 2833 if (!mat->ops.getcolumnij) *done = PETSC_FALSE; 2834 else { 2835 *done = PETSC_TRUE; 2836 ierr = (*mat->ops.getcolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2837 } 2838 PetscFunctionReturn(0); 2839 } 2840 2841 #undef __FUNC__ 2842 #define __FUNC__ "MatRestoreRowIJ" 2843 /*@C 2844 MatRestoreRowIJ - Call after you are completed with the ia,ja indices obtained with 2845 MatGetRowIJ(). EXPERTS ONLY. 2846 2847 Input Parameters: 2848 . mat - the matrix 2849 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2850 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2851 symmetrized 2852 2853 Output Parameters: 2854 . n - size of (possibly compressed) matrix 2855 . ia - the row indices 2856 . ja - the column indices 2857 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2858 @*/ 2859 int MatRestoreRowIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2860 { 2861 int ierr; 2862 2863 PetscFunctionBegin; 2864 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2865 if (ia) PetscValidIntPointer(ia); 2866 if (ja) PetscValidIntPointer(ja); 2867 PetscValidIntPointer(done); 2868 2869 if (!mat->ops.restorerowij) *done = PETSC_FALSE; 2870 else { 2871 *done = PETSC_TRUE; 2872 ierr = (*mat->ops.restorerowij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2873 } 2874 PetscFunctionReturn(0); 2875 } 2876 2877 #undef __FUNC__ 2878 #define __FUNC__ "MatRestoreColumnIJ" 2879 /*@C 2880 MatRestoreColumnIJ - Call after you are completed with the ia,ja indices obtained with 2881 MatGetColumnIJ(). EXPERTS ONLY. 2882 2883 Input Parameters: 2884 . mat - the matrix 2885 . shift - 1 or zero indicating we want the indices starting at 0 or 1 2886 . symmetric - PETSC_TRUE or PETSC_FALSE indicating the matrix data structure should be 2887 symmetrized 2888 2889 Output Parameters: 2890 . n - size of (possibly compressed) matrix 2891 . ia - the Column indices 2892 . ja - the column indices 2893 . done - PETSC_TRUE or PETSC_FALSE indicated that the values have been returned 2894 @*/ 2895 int MatRestoreColumnIJ(Mat mat,int shift,PetscTruth symmetric,int *n,int **ia,int** ja,PetscTruth *done) 2896 { 2897 int ierr; 2898 2899 PetscFunctionBegin; 2900 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2901 if (ia) PetscValidIntPointer(ia); 2902 if (ja) PetscValidIntPointer(ja); 2903 PetscValidIntPointer(done); 2904 2905 if (!mat->ops.restorecolumnij) *done = PETSC_FALSE; 2906 else { 2907 *done = PETSC_TRUE; 2908 ierr = (*mat->ops.restorecolumnij)(mat,shift,symmetric,n,ia,ja,done); CHKERRQ(ierr); 2909 } 2910 PetscFunctionReturn(0); 2911 } 2912 2913 #undef __FUNC__ 2914 #define __FUNC__ "MatColoringPatch" 2915 /*@C 2916 MatColoringPatch - EXPERTS ONLY, used inside matrix coloring routines that 2917 use matGetRowIJ() and/or MatGetColumnIJ(). 2918 2919 Input Parameters: 2920 . mat - the matrix 2921 . n - number of colors 2922 . colorarray - array indicating color for each column 2923 2924 Output Parameters: 2925 . iscoloring - coloring generated using colorarray information 2926 2927 @*/ 2928 int MatColoringPatch(Mat mat,int n,int *colorarray,ISColoring *iscoloring) 2929 { 2930 int ierr; 2931 2932 PetscFunctionBegin; 2933 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2934 PetscValidIntPointer(colorarray); 2935 2936 if (!mat->ops.coloringpatch) {SETERRQ(PETSC_ERR_SUP,0,"");} 2937 else { 2938 ierr = (*mat->ops.coloringpatch)(mat,n,colorarray,iscoloring); CHKERRQ(ierr); 2939 } 2940 PetscFunctionReturn(0); 2941 } 2942 2943 2944 #undef __FUNC__ 2945 #define __FUNC__ "MatSetUnfactored" 2946 /*@ 2947 MatSetUnfactored - Resets a factored matrix to be treated as unfactored. 2948 2949 Input Paramter: 2950 . mat - the factored matrix to be reset 2951 2952 Notes: 2953 This routine should be used only with factored matrices formed by in-place 2954 factorization via ILU(0) (or by in-place LU factorization for the MATSEQDENSE 2955 format). This option can save memory, for example, when solving nonlinear 2956 systems with a matrix-free Newton-Krylov method and a matrix-based, in-place 2957 ILU(0) preconditioner. 2958 2959 Note that one can specify in-place ILU(0) factorization by calling 2960 $ PCType(pc,PCILU); 2961 $ PCILUSeUseInPlace(pc); 2962 or by using the options -pc_type ilu -pc_ilu_in_place 2963 2964 In-place factorization ILU(0) can also be used as a local 2965 solver for the blocks within the block Jacobi or additive Schwarz 2966 methods (runtime option: -sub_pc_ilu_in_place). See the discussion 2967 of these preconditioners in the users manual for details on setting 2968 local solver options. 2969 2970 .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace() 2971 2972 .keywords: matrix-free, in-place ILU, in-place LU 2973 @*/ 2974 int MatSetUnfactored(Mat mat) 2975 { 2976 int ierr; 2977 2978 PetscFunctionBegin; 2979 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2980 mat->factor = 0; 2981 if (!mat->ops.setunfactored) PetscFunctionReturn(0); 2982 ierr = (*mat->ops.setunfactored)(mat); CHKERRQ(ierr); 2983 PetscFunctionReturn(0); 2984 } 2985 2986 #undef __FUNC__ 2987 #define __FUNC__ "MatGetType" 2988 /*@C 2989 MatGetType - Gets the matrix type and name (as a string) from the matrix. 2990 2991 Input Parameter: 2992 . mat - the matrix 2993 2994 Output Parameter: 2995 . type - the matrix type (or use PETSC_NULL) 2996 . name - name of matrix type (or use PETSC_NULL) 2997 2998 .keywords: matrix, get, type, name 2999 @*/ 3000 int MatGetType(Mat mat,MatType *type,char **name) 3001 { 3002 int itype = (int)mat->type; 3003 char *matname[10]; 3004 3005 PetscFunctionBegin; 3006 PetscValidHeaderSpecific(mat,MAT_COOKIE); 3007 3008 if (type) *type = (MatType) mat->type; 3009 if (name) { 3010 /* Note: Be sure that this list corresponds to the enum in mat.h */ 3011 matname[0] = "MATSEQDENSE"; 3012 matname[1] = "MATSEQAIJ"; 3013 matname[2] = "MATMPIAIJ"; 3014 matname[3] = "MATSHELL"; 3015 matname[4] = "MATMPIROWBS"; 3016 matname[5] = "MATSEQBDIAG"; 3017 matname[6] = "MATMPIBDIAG"; 3018 matname[7] = "MATMPIDENSE"; 3019 matname[8] = "MATSEQBAIJ"; 3020 matname[9] = "MATMPIBAIJ"; 3021 3022 if (itype < 0 || itype > 9) *name = "Unknown matrix type"; 3023 else *name = matname[itype]; 3024 } 3025 PetscFunctionReturn(0); 3026 } 3027 3028 /*MC 3029 MatGetArrayF90 - Accesses a matrix array from Fortran90. 3030 3031 Input Parameter: 3032 . x - matrix 3033 3034 Output Parameter: 3035 . xx_v - the Fortran90 pointer to the array 3036 . ierr - error code 3037 3038 Synopsis: 3039 MatGetArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 3040 3041 Usage: 3042 $ Scalar, pointer xx_v(:) 3043 $ .... 3044 $ call MatGetArrayF90(x,xx_v,ierr) 3045 $ a = xx_v(3) 3046 $ call MatRestoreArrayF90(x,xx_v,ierr) 3047 3048 Notes: 3049 Currently only supported using the NAG F90 compiler. 3050 3051 .seealso: MatRestoreArrayF90(), MatGetArray(), MatRestoreArray() 3052 3053 .keywords: matrix, array, f90 3054 M*/ 3055 3056 /*MC 3057 MatRestoreArrayF90 - Restores a matrix array that has been 3058 accessed with MatGetArrayF90(). 3059 3060 Input Parameters: 3061 . x - matrix 3062 . xx_v - the Fortran90 pointer to the array 3063 3064 Output Parameter: 3065 . ierr - error code 3066 3067 Synopsis: 3068 MatRestoreArrayF90(Mat x,{Scalar, pointer :: xx_v(:)},integer ierr) 3069 3070 Example of Usage: 3071 $ Scalar, pointer xx_v(:) 3072 $ .... 3073 $ call MatGetArrayF90(x,xx_v,ierr) 3074 $ a = xx_v(3) 3075 $ call MatRestoreArrayF90(x,xx_v,ierr) 3076 3077 Notes: 3078 Currently only supported using the NAG F90 compiler. 3079 3080 .seealso: MatGetArrayF90(), MatGetArray(), MatRestoreArray() 3081 3082 .keywords: matrix, array, f90 3083 M*/ 3084 3085 3086 #undef __FUNC__ 3087 #define __FUNC__ "MatGetSubMatrix" 3088 /*@ 3089 MatGetSubMatrix - Gets a single submatrix on the same number of processors 3090 as the original matrix. 3091 3092 Input Parameters: 3093 . mat - the original matrix 3094 . isrow - rows this processor should obtain 3095 . iscol - columns for all processors you wish kept 3096 . cll - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 3097 3098 Output Parameters: 3099 . newmat - the new submatrix, of the same type as the old 3100 3101 .seealso: MatGetSubMatrices() 3102 3103 @*/ 3104 int MatGetSubMatrix(Mat mat,IS isrow,IS iscol,MatGetSubMatrixCall cll,Mat *newmat) 3105 { 3106 int ierr, size; 3107 Mat *local; 3108 3109 PetscFunctionBegin; 3110 MPI_Comm_size(mat->comm,&size); 3111 3112 /* if original matrix is on just one processor then use submatrix generated */ 3113 if (size == 1 && cll == MAT_REUSE_MATRIX) { 3114 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_REUSE_MATRIX,&newmat);CHKERRQ(ierr); 3115 PetscFunctionReturn(0); 3116 } else if (size == 1) { 3117 ierr = MatGetSubMatrices(mat,1,&isrow,&iscol,MAT_INITIAL_MATRIX,&local);CHKERRQ(ierr); 3118 *newmat = *local; 3119 PetscFree(local); 3120 PetscFunctionReturn(0); 3121 } 3122 3123 if (!mat->ops.getsubmatrix) SETERRQ(PETSC_ERR_SUP,0,"Not currently implemented"); 3124 ierr = (*mat->ops.getsubmatrix)(mat,isrow,iscol,cll,newmat);CHKERRQ(ierr); 3125 PetscFunctionReturn(0); 3126 } 3127 3128 3129 3130