1 2 /* 3 Provides an interface to the MUMPS sparse solver 4 */ 5 6 #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I "petscmat.h" I*/ 7 #include <../src/mat/impls/sbaij/mpi/mpisbaij.h> 8 9 EXTERN_C_BEGIN 10 #if defined(PETSC_USE_COMPLEX) 11 #if defined(PETSC_USE_REAL_SINGLE) 12 #include <cmumps_c.h> 13 #else 14 #include <zmumps_c.h> 15 #endif 16 #else 17 #if defined(PETSC_USE_REAL_SINGLE) 18 #include <smumps_c.h> 19 #else 20 #include <dmumps_c.h> 21 #endif 22 #endif 23 EXTERN_C_END 24 #define JOB_INIT -1 25 #define JOB_FACTSYMBOLIC 1 26 #define JOB_FACTNUMERIC 2 27 #define JOB_SOLVE 3 28 #define JOB_END -2 29 30 /* calls to MUMPS */ 31 #if defined(PETSC_USE_COMPLEX) 32 #if defined(PETSC_USE_REAL_SINGLE) 33 #define PetscMUMPS_c cmumps_c 34 #else 35 #define PetscMUMPS_c zmumps_c 36 #endif 37 #else 38 #if defined(PETSC_USE_REAL_SINGLE) 39 #define PetscMUMPS_c smumps_c 40 #else 41 #define PetscMUMPS_c dmumps_c 42 #endif 43 #endif 44 45 46 /* macros s.t. indices match MUMPS documentation */ 47 #define ICNTL(I) icntl[(I)-1] 48 #define CNTL(I) cntl[(I)-1] 49 #define INFOG(I) infog[(I)-1] 50 #define INFO(I) info[(I)-1] 51 #define RINFOG(I) rinfog[(I)-1] 52 #define RINFO(I) rinfo[(I)-1] 53 54 typedef struct { 55 #if defined(PETSC_USE_COMPLEX) 56 #if defined(PETSC_USE_REAL_SINGLE) 57 CMUMPS_STRUC_C id; 58 #else 59 ZMUMPS_STRUC_C id; 60 #endif 61 #else 62 #if defined(PETSC_USE_REAL_SINGLE) 63 SMUMPS_STRUC_C id; 64 #else 65 DMUMPS_STRUC_C id; 66 #endif 67 #endif 68 69 MatStructure matstruc; 70 PetscMPIInt myid,size; 71 PetscInt *irn,*jcn,nz,sym; 72 PetscScalar *val; 73 MPI_Comm comm_mumps; 74 VecScatter scat_rhs, scat_sol; 75 PetscBool isAIJ,CleanUpMUMPS; 76 Vec b_seq,x_seq; 77 PetscInt ICNTL9_pre; /* check if ICNTL(9) is changed from previous MatSolve */ 78 79 PetscErrorCode (*Destroy)(Mat); 80 PetscErrorCode (*ConvertToTriples)(Mat, int, MatReuse, int*, int**, int**, PetscScalar**); 81 } Mat_MUMPS; 82 83 extern PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*); 84 85 86 /* MatConvertToTriples_A_B */ 87 /*convert Petsc matrix to triples: row[nz], col[nz], val[nz] */ 88 /* 89 input: 90 A - matrix in aij,baij or sbaij (bs=1) format 91 shift - 0: C style output triple; 1: Fortran style output triple. 92 reuse - MAT_INITIAL_MATRIX: spaces are allocated and values are set for the triple 93 MAT_REUSE_MATRIX: only the values in v array are updated 94 output: 95 nnz - dim of r, c, and v (number of local nonzero entries of A) 96 r, c, v - row and col index, matrix values (matrix triples) 97 98 The returned values r, c, and sometimes v are obtained in a single PetscMalloc(). Then in MatDestroy_MUMPS() it is 99 freed with PetscFree((mumps->irn); This is not ideal code, the fact that v is ONLY sometimes part of mumps->irn means 100 that the PetscMalloc() cannot easily be replaced with a PetscMalloc3(). 101 102 */ 103 104 #undef __FUNCT__ 105 #define __FUNCT__ "MatConvertToTriples_seqaij_seqaij" 106 PetscErrorCode MatConvertToTriples_seqaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 107 { 108 const PetscInt *ai,*aj,*ajj,M=A->rmap->n; 109 PetscInt nz,rnz,i,j; 110 PetscErrorCode ierr; 111 PetscInt *row,*col; 112 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 113 114 PetscFunctionBegin; 115 *v=aa->a; 116 if (reuse == MAT_INITIAL_MATRIX) { 117 nz = aa->nz; 118 ai = aa->i; 119 aj = aa->j; 120 *nnz = nz; 121 ierr = PetscMalloc1(2*nz, &row);CHKERRQ(ierr); 122 col = row + nz; 123 124 nz = 0; 125 for (i=0; i<M; i++) { 126 rnz = ai[i+1] - ai[i]; 127 ajj = aj + ai[i]; 128 for (j=0; j<rnz; j++) { 129 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 130 } 131 } 132 *r = row; *c = col; 133 } 134 PetscFunctionReturn(0); 135 } 136 137 #undef __FUNCT__ 138 #define __FUNCT__ "MatConvertToTriples_seqbaij_seqaij" 139 PetscErrorCode MatConvertToTriples_seqbaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 140 { 141 Mat_SeqBAIJ *aa=(Mat_SeqBAIJ*)A->data; 142 const PetscInt *ai,*aj,*ajj,bs2 = aa->bs2; 143 PetscInt bs,M,nz,idx=0,rnz,i,j,k,m; 144 PetscErrorCode ierr; 145 PetscInt *row,*col; 146 147 PetscFunctionBegin; 148 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 149 M = A->rmap->N/bs; 150 *v = aa->a; 151 if (reuse == MAT_INITIAL_MATRIX) { 152 ai = aa->i; aj = aa->j; 153 nz = bs2*aa->nz; 154 *nnz = nz; 155 ierr = PetscMalloc1(2*nz, &row);CHKERRQ(ierr); 156 col = row + nz; 157 158 for (i=0; i<M; i++) { 159 ajj = aj + ai[i]; 160 rnz = ai[i+1] - ai[i]; 161 for (k=0; k<rnz; k++) { 162 for (j=0; j<bs; j++) { 163 for (m=0; m<bs; m++) { 164 row[idx] = i*bs + m + shift; 165 col[idx++] = bs*(ajj[k]) + j + shift; 166 } 167 } 168 } 169 } 170 *r = row; *c = col; 171 } 172 PetscFunctionReturn(0); 173 } 174 175 #undef __FUNCT__ 176 #define __FUNCT__ "MatConvertToTriples_seqsbaij_seqsbaij" 177 PetscErrorCode MatConvertToTriples_seqsbaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 178 { 179 const PetscInt *ai, *aj,*ajj,M=A->rmap->n; 180 PetscInt nz,rnz,i,j; 181 PetscErrorCode ierr; 182 PetscInt *row,*col; 183 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data; 184 185 PetscFunctionBegin; 186 *v = aa->a; 187 if (reuse == MAT_INITIAL_MATRIX) { 188 nz = aa->nz; 189 ai = aa->i; 190 aj = aa->j; 191 *v = aa->a; 192 *nnz = nz; 193 ierr = PetscMalloc1(2*nz, &row);CHKERRQ(ierr); 194 col = row + nz; 195 196 nz = 0; 197 for (i=0; i<M; i++) { 198 rnz = ai[i+1] - ai[i]; 199 ajj = aj + ai[i]; 200 for (j=0; j<rnz; j++) { 201 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 202 } 203 } 204 *r = row; *c = col; 205 } 206 PetscFunctionReturn(0); 207 } 208 209 #undef __FUNCT__ 210 #define __FUNCT__ "MatConvertToTriples_seqaij_seqsbaij" 211 PetscErrorCode MatConvertToTriples_seqaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 212 { 213 const PetscInt *ai,*aj,*ajj,*adiag,M=A->rmap->n; 214 PetscInt nz,rnz,i,j; 215 const PetscScalar *av,*v1; 216 PetscScalar *val; 217 PetscErrorCode ierr; 218 PetscInt *row,*col; 219 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 220 221 PetscFunctionBegin; 222 ai =aa->i; aj=aa->j;av=aa->a; 223 adiag=aa->diag; 224 if (reuse == MAT_INITIAL_MATRIX) { 225 /* count nz in the uppper triangular part of A */ 226 nz = 0; 227 for (i=0; i<M; i++) nz += ai[i+1] - adiag[i]; 228 *nnz = nz; 229 230 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 231 col = row + nz; 232 val = (PetscScalar*)(col + nz); 233 234 nz = 0; 235 for (i=0; i<M; i++) { 236 rnz = ai[i+1] - adiag[i]; 237 ajj = aj + adiag[i]; 238 v1 = av + adiag[i]; 239 for (j=0; j<rnz; j++) { 240 row[nz] = i+shift; col[nz] = ajj[j] + shift; val[nz++] = v1[j]; 241 } 242 } 243 *r = row; *c = col; *v = val; 244 } else { 245 nz = 0; val = *v; 246 for (i=0; i <M; i++) { 247 rnz = ai[i+1] - adiag[i]; 248 ajj = aj + adiag[i]; 249 v1 = av + adiag[i]; 250 for (j=0; j<rnz; j++) { 251 val[nz++] = v1[j]; 252 } 253 } 254 } 255 PetscFunctionReturn(0); 256 } 257 258 #undef __FUNCT__ 259 #define __FUNCT__ "MatConvertToTriples_mpisbaij_mpisbaij" 260 PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 261 { 262 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 263 PetscErrorCode ierr; 264 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 265 PetscInt *row,*col; 266 const PetscScalar *av, *bv,*v1,*v2; 267 PetscScalar *val; 268 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 269 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)(mat->A)->data; 270 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 271 272 PetscFunctionBegin; 273 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 274 av=aa->a; bv=bb->a; 275 276 garray = mat->garray; 277 278 if (reuse == MAT_INITIAL_MATRIX) { 279 nz = aa->nz + bb->nz; 280 *nnz = nz; 281 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 282 col = row + nz; 283 val = (PetscScalar*)(col + nz); 284 285 *r = row; *c = col; *v = val; 286 } else { 287 row = *r; col = *c; val = *v; 288 } 289 290 jj = 0; irow = rstart; 291 for (i=0; i<m; i++) { 292 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 293 countA = ai[i+1] - ai[i]; 294 countB = bi[i+1] - bi[i]; 295 bjj = bj + bi[i]; 296 v1 = av + ai[i]; 297 v2 = bv + bi[i]; 298 299 /* A-part */ 300 for (j=0; j<countA; j++) { 301 if (reuse == MAT_INITIAL_MATRIX) { 302 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 303 } 304 val[jj++] = v1[j]; 305 } 306 307 /* B-part */ 308 for (j=0; j < countB; j++) { 309 if (reuse == MAT_INITIAL_MATRIX) { 310 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 311 } 312 val[jj++] = v2[j]; 313 } 314 irow++; 315 } 316 PetscFunctionReturn(0); 317 } 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpiaij" 321 PetscErrorCode MatConvertToTriples_mpiaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 322 { 323 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 324 PetscErrorCode ierr; 325 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 326 PetscInt *row,*col; 327 const PetscScalar *av, *bv,*v1,*v2; 328 PetscScalar *val; 329 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 330 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(mat->A)->data; 331 Mat_SeqAIJ *bb = (Mat_SeqAIJ*)(mat->B)->data; 332 333 PetscFunctionBegin; 334 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 335 av=aa->a; bv=bb->a; 336 337 garray = mat->garray; 338 339 if (reuse == MAT_INITIAL_MATRIX) { 340 nz = aa->nz + bb->nz; 341 *nnz = nz; 342 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 343 col = row + nz; 344 val = (PetscScalar*)(col + nz); 345 346 *r = row; *c = col; *v = val; 347 } else { 348 row = *r; col = *c; val = *v; 349 } 350 351 jj = 0; irow = rstart; 352 for (i=0; i<m; i++) { 353 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 354 countA = ai[i+1] - ai[i]; 355 countB = bi[i+1] - bi[i]; 356 bjj = bj + bi[i]; 357 v1 = av + ai[i]; 358 v2 = bv + bi[i]; 359 360 /* A-part */ 361 for (j=0; j<countA; j++) { 362 if (reuse == MAT_INITIAL_MATRIX) { 363 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 364 } 365 val[jj++] = v1[j]; 366 } 367 368 /* B-part */ 369 for (j=0; j < countB; j++) { 370 if (reuse == MAT_INITIAL_MATRIX) { 371 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 372 } 373 val[jj++] = v2[j]; 374 } 375 irow++; 376 } 377 PetscFunctionReturn(0); 378 } 379 380 #undef __FUNCT__ 381 #define __FUNCT__ "MatConvertToTriples_mpibaij_mpiaij" 382 PetscErrorCode MatConvertToTriples_mpibaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 383 { 384 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data; 385 Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)(mat->A)->data; 386 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 387 const PetscInt *ai = aa->i, *bi = bb->i, *aj = aa->j, *bj = bb->j,*ajj, *bjj; 388 const PetscInt *garray = mat->garray,mbs=mat->mbs,rstart=A->rmap->rstart; 389 const PetscInt bs2=mat->bs2; 390 PetscErrorCode ierr; 391 PetscInt bs,nz,i,j,k,n,jj,irow,countA,countB,idx; 392 PetscInt *row,*col; 393 const PetscScalar *av=aa->a, *bv=bb->a,*v1,*v2; 394 PetscScalar *val; 395 396 PetscFunctionBegin; 397 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 398 if (reuse == MAT_INITIAL_MATRIX) { 399 nz = bs2*(aa->nz + bb->nz); 400 *nnz = nz; 401 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 402 col = row + nz; 403 val = (PetscScalar*)(col + nz); 404 405 *r = row; *c = col; *v = val; 406 } else { 407 row = *r; col = *c; val = *v; 408 } 409 410 jj = 0; irow = rstart; 411 for (i=0; i<mbs; i++) { 412 countA = ai[i+1] - ai[i]; 413 countB = bi[i+1] - bi[i]; 414 ajj = aj + ai[i]; 415 bjj = bj + bi[i]; 416 v1 = av + bs2*ai[i]; 417 v2 = bv + bs2*bi[i]; 418 419 idx = 0; 420 /* A-part */ 421 for (k=0; k<countA; k++) { 422 for (j=0; j<bs; j++) { 423 for (n=0; n<bs; n++) { 424 if (reuse == MAT_INITIAL_MATRIX) { 425 row[jj] = irow + n + shift; 426 col[jj] = rstart + bs*ajj[k] + j + shift; 427 } 428 val[jj++] = v1[idx++]; 429 } 430 } 431 } 432 433 idx = 0; 434 /* B-part */ 435 for (k=0; k<countB; k++) { 436 for (j=0; j<bs; j++) { 437 for (n=0; n<bs; n++) { 438 if (reuse == MAT_INITIAL_MATRIX) { 439 row[jj] = irow + n + shift; 440 col[jj] = bs*garray[bjj[k]] + j + shift; 441 } 442 val[jj++] = v2[idx++]; 443 } 444 } 445 } 446 irow += bs; 447 } 448 PetscFunctionReturn(0); 449 } 450 451 #undef __FUNCT__ 452 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpisbaij" 453 PetscErrorCode MatConvertToTriples_mpiaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 454 { 455 const PetscInt *ai, *aj,*adiag, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 456 PetscErrorCode ierr; 457 PetscInt rstart,nz,nza,nzb,i,j,jj,irow,countA,countB; 458 PetscInt *row,*col; 459 const PetscScalar *av, *bv,*v1,*v2; 460 PetscScalar *val; 461 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 462 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(mat->A)->data; 463 Mat_SeqAIJ *bb =(Mat_SeqAIJ*)(mat->B)->data; 464 465 PetscFunctionBegin; 466 ai=aa->i; aj=aa->j; adiag=aa->diag; 467 bi=bb->i; bj=bb->j; garray = mat->garray; 468 av=aa->a; bv=bb->a; 469 470 rstart = A->rmap->rstart; 471 472 if (reuse == MAT_INITIAL_MATRIX) { 473 nza = 0; /* num of upper triangular entries in mat->A, including diagonals */ 474 nzb = 0; /* num of upper triangular entries in mat->B */ 475 for (i=0; i<m; i++) { 476 nza += (ai[i+1] - adiag[i]); 477 countB = bi[i+1] - bi[i]; 478 bjj = bj + bi[i]; 479 for (j=0; j<countB; j++) { 480 if (garray[bjj[j]] > rstart) nzb++; 481 } 482 } 483 484 nz = nza + nzb; /* total nz of upper triangular part of mat */ 485 *nnz = nz; 486 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 487 col = row + nz; 488 val = (PetscScalar*)(col + nz); 489 490 *r = row; *c = col; *v = val; 491 } else { 492 row = *r; col = *c; val = *v; 493 } 494 495 jj = 0; irow = rstart; 496 for (i=0; i<m; i++) { 497 ajj = aj + adiag[i]; /* ptr to the beginning of the diagonal of this row */ 498 v1 = av + adiag[i]; 499 countA = ai[i+1] - adiag[i]; 500 countB = bi[i+1] - bi[i]; 501 bjj = bj + bi[i]; 502 v2 = bv + bi[i]; 503 504 /* A-part */ 505 for (j=0; j<countA; j++) { 506 if (reuse == MAT_INITIAL_MATRIX) { 507 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 508 } 509 val[jj++] = v1[j]; 510 } 511 512 /* B-part */ 513 for (j=0; j < countB; j++) { 514 if (garray[bjj[j]] > rstart) { 515 if (reuse == MAT_INITIAL_MATRIX) { 516 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 517 } 518 val[jj++] = v2[j]; 519 } 520 } 521 irow++; 522 } 523 PetscFunctionReturn(0); 524 } 525 526 #undef __FUNCT__ 527 #define __FUNCT__ "MatGetDiagonal_MUMPS" 528 PetscErrorCode MatGetDiagonal_MUMPS(Mat A,Vec v) 529 { 530 PetscFunctionBegin; 531 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: MUMPS factor"); 532 PetscFunctionReturn(0); 533 } 534 535 #undef __FUNCT__ 536 #define __FUNCT__ "MatDestroy_MUMPS" 537 PetscErrorCode MatDestroy_MUMPS(Mat A) 538 { 539 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 540 PetscErrorCode ierr; 541 542 PetscFunctionBegin; 543 if (mumps->CleanUpMUMPS) { 544 /* Terminate instance, deallocate memories */ 545 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 546 ierr = VecScatterDestroy(&mumps->scat_rhs);CHKERRQ(ierr); 547 ierr = VecDestroy(&mumps->b_seq);CHKERRQ(ierr); 548 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 549 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 550 ierr = PetscFree(mumps->id.perm_in);CHKERRQ(ierr); 551 ierr = PetscFree(mumps->irn);CHKERRQ(ierr); 552 ierr = PetscFree2(mumps->id.listvar_schur,mumps->id.schur);CHKERRQ(ierr); 553 554 mumps->id.job = JOB_END; 555 PetscMUMPS_c(&mumps->id); 556 ierr = MPI_Comm_free(&(mumps->comm_mumps));CHKERRQ(ierr); 557 } 558 if (mumps->Destroy) { 559 ierr = (mumps->Destroy)(A);CHKERRQ(ierr); 560 } 561 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 562 563 /* clear composed functions */ 564 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);CHKERRQ(ierr); 565 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetIcntl_C",NULL);CHKERRQ(ierr); 566 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetIcntl_C",NULL);CHKERRQ(ierr); 567 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetCntl_C",NULL);CHKERRQ(ierr); 568 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetCntl_C",NULL);CHKERRQ(ierr); 569 570 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfo_C",NULL);CHKERRQ(ierr); 571 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfog_C",NULL);CHKERRQ(ierr); 572 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfo_C",NULL);CHKERRQ(ierr); 573 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfog_C",NULL);CHKERRQ(ierr); 574 575 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetSchurIndices_C",NULL);CHKERRQ(ierr); 576 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetSchurComplement_C",NULL);CHKERRQ(ierr); 577 PetscFunctionReturn(0); 578 } 579 580 #undef __FUNCT__ 581 #define __FUNCT__ "MatSolve_MUMPS" 582 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 583 { 584 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 585 PetscScalar *array; 586 Vec b_seq; 587 IS is_iden,is_petsc; 588 PetscErrorCode ierr; 589 PetscInt i; 590 static PetscBool cite1 = PETSC_FALSE,cite2 = PETSC_FALSE; 591 592 PetscFunctionBegin; 593 if (mumps->id.size_schur) { 594 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use MatSolve when Schur complement has been requested\n"); 595 } 596 ierr = PetscCitationsRegister("@article{MUMPS01,\n author = {P.~R. Amestoy and I.~S. Duff and J.-Y. L'Excellent and J. Koster},\n title = {A fully asynchronous multifrontal solver using distributed dynamic scheduling},\n journal = {SIAM Journal on Matrix Analysis and Applications},\n volume = {23},\n number = {1},\n pages = {15--41},\n year = {2001}\n}\n",&cite1);CHKERRQ(ierr); 597 ierr = PetscCitationsRegister("@article{MUMPS02,\n author = {P.~R. Amestoy and A. Guermouche and J.-Y. L'Excellent and S. Pralet},\n title = {Hybrid scheduling for the parallel solution of linear systems},\n journal = {Parallel Computing},\n volume = {32},\n number = {2},\n pages = {136--156},\n year = {2006}\n}\n",&cite2);CHKERRQ(ierr); 598 mumps->id.nrhs = 1; 599 b_seq = mumps->b_seq; 600 if (mumps->size > 1) { 601 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 602 ierr = VecScatterBegin(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 603 ierr = VecScatterEnd(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 604 if (!mumps->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 605 } else { /* size == 1 */ 606 ierr = VecCopy(b,x);CHKERRQ(ierr); 607 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 608 } 609 if (!mumps->myid) { /* define rhs on the host */ 610 mumps->id.nrhs = 1; 611 #if defined(PETSC_USE_COMPLEX) 612 #if defined(PETSC_USE_REAL_SINGLE) 613 mumps->id.rhs = (mumps_complex*)array; 614 #else 615 mumps->id.rhs = (mumps_double_complex*)array; 616 #endif 617 #else 618 mumps->id.rhs = array; 619 #endif 620 } 621 622 /* solve phase */ 623 /*-------------*/ 624 mumps->id.job = JOB_SOLVE; 625 PetscMUMPS_c(&mumps->id); 626 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 627 628 if (mumps->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 629 if (mumps->scat_sol && mumps->ICNTL9_pre != mumps->id.ICNTL(9)) { 630 /* when id.ICNTL(9) changes, the contents of lsol_loc may change (not its size, lsol_loc), recreates scat_sol */ 631 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 632 } 633 if (!mumps->scat_sol) { /* create scatter scat_sol */ 634 ierr = ISCreateStride(PETSC_COMM_SELF,mumps->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 635 for (i=0; i<mumps->id.lsol_loc; i++) { 636 mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 637 } 638 ierr = ISCreateGeneral(PETSC_COMM_SELF,mumps->id.lsol_loc,mumps->id.isol_loc,PETSC_COPY_VALUES,&is_petsc);CHKERRQ(ierr); /* to */ 639 ierr = VecScatterCreate(mumps->x_seq,is_iden,x,is_petsc,&mumps->scat_sol);CHKERRQ(ierr); 640 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 641 ierr = ISDestroy(&is_petsc);CHKERRQ(ierr); 642 643 mumps->ICNTL9_pre = mumps->id.ICNTL(9); /* save current value of id.ICNTL(9) */ 644 } 645 646 ierr = VecScatterBegin(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 647 ierr = VecScatterEnd(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 648 } 649 PetscFunctionReturn(0); 650 } 651 652 #undef __FUNCT__ 653 #define __FUNCT__ "MatSolveTranspose_MUMPS" 654 PetscErrorCode MatSolveTranspose_MUMPS(Mat A,Vec b,Vec x) 655 { 656 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 657 PetscErrorCode ierr; 658 659 PetscFunctionBegin; 660 mumps->id.ICNTL(9) = 0; 661 662 ierr = MatSolve_MUMPS(A,b,x);CHKERRQ(ierr); 663 664 mumps->id.ICNTL(9) = 1; 665 PetscFunctionReturn(0); 666 } 667 668 #undef __FUNCT__ 669 #define __FUNCT__ "MatMatSolve_MUMPS" 670 PetscErrorCode MatMatSolve_MUMPS(Mat A,Mat B,Mat X) 671 { 672 PetscErrorCode ierr; 673 PetscBool flg; 674 675 PetscFunctionBegin; 676 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 677 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 678 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 679 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 680 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_MUMPS() is not implemented yet"); 681 PetscFunctionReturn(0); 682 } 683 684 #if !defined(PETSC_USE_COMPLEX) 685 /* 686 input: 687 F: numeric factor 688 output: 689 nneg: total number of negative pivots 690 nzero: 0 691 npos: (global dimension of F) - nneg 692 */ 693 694 #undef __FUNCT__ 695 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 696 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 697 { 698 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 699 PetscErrorCode ierr; 700 PetscMPIInt size; 701 702 PetscFunctionBegin; 703 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);CHKERRQ(ierr); 704 /* MUMPS 4.3.1 calls ScaLAPACK when ICNTL(13)=0 (default), which does not offer the possibility to compute the inertia of a dense matrix. Set ICNTL(13)=1 to skip ScaLAPACK */ 705 if (size > 1 && mumps->id.ICNTL(13) != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"ICNTL(13)=%d. -mat_mumps_icntl_13 must be set as 1 for correct global matrix inertia\n",mumps->id.INFOG(13)); 706 707 if (nneg) *nneg = mumps->id.INFOG(12); 708 if (nzero || npos) { 709 if (mumps->id.ICNTL(24) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"-mat_mumps_icntl_24 must be set as 1 for null pivot row detection"); 710 if (nzero) *nzero = mumps->id.INFOG(28); 711 if (npos) *npos = F->rmap->N - (mumps->id.INFOG(12) + mumps->id.INFOG(28)); 712 } 713 PetscFunctionReturn(0); 714 } 715 #endif /* !defined(PETSC_USE_COMPLEX) */ 716 717 #undef __FUNCT__ 718 #define __FUNCT__ "MatFactorNumeric_MUMPS" 719 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 720 { 721 Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr; 722 PetscErrorCode ierr; 723 Mat F_diag; 724 PetscBool isMPIAIJ; 725 726 PetscFunctionBegin; 727 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 728 729 /* numerical factorization phase */ 730 /*-------------------------------*/ 731 mumps->id.job = JOB_FACTNUMERIC; 732 if (!mumps->id.ICNTL(18)) { 733 if (!mumps->myid) { 734 #if defined(PETSC_USE_COMPLEX) 735 #if defined(PETSC_USE_REAL_SINGLE) 736 mumps->id.a = (mumps_complex*)mumps->val; 737 #else 738 mumps->id.a = (mumps_double_complex*)mumps->val; 739 #endif 740 #else 741 mumps->id.a = mumps->val; 742 #endif 743 } 744 } else { 745 #if defined(PETSC_USE_COMPLEX) 746 #if defined(PETSC_USE_REAL_SINGLE) 747 mumps->id.a_loc = (mumps_complex*)mumps->val; 748 #else 749 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 750 #endif 751 #else 752 mumps->id.a_loc = mumps->val; 753 #endif 754 } 755 PetscMUMPS_c(&mumps->id); 756 if (mumps->id.INFOG(1) < 0) { 757 if (mumps->id.INFO(1) == -13) { 758 if (mumps->id.INFO(2) < 0) { 759 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d megabytes\n",-mumps->id.INFO(2)); 760 } else { 761 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d bytes\n",mumps->id.INFO(2)); 762 } 763 } else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: INFO(1)=%d, INFO(2)=%d\n",mumps->id.INFO(1),mumps->id.INFO(2)); 764 } 765 if (!mumps->myid && mumps->id.ICNTL(16) > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB," mumps->id.ICNTL(16):=%d\n",mumps->id.INFOG(16)); 766 767 (F)->assembled = PETSC_TRUE; 768 mumps->matstruc = SAME_NONZERO_PATTERN; 769 mumps->CleanUpMUMPS = PETSC_TRUE; 770 771 if (mumps->size > 1) { 772 PetscInt lsol_loc; 773 PetscScalar *sol_loc; 774 775 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 776 if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A; 777 else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A; 778 F_diag->assembled = PETSC_TRUE; 779 780 /* distributed solution; Create x_seq=sol_loc for repeated use */ 781 if (mumps->x_seq) { 782 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 783 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 784 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 785 } 786 lsol_loc = mumps->id.INFO(23); /* length of sol_loc */ 787 ierr = PetscMalloc2(lsol_loc,&sol_loc,lsol_loc,&mumps->id.isol_loc);CHKERRQ(ierr); 788 mumps->id.lsol_loc = lsol_loc; 789 #if defined(PETSC_USE_COMPLEX) 790 #if defined(PETSC_USE_REAL_SINGLE) 791 mumps->id.sol_loc = (mumps_complex*)sol_loc; 792 #else 793 mumps->id.sol_loc = (mumps_double_complex*)sol_loc; 794 #endif 795 #else 796 mumps->id.sol_loc = sol_loc; 797 #endif 798 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 799 } 800 PetscFunctionReturn(0); 801 } 802 803 /* Sets MUMPS options from the options database */ 804 #undef __FUNCT__ 805 #define __FUNCT__ "PetscSetMUMPSFromOptions" 806 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 807 { 808 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 809 PetscErrorCode ierr; 810 PetscInt icntl; 811 PetscBool flg; 812 813 PetscFunctionBegin; 814 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 815 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 816 if (flg) mumps->id.ICNTL(1) = icntl; 817 ierr = PetscOptionsInt("-mat_mumps_icntl_2","ICNTL(2): output stream for diagnostic printing, statistics, and warning","None",mumps->id.ICNTL(2),&icntl,&flg);CHKERRQ(ierr); 818 if (flg) mumps->id.ICNTL(2) = icntl; 819 ierr = PetscOptionsInt("-mat_mumps_icntl_3","ICNTL(3): output stream for global information, collected on the host","None",mumps->id.ICNTL(3),&icntl,&flg);CHKERRQ(ierr); 820 if (flg) mumps->id.ICNTL(3) = icntl; 821 822 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 823 if (flg) mumps->id.ICNTL(4) = icntl; 824 if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 825 826 ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): permuting and/or scaling the matrix (0 to 7)","None",mumps->id.ICNTL(6),&icntl,&flg);CHKERRQ(ierr); 827 if (flg) mumps->id.ICNTL(6) = icntl; 828 829 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7). 3=Scotch, 4=PORD, 5=Metis","None",mumps->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 830 if (flg) { 831 if (icntl== 1 && mumps->size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n"); 832 else mumps->id.ICNTL(7) = icntl; 833 } 834 835 ierr = PetscOptionsInt("-mat_mumps_icntl_8","ICNTL(8): scaling strategy (-2 to 8 or 77)","None",mumps->id.ICNTL(8),&mumps->id.ICNTL(8),NULL);CHKERRQ(ierr); 836 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);CHKERRQ(ierr); 837 ierr = PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): statistics related to the linear system solved (via -ksp_view)","None",mumps->id.ICNTL(11),&mumps->id.ICNTL(11),NULL);CHKERRQ(ierr); 838 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control: defines the ordering strategy with scaling constraints (0 to 3)","None",mumps->id.ICNTL(12),&mumps->id.ICNTL(12),NULL);CHKERRQ(ierr); 839 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control: with or without ScaLAPACK","None",mumps->id.ICNTL(13),&mumps->id.ICNTL(13),NULL);CHKERRQ(ierr); 840 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): percentage of estimated workspace increase","None",mumps->id.ICNTL(14),&mumps->id.ICNTL(14),NULL);CHKERRQ(ierr); 841 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);CHKERRQ(ierr); 842 843 ierr = PetscOptionsInt("-mat_mumps_icntl_22","ICNTL(22): in-core/out-of-core facility (0 or 1)","None",mumps->id.ICNTL(22),&mumps->id.ICNTL(22),NULL);CHKERRQ(ierr); 844 ierr = PetscOptionsInt("-mat_mumps_icntl_23","ICNTL(23): max size of the working memory (MB) that can allocate per processor","None",mumps->id.ICNTL(23),&mumps->id.ICNTL(23),NULL);CHKERRQ(ierr); 845 ierr = PetscOptionsInt("-mat_mumps_icntl_24","ICNTL(24): detection of null pivot rows (0 or 1)","None",mumps->id.ICNTL(24),&mumps->id.ICNTL(24),NULL);CHKERRQ(ierr); 846 if (mumps->id.ICNTL(24)) { 847 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 848 } 849 850 ierr = PetscOptionsInt("-mat_mumps_icntl_25","ICNTL(25): computation of a null space basis","None",mumps->id.ICNTL(25),&mumps->id.ICNTL(25),NULL);CHKERRQ(ierr); 851 ierr = PetscOptionsInt("-mat_mumps_icntl_26","ICNTL(26): Schur options for right-hand side or solution vector","None",mumps->id.ICNTL(26),&mumps->id.ICNTL(26),NULL);CHKERRQ(ierr); 852 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);CHKERRQ(ierr); 853 ierr = PetscOptionsInt("-mat_mumps_icntl_28","ICNTL(28): use 1 for sequential analysis and ictnl(7) ordering, or 2 for parallel analysis and ictnl(29) ordering","None",mumps->id.ICNTL(28),&mumps->id.ICNTL(28),NULL);CHKERRQ(ierr); 854 ierr = PetscOptionsInt("-mat_mumps_icntl_29","ICNTL(29): parallel ordering 1 = ptscotch 2 = parmetis","None",mumps->id.ICNTL(29),&mumps->id.ICNTL(29),NULL);CHKERRQ(ierr); 855 ierr = PetscOptionsInt("-mat_mumps_icntl_30","ICNTL(30): compute user-specified set of entries in inv(A)","None",mumps->id.ICNTL(30),&mumps->id.ICNTL(30),NULL);CHKERRQ(ierr); 856 ierr = PetscOptionsInt("-mat_mumps_icntl_31","ICNTL(31): factors can be discarded in the solve phase","None",mumps->id.ICNTL(31),&mumps->id.ICNTL(31),NULL);CHKERRQ(ierr); 857 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);CHKERRQ(ierr); 858 859 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);CHKERRQ(ierr); 860 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);CHKERRQ(ierr); 861 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);CHKERRQ(ierr); 862 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);CHKERRQ(ierr); 863 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);CHKERRQ(ierr); 864 865 ierr = PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL); 866 PetscOptionsEnd(); 867 PetscFunctionReturn(0); 868 } 869 870 #undef __FUNCT__ 871 #define __FUNCT__ "PetscInitializeMUMPS" 872 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps) 873 { 874 PetscErrorCode ierr; 875 876 PetscFunctionBegin; 877 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid); 878 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);CHKERRQ(ierr); 879 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));CHKERRQ(ierr); 880 881 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 882 883 mumps->id.job = JOB_INIT; 884 mumps->id.par = 1; /* host participates factorizaton and solve */ 885 mumps->id.sym = mumps->sym; 886 PetscMUMPS_c(&mumps->id); 887 888 mumps->CleanUpMUMPS = PETSC_FALSE; 889 mumps->scat_rhs = NULL; 890 mumps->scat_sol = NULL; 891 892 /* set PETSc-MUMPS default options - override MUMPS default */ 893 mumps->id.ICNTL(3) = 0; 894 mumps->id.ICNTL(4) = 0; 895 if (mumps->size == 1) { 896 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 897 } else { 898 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 899 mumps->id.ICNTL(21) = 1; /* distributed solution */ 900 } 901 902 /* schur */ 903 mumps->id.size_schur = 0; 904 mumps->id.listvar_schur = NULL; 905 mumps->id.schur = NULL; 906 PetscFunctionReturn(0); 907 } 908 909 /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */ 910 #undef __FUNCT__ 911 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 912 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 913 { 914 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 915 PetscErrorCode ierr; 916 Vec b; 917 IS is_iden; 918 const PetscInt M = A->rmap->N; 919 920 PetscFunctionBegin; 921 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 922 923 /* Set MUMPS options from the options database */ 924 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 925 926 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 927 928 /* analysis phase */ 929 /*----------------*/ 930 mumps->id.job = JOB_FACTSYMBOLIC; 931 mumps->id.n = M; 932 switch (mumps->id.ICNTL(18)) { 933 case 0: /* centralized assembled matrix input */ 934 if (!mumps->myid) { 935 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 936 if (mumps->id.ICNTL(6)>1) { 937 #if defined(PETSC_USE_COMPLEX) 938 #if defined(PETSC_USE_REAL_SINGLE) 939 mumps->id.a = (mumps_complex*)mumps->val; 940 #else 941 mumps->id.a = (mumps_double_complex*)mumps->val; 942 #endif 943 #else 944 mumps->id.a = mumps->val; 945 #endif 946 } 947 if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */ 948 /* 949 PetscBool flag; 950 ierr = ISEqual(r,c,&flag);CHKERRQ(ierr); 951 if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm"); 952 ierr = ISView(r,PETSC_VIEWER_STDOUT_SELF); 953 */ 954 if (!mumps->myid) { 955 const PetscInt *idx; 956 PetscInt i,*perm_in; 957 958 ierr = PetscMalloc1(M,&perm_in);CHKERRQ(ierr); 959 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 960 961 mumps->id.perm_in = perm_in; 962 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 963 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 964 } 965 } 966 } 967 break; 968 case 3: /* distributed assembled matrix input (size>1) */ 969 mumps->id.nz_loc = mumps->nz; 970 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 971 if (mumps->id.ICNTL(6)>1) { 972 #if defined(PETSC_USE_COMPLEX) 973 #if defined(PETSC_USE_REAL_SINGLE) 974 mumps->id.a_loc = (mumps_complex*)mumps->val; 975 #else 976 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 977 #endif 978 #else 979 mumps->id.a_loc = mumps->val; 980 #endif 981 } 982 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 983 if (!mumps->myid) { 984 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 985 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 986 } else { 987 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 988 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 989 } 990 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 991 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 992 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 993 ierr = VecDestroy(&b);CHKERRQ(ierr); 994 break; 995 } 996 PetscMUMPS_c(&mumps->id); 997 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 998 999 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1000 F->ops->solve = MatSolve_MUMPS; 1001 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1002 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 1003 PetscFunctionReturn(0); 1004 } 1005 1006 /* Note the Petsc r and c permutations are ignored */ 1007 #undef __FUNCT__ 1008 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 1009 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1010 { 1011 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1012 PetscErrorCode ierr; 1013 Vec b; 1014 IS is_iden; 1015 const PetscInt M = A->rmap->N; 1016 1017 PetscFunctionBegin; 1018 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1019 1020 /* Set MUMPS options from the options database */ 1021 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1022 1023 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1024 1025 /* analysis phase */ 1026 /*----------------*/ 1027 mumps->id.job = JOB_FACTSYMBOLIC; 1028 mumps->id.n = M; 1029 switch (mumps->id.ICNTL(18)) { 1030 case 0: /* centralized assembled matrix input */ 1031 if (!mumps->myid) { 1032 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1033 if (mumps->id.ICNTL(6)>1) { 1034 #if defined(PETSC_USE_COMPLEX) 1035 #if defined(PETSC_USE_REAL_SINGLE) 1036 mumps->id.a = (mumps_complex*)mumps->val; 1037 #else 1038 mumps->id.a = (mumps_double_complex*)mumps->val; 1039 #endif 1040 #else 1041 mumps->id.a = mumps->val; 1042 #endif 1043 } 1044 } 1045 break; 1046 case 3: /* distributed assembled matrix input (size>1) */ 1047 mumps->id.nz_loc = mumps->nz; 1048 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1049 if (mumps->id.ICNTL(6)>1) { 1050 #if defined(PETSC_USE_COMPLEX) 1051 #if defined(PETSC_USE_REAL_SINGLE) 1052 mumps->id.a_loc = (mumps_complex*)mumps->val; 1053 #else 1054 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1055 #endif 1056 #else 1057 mumps->id.a_loc = mumps->val; 1058 #endif 1059 } 1060 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1061 if (!mumps->myid) { 1062 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1063 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1064 } else { 1065 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1066 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1067 } 1068 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1069 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1070 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1071 ierr = VecDestroy(&b);CHKERRQ(ierr); 1072 break; 1073 } 1074 PetscMUMPS_c(&mumps->id); 1075 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 1076 1077 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1078 F->ops->solve = MatSolve_MUMPS; 1079 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1080 PetscFunctionReturn(0); 1081 } 1082 1083 /* Note the Petsc r permutation and factor info are ignored */ 1084 #undef __FUNCT__ 1085 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 1086 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 1087 { 1088 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1089 PetscErrorCode ierr; 1090 Vec b; 1091 IS is_iden; 1092 const PetscInt M = A->rmap->N; 1093 1094 PetscFunctionBegin; 1095 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1096 1097 /* Set MUMPS options from the options database */ 1098 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1099 1100 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1101 1102 /* analysis phase */ 1103 /*----------------*/ 1104 mumps->id.job = JOB_FACTSYMBOLIC; 1105 mumps->id.n = M; 1106 switch (mumps->id.ICNTL(18)) { 1107 case 0: /* centralized assembled matrix input */ 1108 if (!mumps->myid) { 1109 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1110 if (mumps->id.ICNTL(6)>1) { 1111 #if defined(PETSC_USE_COMPLEX) 1112 #if defined(PETSC_USE_REAL_SINGLE) 1113 mumps->id.a = (mumps_complex*)mumps->val; 1114 #else 1115 mumps->id.a = (mumps_double_complex*)mumps->val; 1116 #endif 1117 #else 1118 mumps->id.a = mumps->val; 1119 #endif 1120 } 1121 } 1122 break; 1123 case 3: /* distributed assembled matrix input (size>1) */ 1124 mumps->id.nz_loc = mumps->nz; 1125 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1126 if (mumps->id.ICNTL(6)>1) { 1127 #if defined(PETSC_USE_COMPLEX) 1128 #if defined(PETSC_USE_REAL_SINGLE) 1129 mumps->id.a_loc = (mumps_complex*)mumps->val; 1130 #else 1131 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1132 #endif 1133 #else 1134 mumps->id.a_loc = mumps->val; 1135 #endif 1136 } 1137 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1138 if (!mumps->myid) { 1139 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1140 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1141 } else { 1142 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1143 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1144 } 1145 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1146 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1147 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1148 ierr = VecDestroy(&b);CHKERRQ(ierr); 1149 break; 1150 } 1151 PetscMUMPS_c(&mumps->id); 1152 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 1153 1154 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1155 F->ops->solve = MatSolve_MUMPS; 1156 F->ops->solvetranspose = MatSolve_MUMPS; 1157 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 1158 #if !defined(PETSC_USE_COMPLEX) 1159 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1160 #else 1161 F->ops->getinertia = NULL; 1162 #endif 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "MatView_MUMPS" 1168 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1169 { 1170 PetscErrorCode ierr; 1171 PetscBool iascii; 1172 PetscViewerFormat format; 1173 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1174 1175 PetscFunctionBegin; 1176 /* check if matrix is mumps type */ 1177 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1178 1179 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1180 if (iascii) { 1181 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1182 if (format == PETSC_VIEWER_ASCII_INFO) { 1183 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1184 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1185 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1186 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));CHKERRQ(ierr); 1187 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));CHKERRQ(ierr); 1188 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));CHKERRQ(ierr); 1189 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));CHKERRQ(ierr); 1190 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));CHKERRQ(ierr); 1191 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));CHKERRQ(ierr); 1192 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));CHKERRQ(ierr); 1193 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));CHKERRQ(ierr); 1194 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));CHKERRQ(ierr); 1195 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));CHKERRQ(ierr); 1196 if (mumps->id.ICNTL(11)>0) { 1197 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));CHKERRQ(ierr); 1198 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));CHKERRQ(ierr); 1199 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));CHKERRQ(ierr); 1200 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));CHKERRQ(ierr); 1201 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));CHKERRQ(ierr); 1202 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));CHKERRQ(ierr); 1203 } 1204 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));CHKERRQ(ierr); 1205 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));CHKERRQ(ierr); 1206 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));CHKERRQ(ierr); 1207 /* ICNTL(15-17) not used */ 1208 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));CHKERRQ(ierr); 1209 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));CHKERRQ(ierr); 1210 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));CHKERRQ(ierr); 1211 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));CHKERRQ(ierr); 1212 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));CHKERRQ(ierr); 1213 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));CHKERRQ(ierr); 1214 1215 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));CHKERRQ(ierr); 1216 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));CHKERRQ(ierr); 1217 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));CHKERRQ(ierr); 1218 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));CHKERRQ(ierr); 1219 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));CHKERRQ(ierr); 1220 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));CHKERRQ(ierr); 1221 1222 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));CHKERRQ(ierr); 1223 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));CHKERRQ(ierr); 1224 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));CHKERRQ(ierr); 1225 1226 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));CHKERRQ(ierr); 1227 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));CHKERRQ(ierr); 1228 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));CHKERRQ(ierr); 1229 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));CHKERRQ(ierr); 1230 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));CHKERRQ(ierr); 1231 1232 /* infomation local to each processor */ 1233 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1234 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1235 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));CHKERRQ(ierr); 1236 ierr = PetscViewerFlush(viewer); 1237 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1238 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));CHKERRQ(ierr); 1239 ierr = PetscViewerFlush(viewer); 1240 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1241 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));CHKERRQ(ierr); 1242 ierr = PetscViewerFlush(viewer); 1243 1244 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1245 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));CHKERRQ(ierr); 1246 ierr = PetscViewerFlush(viewer); 1247 1248 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1249 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));CHKERRQ(ierr); 1250 ierr = PetscViewerFlush(viewer); 1251 1252 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1253 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));CHKERRQ(ierr); 1254 ierr = PetscViewerFlush(viewer); 1255 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1256 1257 if (!mumps->myid) { /* information from the host */ 1258 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));CHKERRQ(ierr); 1259 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));CHKERRQ(ierr); 1260 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));CHKERRQ(ierr); 1261 ierr = PetscViewerASCIIPrintf(viewer," (RINFOG(12) RINFOG(13))*2^INFOG(34) (determinant): (%g,%g)*(2^%d)\n",mumps->id.RINFOG(12),mumps->id.RINFOG(13),mumps->id.INFOG(34));CHKERRQ(ierr); 1262 1263 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));CHKERRQ(ierr); 1264 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));CHKERRQ(ierr); 1265 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));CHKERRQ(ierr); 1266 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));CHKERRQ(ierr); 1267 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));CHKERRQ(ierr); 1268 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));CHKERRQ(ierr); 1269 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));CHKERRQ(ierr); 1270 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));CHKERRQ(ierr); 1271 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));CHKERRQ(ierr); 1272 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));CHKERRQ(ierr); 1273 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));CHKERRQ(ierr); 1274 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));CHKERRQ(ierr); 1275 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));CHKERRQ(ierr); 1276 ierr = PetscViewerASCIIPrintf(viewer," INFOG(16) (estimated size (in MB) of all MUMPS internal data for factorization after analysis: value on the most memory consuming processor): %d \n",mumps->id.INFOG(16));CHKERRQ(ierr); 1277 ierr = PetscViewerASCIIPrintf(viewer," INFOG(17) (estimated size of all MUMPS internal data for factorization after analysis: sum over all processors): %d \n",mumps->id.INFOG(17));CHKERRQ(ierr); 1278 ierr = PetscViewerASCIIPrintf(viewer," INFOG(18) (size of all MUMPS internal data allocated during factorization: value on the most memory consuming processor): %d \n",mumps->id.INFOG(18));CHKERRQ(ierr); 1279 ierr = PetscViewerASCIIPrintf(viewer," INFOG(19) (size of all MUMPS internal data allocated during factorization: sum over all processors): %d \n",mumps->id.INFOG(19));CHKERRQ(ierr); 1280 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));CHKERRQ(ierr); 1281 ierr = PetscViewerASCIIPrintf(viewer," INFOG(21) (size in MB of memory effectively used during factorization - value on the most memory consuming processor): %d \n",mumps->id.INFOG(21));CHKERRQ(ierr); 1282 ierr = PetscViewerASCIIPrintf(viewer," INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): %d \n",mumps->id.INFOG(22));CHKERRQ(ierr); 1283 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));CHKERRQ(ierr); 1284 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));CHKERRQ(ierr); 1285 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));CHKERRQ(ierr); 1286 ierr = PetscViewerASCIIPrintf(viewer," INFOG(28) (after factorization: number of null pivots encountered): %d\n",mumps->id.INFOG(28));CHKERRQ(ierr); 1287 ierr = PetscViewerASCIIPrintf(viewer," INFOG(29) (after factorization: effective number of entries in the factors (sum over all processors)): %d\n",mumps->id.INFOG(29));CHKERRQ(ierr); 1288 ierr = PetscViewerASCIIPrintf(viewer," INFOG(30, 31) (after solution: size in Mbytes of memory used during solution phase): %d, %d\n",mumps->id.INFOG(30),mumps->id.INFOG(31));CHKERRQ(ierr); 1289 ierr = PetscViewerASCIIPrintf(viewer," INFOG(32) (after analysis: type of analysis done): %d\n",mumps->id.INFOG(32));CHKERRQ(ierr); 1290 ierr = PetscViewerASCIIPrintf(viewer," INFOG(33) (value used for ICNTL(8)): %d\n",mumps->id.INFOG(33));CHKERRQ(ierr); 1291 ierr = PetscViewerASCIIPrintf(viewer," INFOG(34) (exponent of the determinant if determinant is requested): %d\n",mumps->id.INFOG(34));CHKERRQ(ierr); 1292 } 1293 } 1294 } 1295 PetscFunctionReturn(0); 1296 } 1297 1298 #undef __FUNCT__ 1299 #define __FUNCT__ "MatGetInfo_MUMPS" 1300 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1301 { 1302 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1303 1304 PetscFunctionBegin; 1305 info->block_size = 1.0; 1306 info->nz_allocated = mumps->id.INFOG(20); 1307 info->nz_used = mumps->id.INFOG(20); 1308 info->nz_unneeded = 0.0; 1309 info->assemblies = 0.0; 1310 info->mallocs = 0.0; 1311 info->memory = 0.0; 1312 info->fill_ratio_given = 0; 1313 info->fill_ratio_needed = 0; 1314 info->factor_mallocs = 0; 1315 PetscFunctionReturn(0); 1316 } 1317 1318 /* -------------------------------------------------------------------------------------------*/ 1319 #undef __FUNCT__ 1320 #define __FUNCT__ "MatMumpsSetSchurIndices_MUMPS" 1321 PetscErrorCode MatMumpsSetSchurIndices_MUMPS(Mat F,PetscInt size,PetscInt idxs[]) 1322 { 1323 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1324 PetscErrorCode ierr; 1325 1326 PetscFunctionBegin; 1327 if (mumps->size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MUMPS parallel computation of Schur complements not yet supported from PETSc\n"); 1328 if (mumps->id.size_schur != size) { 1329 ierr = PetscFree2(mumps->id.listvar_schur,mumps->id.schur);CHKERRQ(ierr); 1330 mumps->id.size_schur = size; 1331 mumps->id.schur_lld = size; 1332 ierr = PetscMalloc2(size,&mumps->id.listvar_schur,size*size,&mumps->id.schur);CHKERRQ(ierr); 1333 } 1334 ierr = PetscMemcpy(mumps->id.listvar_schur,idxs,size*sizeof(PetscInt));CHKERRQ(ierr); 1335 if (F->factortype == MAT_FACTOR_LU) { 1336 mumps->id.ICNTL(19) = 3; /* return full matrix */ 1337 } else { 1338 mumps->id.ICNTL(19) = 2; /* return lower triangular part */ 1339 } 1340 PetscFunctionReturn(0); 1341 } 1342 1343 #undef __FUNCT__ 1344 #define __FUNCT__ "MatMumpsSetSchurIndices" 1345 /*@ 1346 MatMumpsSetSchurIndices - Set indices defining the Schur complement that MUMPS will compute during the factorization steps 1347 1348 Logically Collective on Mat 1349 1350 Input Parameters: 1351 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1352 . size - size of the Schur complement indices 1353 - idxs[] - array of Schur complement indices 1354 1355 Notes: 1356 The user has to free the array idxs[] since it is copied by the routine. 1357 Currently implemented for sequential matrices 1358 1359 Level: advanced 1360 1361 References: MUMPS Users' Guide 1362 1363 .seealso: MatGetFactor() 1364 @*/ 1365 PetscErrorCode MatMumpsSetSchurIndices(Mat F,PetscInt size,PetscInt idxs[]) 1366 { 1367 PetscErrorCode ierr; 1368 1369 PetscFunctionBegin; 1370 PetscValidIntPointer(idxs,3); 1371 ierr = PetscTryMethod(F,"MatMumpsSetSchurIndices_C",(Mat,PetscInt,PetscInt[]),(F,size,idxs));CHKERRQ(ierr); 1372 PetscFunctionReturn(0); 1373 } 1374 1375 /* -------------------------------------------------------------------------------------------*/ 1376 #undef __FUNCT__ 1377 #define __FUNCT__ "MatMumpsGetSchurComplement_MUMPS" 1378 PetscErrorCode MatMumpsGetSchurComplement_MUMPS(Mat F,Mat* S) 1379 { 1380 Mat St; 1381 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1382 PetscScalar *array; 1383 PetscErrorCode ierr; 1384 1385 PetscFunctionBegin; 1386 if (mumps->id.job != JOB_FACTNUMERIC) { 1387 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_WRONG,"Numerical factorization phase not yet performed! You should call MatFactorSymbolic/Numeric before"); 1388 } else if (mumps->id.size_schur == 0) { 1389 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_WRONG,"Schur indices not set! You should call MatMumpsSetSchurIndices before"); 1390 } 1391 1392 ierr = MatCreate(PetscObjectComm((PetscObject)F),&St);CHKERRQ(ierr); 1393 ierr = MatSetSizes(St,PETSC_DECIDE,PETSC_DECIDE,mumps->id.size_schur,mumps->id.size_schur);CHKERRQ(ierr); 1394 ierr = MatSetType(St,MATDENSE);CHKERRQ(ierr); 1395 ierr = MatSetUp(St);CHKERRQ(ierr); 1396 ierr = MatDenseGetArray(St,&array);CHKERRQ(ierr); 1397 if (mumps->sym == 0) { /* MUMPS always return a full matrix */ 1398 if (mumps->id.ICNTL(19) == 1) { /* stored by rows */ 1399 PetscInt i,j,N=mumps->id.size_schur; 1400 for (i=0;i<N;i++) { 1401 for (j=0;j<N;j++) { 1402 PetscScalar val = mumps->id.schur[i*N+j]; 1403 array[j*N+i] = val; 1404 } 1405 } 1406 } else { /* stored by columns */ 1407 ierr = PetscMemcpy(array,mumps->id.schur,mumps->id.size_schur*mumps->id.size_schur*sizeof(PetscScalar));CHKERRQ(ierr); 1408 } 1409 } else { /* either full or lower-triangular (not packed) */ 1410 if (mumps->id.ICNTL(19) == 2) { /* lower triangular stored by columns */ 1411 PetscInt i,j,N=mumps->id.size_schur; 1412 for (i=0;i<N;i++) { 1413 for (j=i;j<N;j++) { 1414 PetscScalar val = mumps->id.schur[i*N+j]; 1415 array[i*N+j] = val; 1416 } 1417 for (j=i;j<N;j++) { 1418 PetscScalar val = mumps->id.schur[i*N+j]; 1419 array[j*N+i] = val; 1420 } 1421 } 1422 } else if (mumps->id.ICNTL(19) == 3) { /* full matrix */ 1423 ierr = PetscMemcpy(array,mumps->id.schur,mumps->id.size_schur*mumps->id.size_schur*sizeof(PetscScalar));CHKERRQ(ierr); 1424 } else { /* ICNTL(19) == 1 lower triangular stored by rows */ 1425 PetscInt i,j,N=mumps->id.size_schur; 1426 for (i=0;i<N;i++) { 1427 for (j=0;j<i+1;j++) { 1428 PetscScalar val = mumps->id.schur[i*N+j]; 1429 array[i*N+j] = val; 1430 } 1431 for (j=0;j<i+1;j++) { 1432 PetscScalar val = mumps->id.schur[i*N+j]; 1433 array[j*N+i] = val; 1434 } 1435 } 1436 } 1437 } 1438 ierr = MatDenseRestoreArray(St,&array);CHKERRQ(ierr); 1439 *S = St; 1440 PetscFunctionReturn(0); 1441 } 1442 1443 #undef __FUNCT__ 1444 #define __FUNCT__ "MatMumpsGetSchurComplement" 1445 /*@ 1446 MatMumpsGetSchurComplement - Get Schur complement matrix computed by MUMPS during the factorization step 1447 1448 Logically Collective on Mat 1449 1450 Input Parameters: 1451 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1452 . *S - location where to return the Schur complement (MATDENSE) 1453 1454 Notes: 1455 Currently implemented for sequential matrices 1456 1457 Level: advanced 1458 1459 References: MUMPS Users' Guide 1460 1461 .seealso: MatGetFactor() 1462 @*/ 1463 PetscErrorCode MatMumpsGetSchurComplement(Mat F,Mat* S) 1464 { 1465 PetscErrorCode ierr; 1466 1467 PetscFunctionBegin; 1468 ierr = PetscTryMethod(F,"MatMumpsGetSchurComplement_C",(Mat,Mat*),(F,S));CHKERRQ(ierr); 1469 PetscFunctionReturn(0); 1470 } 1471 1472 /* -------------------------------------------------------------------------------------------*/ 1473 #undef __FUNCT__ 1474 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1475 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1476 { 1477 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1478 1479 PetscFunctionBegin; 1480 mumps->id.ICNTL(icntl) = ival; 1481 PetscFunctionReturn(0); 1482 } 1483 1484 #undef __FUNCT__ 1485 #define __FUNCT__ "MatMumpsGetIcntl_MUMPS" 1486 PetscErrorCode MatMumpsGetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt *ival) 1487 { 1488 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1489 1490 PetscFunctionBegin; 1491 *ival = mumps->id.ICNTL(icntl); 1492 PetscFunctionReturn(0); 1493 } 1494 1495 #undef __FUNCT__ 1496 #define __FUNCT__ "MatMumpsSetIcntl" 1497 /*@ 1498 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1499 1500 Logically Collective on Mat 1501 1502 Input Parameters: 1503 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1504 . icntl - index of MUMPS parameter array ICNTL() 1505 - ival - value of MUMPS ICNTL(icntl) 1506 1507 Options Database: 1508 . -mat_mumps_icntl_<icntl> <ival> 1509 1510 Level: beginner 1511 1512 References: MUMPS Users' Guide 1513 1514 .seealso: MatGetFactor() 1515 @*/ 1516 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1517 { 1518 PetscErrorCode ierr; 1519 1520 PetscFunctionBegin; 1521 PetscValidLogicalCollectiveInt(F,icntl,2); 1522 PetscValidLogicalCollectiveInt(F,ival,3); 1523 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1524 PetscFunctionReturn(0); 1525 } 1526 1527 #undef __FUNCT__ 1528 #define __FUNCT__ "MatMumpsGetIcntl" 1529 /*@ 1530 MatMumpsGetIcntl - Get MUMPS parameter ICNTL() 1531 1532 Logically Collective on Mat 1533 1534 Input Parameters: 1535 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1536 - icntl - index of MUMPS parameter array ICNTL() 1537 1538 Output Parameter: 1539 . ival - value of MUMPS ICNTL(icntl) 1540 1541 Level: beginner 1542 1543 References: MUMPS Users' Guide 1544 1545 .seealso: MatGetFactor() 1546 @*/ 1547 PetscErrorCode MatMumpsGetIcntl(Mat F,PetscInt icntl,PetscInt *ival) 1548 { 1549 PetscErrorCode ierr; 1550 1551 PetscFunctionBegin; 1552 PetscValidLogicalCollectiveInt(F,icntl,2); 1553 PetscValidIntPointer(ival,3); 1554 ierr = PetscTryMethod(F,"MatMumpsGetIcntl_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1555 PetscFunctionReturn(0); 1556 } 1557 1558 /* -------------------------------------------------------------------------------------------*/ 1559 #undef __FUNCT__ 1560 #define __FUNCT__ "MatMumpsSetCntl_MUMPS" 1561 PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val) 1562 { 1563 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1564 1565 PetscFunctionBegin; 1566 mumps->id.CNTL(icntl) = val; 1567 PetscFunctionReturn(0); 1568 } 1569 1570 #undef __FUNCT__ 1571 #define __FUNCT__ "MatMumpsGetCntl_MUMPS" 1572 PetscErrorCode MatMumpsGetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal *val) 1573 { 1574 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1575 1576 PetscFunctionBegin; 1577 *val = mumps->id.CNTL(icntl); 1578 PetscFunctionReturn(0); 1579 } 1580 1581 #undef __FUNCT__ 1582 #define __FUNCT__ "MatMumpsSetCntl" 1583 /*@ 1584 MatMumpsSetCntl - Set MUMPS parameter CNTL() 1585 1586 Logically Collective on Mat 1587 1588 Input Parameters: 1589 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1590 . icntl - index of MUMPS parameter array CNTL() 1591 - val - value of MUMPS CNTL(icntl) 1592 1593 Options Database: 1594 . -mat_mumps_cntl_<icntl> <val> 1595 1596 Level: beginner 1597 1598 References: MUMPS Users' Guide 1599 1600 .seealso: MatGetFactor() 1601 @*/ 1602 PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val) 1603 { 1604 PetscErrorCode ierr; 1605 1606 PetscFunctionBegin; 1607 PetscValidLogicalCollectiveInt(F,icntl,2); 1608 PetscValidLogicalCollectiveReal(F,val,3); 1609 ierr = PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));CHKERRQ(ierr); 1610 PetscFunctionReturn(0); 1611 } 1612 1613 #undef __FUNCT__ 1614 #define __FUNCT__ "MatMumpsGetCntl" 1615 /*@ 1616 MatMumpsGetCntl - Get MUMPS parameter CNTL() 1617 1618 Logically Collective on Mat 1619 1620 Input Parameters: 1621 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1622 - icntl - index of MUMPS parameter array CNTL() 1623 1624 Output Parameter: 1625 . val - value of MUMPS CNTL(icntl) 1626 1627 Level: beginner 1628 1629 References: MUMPS Users' Guide 1630 1631 .seealso: MatGetFactor() 1632 @*/ 1633 PetscErrorCode MatMumpsGetCntl(Mat F,PetscInt icntl,PetscReal *val) 1634 { 1635 PetscErrorCode ierr; 1636 1637 PetscFunctionBegin; 1638 PetscValidLogicalCollectiveInt(F,icntl,2); 1639 PetscValidRealPointer(val,3); 1640 ierr = PetscTryMethod(F,"MatMumpsGetCntl_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1641 PetscFunctionReturn(0); 1642 } 1643 1644 #undef __FUNCT__ 1645 #define __FUNCT__ "MatMumpsGetInfo_MUMPS" 1646 PetscErrorCode MatMumpsGetInfo_MUMPS(Mat F,PetscInt icntl,PetscInt *info) 1647 { 1648 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1649 1650 PetscFunctionBegin; 1651 *info = mumps->id.INFO(icntl); 1652 PetscFunctionReturn(0); 1653 } 1654 1655 #undef __FUNCT__ 1656 #define __FUNCT__ "MatMumpsGetInfog_MUMPS" 1657 PetscErrorCode MatMumpsGetInfog_MUMPS(Mat F,PetscInt icntl,PetscInt *infog) 1658 { 1659 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1660 1661 PetscFunctionBegin; 1662 *infog = mumps->id.INFOG(icntl); 1663 PetscFunctionReturn(0); 1664 } 1665 1666 #undef __FUNCT__ 1667 #define __FUNCT__ "MatMumpsGetRinfo_MUMPS" 1668 PetscErrorCode MatMumpsGetRinfo_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfo) 1669 { 1670 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1671 1672 PetscFunctionBegin; 1673 *rinfo = mumps->id.RINFO(icntl); 1674 PetscFunctionReturn(0); 1675 } 1676 1677 #undef __FUNCT__ 1678 #define __FUNCT__ "MatMumpsGetRinfog_MUMPS" 1679 PetscErrorCode MatMumpsGetRinfog_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfog) 1680 { 1681 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1682 1683 PetscFunctionBegin; 1684 *rinfog = mumps->id.RINFOG(icntl); 1685 PetscFunctionReturn(0); 1686 } 1687 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "MatMumpsGetInfo" 1690 /*@ 1691 MatMumpsGetInfo - Get MUMPS parameter INFO() 1692 1693 Logically Collective on Mat 1694 1695 Input Parameters: 1696 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1697 - icntl - index of MUMPS parameter array INFO() 1698 1699 Output Parameter: 1700 . ival - value of MUMPS INFO(icntl) 1701 1702 Level: beginner 1703 1704 References: MUMPS Users' Guide 1705 1706 .seealso: MatGetFactor() 1707 @*/ 1708 PetscErrorCode MatMumpsGetInfo(Mat F,PetscInt icntl,PetscInt *ival) 1709 { 1710 PetscErrorCode ierr; 1711 1712 PetscFunctionBegin; 1713 PetscValidIntPointer(ival,3); 1714 ierr = PetscTryMethod(F,"MatMumpsGetInfo_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1715 PetscFunctionReturn(0); 1716 } 1717 1718 #undef __FUNCT__ 1719 #define __FUNCT__ "MatMumpsGetInfog" 1720 /*@ 1721 MatMumpsGetInfog - Get MUMPS parameter INFOG() 1722 1723 Logically Collective on Mat 1724 1725 Input Parameters: 1726 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1727 - icntl - index of MUMPS parameter array INFOG() 1728 1729 Output Parameter: 1730 . ival - value of MUMPS INFOG(icntl) 1731 1732 Level: beginner 1733 1734 References: MUMPS Users' Guide 1735 1736 .seealso: MatGetFactor() 1737 @*/ 1738 PetscErrorCode MatMumpsGetInfog(Mat F,PetscInt icntl,PetscInt *ival) 1739 { 1740 PetscErrorCode ierr; 1741 1742 PetscFunctionBegin; 1743 PetscValidIntPointer(ival,3); 1744 ierr = PetscTryMethod(F,"MatMumpsGetInfog_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 #undef __FUNCT__ 1749 #define __FUNCT__ "MatMumpsGetRinfo" 1750 /*@ 1751 MatMumpsGetRinfo - Get MUMPS parameter RINFO() 1752 1753 Logically Collective on Mat 1754 1755 Input Parameters: 1756 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1757 - icntl - index of MUMPS parameter array RINFO() 1758 1759 Output Parameter: 1760 . val - value of MUMPS RINFO(icntl) 1761 1762 Level: beginner 1763 1764 References: MUMPS Users' Guide 1765 1766 .seealso: MatGetFactor() 1767 @*/ 1768 PetscErrorCode MatMumpsGetRinfo(Mat F,PetscInt icntl,PetscReal *val) 1769 { 1770 PetscErrorCode ierr; 1771 1772 PetscFunctionBegin; 1773 PetscValidRealPointer(val,3); 1774 ierr = PetscTryMethod(F,"MatMumpsGetRinfo_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1775 PetscFunctionReturn(0); 1776 } 1777 1778 #undef __FUNCT__ 1779 #define __FUNCT__ "MatMumpsGetRinfog" 1780 /*@ 1781 MatMumpsGetRinfog - Get MUMPS parameter RINFOG() 1782 1783 Logically Collective on Mat 1784 1785 Input Parameters: 1786 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1787 - icntl - index of MUMPS parameter array RINFOG() 1788 1789 Output Parameter: 1790 . val - value of MUMPS RINFOG(icntl) 1791 1792 Level: beginner 1793 1794 References: MUMPS Users' Guide 1795 1796 .seealso: MatGetFactor() 1797 @*/ 1798 PetscErrorCode MatMumpsGetRinfog(Mat F,PetscInt icntl,PetscReal *val) 1799 { 1800 PetscErrorCode ierr; 1801 1802 PetscFunctionBegin; 1803 PetscValidRealPointer(val,3); 1804 ierr = PetscTryMethod(F,"MatMumpsGetRinfog_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1805 PetscFunctionReturn(0); 1806 } 1807 1808 /*MC 1809 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1810 distributed and sequential matrices via the external package MUMPS. 1811 1812 Works with MATAIJ and MATSBAIJ matrices 1813 1814 Options Database Keys: 1815 + -mat_mumps_icntl_4 <0,...,4> - print level 1816 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1817 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec) 1818 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1819 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1820 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1821 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1822 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1823 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1824 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1825 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1826 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1827 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1828 1829 Level: beginner 1830 1831 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1832 1833 M*/ 1834 1835 #undef __FUNCT__ 1836 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1837 static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1838 { 1839 PetscFunctionBegin; 1840 *type = MATSOLVERMUMPS; 1841 PetscFunctionReturn(0); 1842 } 1843 1844 /* MatGetFactor for Seq and MPI AIJ matrices */ 1845 #undef __FUNCT__ 1846 #define __FUNCT__ "MatGetFactor_aij_mumps" 1847 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1848 { 1849 Mat B; 1850 PetscErrorCode ierr; 1851 Mat_MUMPS *mumps; 1852 PetscBool isSeqAIJ; 1853 1854 PetscFunctionBegin; 1855 /* Create the factorization matrix */ 1856 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1857 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1858 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1859 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1860 if (isSeqAIJ) { 1861 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1862 } else { 1863 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1864 } 1865 1866 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1867 1868 B->ops->view = MatView_MUMPS; 1869 B->ops->getinfo = MatGetInfo_MUMPS; 1870 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1871 1872 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1873 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1874 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1875 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1876 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1877 1878 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1879 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1880 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1881 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1882 1883 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 1884 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 1885 if (ftype == MAT_FACTOR_LU) { 1886 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1887 B->factortype = MAT_FACTOR_LU; 1888 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1889 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1890 mumps->sym = 0; 1891 } else { 1892 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1893 B->factortype = MAT_FACTOR_CHOLESKY; 1894 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1895 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1896 if (A->spd_set && A->spd) mumps->sym = 1; 1897 else mumps->sym = 2; 1898 } 1899 1900 mumps->isAIJ = PETSC_TRUE; 1901 mumps->Destroy = B->ops->destroy; 1902 B->ops->destroy = MatDestroy_MUMPS; 1903 B->spptr = (void*)mumps; 1904 1905 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1906 1907 *F = B; 1908 PetscFunctionReturn(0); 1909 } 1910 1911 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1912 #undef __FUNCT__ 1913 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1914 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1915 { 1916 Mat B; 1917 PetscErrorCode ierr; 1918 Mat_MUMPS *mumps; 1919 PetscBool isSeqSBAIJ; 1920 1921 PetscFunctionBegin; 1922 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1923 if (A->rmap->bs > 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with block size > 1 with MUMPS Cholesky, use AIJ matrix instead"); 1924 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1925 /* Create the factorization matrix */ 1926 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1927 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1928 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1929 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1930 if (isSeqSBAIJ) { 1931 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 1932 1933 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1934 } else { 1935 ierr = MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);CHKERRQ(ierr); 1936 1937 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1938 } 1939 1940 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1941 B->ops->view = MatView_MUMPS; 1942 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1943 1944 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1945 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1946 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1947 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1948 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1949 1950 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1951 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1952 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1953 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1954 1955 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 1956 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 1957 1958 B->factortype = MAT_FACTOR_CHOLESKY; 1959 if (A->spd_set && A->spd) mumps->sym = 1; 1960 else mumps->sym = 2; 1961 1962 mumps->isAIJ = PETSC_FALSE; 1963 mumps->Destroy = B->ops->destroy; 1964 B->ops->destroy = MatDestroy_MUMPS; 1965 B->spptr = (void*)mumps; 1966 1967 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1968 1969 *F = B; 1970 PetscFunctionReturn(0); 1971 } 1972 1973 #undef __FUNCT__ 1974 #define __FUNCT__ "MatGetFactor_baij_mumps" 1975 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 1976 { 1977 Mat B; 1978 PetscErrorCode ierr; 1979 Mat_MUMPS *mumps; 1980 PetscBool isSeqBAIJ; 1981 1982 PetscFunctionBegin; 1983 /* Create the factorization matrix */ 1984 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1985 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1986 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1987 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1988 if (isSeqBAIJ) { 1989 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);CHKERRQ(ierr); 1990 } else { 1991 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 1992 } 1993 1994 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1995 if (ftype == MAT_FACTOR_LU) { 1996 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1997 B->factortype = MAT_FACTOR_LU; 1998 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 1999 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 2000 mumps->sym = 0; 2001 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 2002 2003 B->ops->view = MatView_MUMPS; 2004 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 2005 2006 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 2007 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 2008 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 2009 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 2010 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 2011 2012 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 2013 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 2014 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 2015 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 2016 2017 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 2018 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 2019 2020 mumps->isAIJ = PETSC_TRUE; 2021 mumps->Destroy = B->ops->destroy; 2022 B->ops->destroy = MatDestroy_MUMPS; 2023 B->spptr = (void*)mumps; 2024 2025 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 2026 2027 *F = B; 2028 PetscFunctionReturn(0); 2029 } 2030 2031 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 2032 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat,MatFactorType,Mat*); 2033 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 2034 2035 #undef __FUNCT__ 2036 #define __FUNCT__ "MatSolverPackageRegister_MUMPS" 2037 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MUMPS(void) 2038 { 2039 PetscErrorCode ierr; 2040 2041 PetscFunctionBegin; 2042 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2043 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2044 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2045 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2046 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPISBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2047 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2048 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2049 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2050 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2051 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2052 PetscFunctionReturn(0); 2053 } 2054 2055