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