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