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 553 mumps->id.job = JOB_END; 554 PetscMUMPS_c(&mumps->id); 555 ierr = MPI_Comm_free(&(mumps->comm_mumps));CHKERRQ(ierr); 556 } 557 if (mumps->Destroy) { 558 ierr = (mumps->Destroy)(A);CHKERRQ(ierr); 559 } 560 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 561 562 /* clear composed functions */ 563 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);CHKERRQ(ierr); 564 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetIcntl_C",NULL);CHKERRQ(ierr); 565 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetIcntl_C",NULL);CHKERRQ(ierr); 566 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetCntl_C",NULL);CHKERRQ(ierr); 567 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetCntl_C",NULL);CHKERRQ(ierr); 568 569 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfo_C",NULL);CHKERRQ(ierr); 570 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfog_C",NULL);CHKERRQ(ierr); 571 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfo_C",NULL);CHKERRQ(ierr); 572 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfog_C",NULL);CHKERRQ(ierr); 573 574 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetSchurIndices_C",NULL);CHKERRQ(ierr); 575 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetSchurComplement_C",NULL);CHKERRQ(ierr); 576 PetscFunctionReturn(0); 577 } 578 579 #undef __FUNCT__ 580 #define __FUNCT__ "MatSolve_MUMPS" 581 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 582 { 583 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 584 PetscScalar *array; 585 Vec b_seq; 586 IS is_iden,is_petsc; 587 PetscErrorCode ierr; 588 PetscInt i; 589 static PetscBool cite1 = PETSC_FALSE,cite2 = PETSC_FALSE; 590 591 PetscFunctionBegin; 592 if (mumps->id.size_schur) { 593 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use MatSolve when Schur complement has been requested\n"); 594 } 595 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); 596 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); 597 mumps->id.nrhs = 1; 598 b_seq = mumps->b_seq; 599 if (mumps->size > 1) { 600 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 601 ierr = VecScatterBegin(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 602 ierr = VecScatterEnd(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 603 if (!mumps->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 604 } else { /* size == 1 */ 605 ierr = VecCopy(b,x);CHKERRQ(ierr); 606 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 607 } 608 if (!mumps->myid) { /* define rhs on the host */ 609 mumps->id.nrhs = 1; 610 #if defined(PETSC_USE_COMPLEX) 611 #if defined(PETSC_USE_REAL_SINGLE) 612 mumps->id.rhs = (mumps_complex*)array; 613 #else 614 mumps->id.rhs = (mumps_double_complex*)array; 615 #endif 616 #else 617 mumps->id.rhs = array; 618 #endif 619 } 620 621 /* solve phase */ 622 /*-------------*/ 623 mumps->id.job = JOB_SOLVE; 624 PetscMUMPS_c(&mumps->id); 625 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)); 626 627 if (mumps->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 628 if (mumps->scat_sol && mumps->ICNTL9_pre != mumps->id.ICNTL(9)) { 629 /* when id.ICNTL(9) changes, the contents of lsol_loc may change (not its size, lsol_loc), recreates scat_sol */ 630 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 631 } 632 if (!mumps->scat_sol) { /* create scatter scat_sol */ 633 ierr = ISCreateStride(PETSC_COMM_SELF,mumps->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 634 for (i=0; i<mumps->id.lsol_loc; i++) { 635 mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 636 } 637 ierr = ISCreateGeneral(PETSC_COMM_SELF,mumps->id.lsol_loc,mumps->id.isol_loc,PETSC_COPY_VALUES,&is_petsc);CHKERRQ(ierr); /* to */ 638 ierr = VecScatterCreate(mumps->x_seq,is_iden,x,is_petsc,&mumps->scat_sol);CHKERRQ(ierr); 639 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 640 ierr = ISDestroy(&is_petsc);CHKERRQ(ierr); 641 642 mumps->ICNTL9_pre = mumps->id.ICNTL(9); /* save current value of id.ICNTL(9) */ 643 } 644 645 ierr = VecScatterBegin(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 646 ierr = VecScatterEnd(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 647 } 648 PetscFunctionReturn(0); 649 } 650 651 #undef __FUNCT__ 652 #define __FUNCT__ "MatSolveTranspose_MUMPS" 653 PetscErrorCode MatSolveTranspose_MUMPS(Mat A,Vec b,Vec x) 654 { 655 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 656 PetscErrorCode ierr; 657 658 PetscFunctionBegin; 659 mumps->id.ICNTL(9) = 0; 660 661 ierr = MatSolve_MUMPS(A,b,x);CHKERRQ(ierr); 662 663 mumps->id.ICNTL(9) = 1; 664 PetscFunctionReturn(0); 665 } 666 667 #undef __FUNCT__ 668 #define __FUNCT__ "MatMatSolve_MUMPS" 669 PetscErrorCode MatMatSolve_MUMPS(Mat A,Mat B,Mat X) 670 { 671 PetscErrorCode ierr; 672 PetscBool flg; 673 674 PetscFunctionBegin; 675 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 676 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 677 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 678 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 679 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_MUMPS() is not implemented yet"); 680 PetscFunctionReturn(0); 681 } 682 683 #if !defined(PETSC_USE_COMPLEX) 684 /* 685 input: 686 F: numeric factor 687 output: 688 nneg: total number of negative pivots 689 nzero: 0 690 npos: (global dimension of F) - nneg 691 */ 692 693 #undef __FUNCT__ 694 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 695 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 696 { 697 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 698 PetscErrorCode ierr; 699 PetscMPIInt size; 700 701 PetscFunctionBegin; 702 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);CHKERRQ(ierr); 703 /* 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 */ 704 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)); 705 706 if (nneg) *nneg = mumps->id.INFOG(12); 707 if (nzero || npos) { 708 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"); 709 if (nzero) *nzero = mumps->id.INFOG(28); 710 if (npos) *npos = F->rmap->N - (mumps->id.INFOG(12) + mumps->id.INFOG(28)); 711 } 712 PetscFunctionReturn(0); 713 } 714 #endif /* !defined(PETSC_USE_COMPLEX) */ 715 716 #undef __FUNCT__ 717 #define __FUNCT__ "MatFactorNumeric_MUMPS" 718 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 719 { 720 Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr; 721 PetscErrorCode ierr; 722 Mat F_diag; 723 PetscBool isMPIAIJ; 724 725 PetscFunctionBegin; 726 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 727 728 /* numerical factorization phase */ 729 /*-------------------------------*/ 730 mumps->id.job = JOB_FACTNUMERIC; 731 if (!mumps->id.ICNTL(18)) { 732 if (!mumps->myid) { 733 #if defined(PETSC_USE_COMPLEX) 734 #if defined(PETSC_USE_REAL_SINGLE) 735 mumps->id.a = (mumps_complex*)mumps->val; 736 #else 737 mumps->id.a = (mumps_double_complex*)mumps->val; 738 #endif 739 #else 740 mumps->id.a = mumps->val; 741 #endif 742 } 743 } else { 744 #if defined(PETSC_USE_COMPLEX) 745 #if defined(PETSC_USE_REAL_SINGLE) 746 mumps->id.a_loc = (mumps_complex*)mumps->val; 747 #else 748 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 749 #endif 750 #else 751 mumps->id.a_loc = mumps->val; 752 #endif 753 } 754 PetscMUMPS_c(&mumps->id); 755 if (mumps->id.INFOG(1) < 0) { 756 if (mumps->id.INFO(1) == -13) { 757 if (mumps->id.INFO(2) < 0) { 758 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)); 759 } else { 760 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)); 761 } 762 } 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)); 763 } 764 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)); 765 766 (F)->assembled = PETSC_TRUE; 767 mumps->matstruc = SAME_NONZERO_PATTERN; 768 mumps->CleanUpMUMPS = PETSC_TRUE; 769 770 if (mumps->size > 1) { 771 PetscInt lsol_loc; 772 PetscScalar *sol_loc; 773 774 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 775 if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A; 776 else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A; 777 F_diag->assembled = PETSC_TRUE; 778 779 /* distributed solution; Create x_seq=sol_loc for repeated use */ 780 if (mumps->x_seq) { 781 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 782 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 783 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 784 } 785 lsol_loc = mumps->id.INFO(23); /* length of sol_loc */ 786 ierr = PetscMalloc2(lsol_loc,&sol_loc,lsol_loc,&mumps->id.isol_loc);CHKERRQ(ierr); 787 mumps->id.lsol_loc = lsol_loc; 788 #if defined(PETSC_USE_COMPLEX) 789 #if defined(PETSC_USE_REAL_SINGLE) 790 mumps->id.sol_loc = (mumps_complex*)sol_loc; 791 #else 792 mumps->id.sol_loc = (mumps_double_complex*)sol_loc; 793 #endif 794 #else 795 mumps->id.sol_loc = sol_loc; 796 #endif 797 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 798 } 799 PetscFunctionReturn(0); 800 } 801 802 /* Sets MUMPS options from the options database */ 803 #undef __FUNCT__ 804 #define __FUNCT__ "PetscSetMUMPSFromOptions" 805 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 806 { 807 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 808 PetscErrorCode ierr; 809 PetscInt icntl; 810 PetscBool flg; 811 812 PetscFunctionBegin; 813 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 814 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 815 if (flg) mumps->id.ICNTL(1) = icntl; 816 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); 817 if (flg) mumps->id.ICNTL(2) = icntl; 818 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); 819 if (flg) mumps->id.ICNTL(3) = icntl; 820 821 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 822 if (flg) mumps->id.ICNTL(4) = icntl; 823 if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 824 825 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); 826 if (flg) mumps->id.ICNTL(6) = icntl; 827 828 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); 829 if (flg) { 830 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"); 831 else mumps->id.ICNTL(7) = icntl; 832 } 833 834 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); 835 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);CHKERRQ(ierr); 836 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); 837 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); 838 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); 839 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); 840 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);CHKERRQ(ierr); 841 842 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); 843 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); 844 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); 845 if (mumps->id.ICNTL(24)) { 846 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 847 } 848 849 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); 850 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); 851 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);CHKERRQ(ierr); 852 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); 853 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); 854 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); 855 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); 856 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);CHKERRQ(ierr); 857 858 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);CHKERRQ(ierr); 859 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);CHKERRQ(ierr); 860 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);CHKERRQ(ierr); 861 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);CHKERRQ(ierr); 862 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);CHKERRQ(ierr); 863 864 ierr = PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL); 865 PetscOptionsEnd(); 866 PetscFunctionReturn(0); 867 } 868 869 #undef __FUNCT__ 870 #define __FUNCT__ "PetscInitializeMUMPS" 871 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps) 872 { 873 PetscErrorCode ierr; 874 875 PetscFunctionBegin; 876 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid); 877 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);CHKERRQ(ierr); 878 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));CHKERRQ(ierr); 879 880 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 881 882 mumps->id.job = JOB_INIT; 883 mumps->id.par = 1; /* host participates factorizaton and solve */ 884 mumps->id.sym = mumps->sym; 885 PetscMUMPS_c(&mumps->id); 886 887 mumps->CleanUpMUMPS = PETSC_FALSE; 888 mumps->scat_rhs = NULL; 889 mumps->scat_sol = NULL; 890 891 /* set PETSc-MUMPS default options - override MUMPS default */ 892 mumps->id.ICNTL(3) = 0; 893 mumps->id.ICNTL(4) = 0; 894 if (mumps->size == 1) { 895 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 896 } else { 897 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 898 mumps->id.ICNTL(21) = 1; /* distributed solution */ 899 } 900 901 /* schur */ 902 mumps->id.size_schur = 0; 903 mumps->id.listvar_schur = NULL; 904 mumps->id.schur = NULL; 905 PetscFunctionReturn(0); 906 } 907 908 /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */ 909 #undef __FUNCT__ 910 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 911 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 912 { 913 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 914 PetscErrorCode ierr; 915 Vec b; 916 IS is_iden; 917 const PetscInt M = A->rmap->N; 918 919 PetscFunctionBegin; 920 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 921 922 /* Set MUMPS options from the options database */ 923 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 924 925 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 926 927 /* analysis phase */ 928 /*----------------*/ 929 mumps->id.job = JOB_FACTSYMBOLIC; 930 mumps->id.n = M; 931 switch (mumps->id.ICNTL(18)) { 932 case 0: /* centralized assembled matrix input */ 933 if (!mumps->myid) { 934 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 935 if (mumps->id.ICNTL(6)>1) { 936 #if defined(PETSC_USE_COMPLEX) 937 #if defined(PETSC_USE_REAL_SINGLE) 938 mumps->id.a = (mumps_complex*)mumps->val; 939 #else 940 mumps->id.a = (mumps_double_complex*)mumps->val; 941 #endif 942 #else 943 mumps->id.a = mumps->val; 944 #endif 945 } 946 if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */ 947 /* 948 PetscBool flag; 949 ierr = ISEqual(r,c,&flag);CHKERRQ(ierr); 950 if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm"); 951 ierr = ISView(r,PETSC_VIEWER_STDOUT_SELF); 952 */ 953 if (!mumps->myid) { 954 const PetscInt *idx; 955 PetscInt i,*perm_in; 956 957 ierr = PetscMalloc1(M,&perm_in);CHKERRQ(ierr); 958 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 959 960 mumps->id.perm_in = perm_in; 961 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 962 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 963 } 964 } 965 } 966 break; 967 case 3: /* distributed assembled matrix input (size>1) */ 968 mumps->id.nz_loc = mumps->nz; 969 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 970 if (mumps->id.ICNTL(6)>1) { 971 #if defined(PETSC_USE_COMPLEX) 972 #if defined(PETSC_USE_REAL_SINGLE) 973 mumps->id.a_loc = (mumps_complex*)mumps->val; 974 #else 975 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 976 #endif 977 #else 978 mumps->id.a_loc = mumps->val; 979 #endif 980 } 981 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 982 if (!mumps->myid) { 983 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 984 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 985 } else { 986 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 987 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 988 } 989 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 990 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 991 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 992 ierr = VecDestroy(&b);CHKERRQ(ierr); 993 break; 994 } 995 PetscMUMPS_c(&mumps->id); 996 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)); 997 998 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 999 F->ops->solve = MatSolve_MUMPS; 1000 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1001 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 1002 PetscFunctionReturn(0); 1003 } 1004 1005 /* Note the Petsc r and c permutations are ignored */ 1006 #undef __FUNCT__ 1007 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 1008 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1009 { 1010 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1011 PetscErrorCode ierr; 1012 Vec b; 1013 IS is_iden; 1014 const PetscInt M = A->rmap->N; 1015 1016 PetscFunctionBegin; 1017 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1018 1019 /* Set MUMPS options from the options database */ 1020 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1021 1022 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1023 1024 /* analysis phase */ 1025 /*----------------*/ 1026 mumps->id.job = JOB_FACTSYMBOLIC; 1027 mumps->id.n = M; 1028 switch (mumps->id.ICNTL(18)) { 1029 case 0: /* centralized assembled matrix input */ 1030 if (!mumps->myid) { 1031 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1032 if (mumps->id.ICNTL(6)>1) { 1033 #if defined(PETSC_USE_COMPLEX) 1034 #if defined(PETSC_USE_REAL_SINGLE) 1035 mumps->id.a = (mumps_complex*)mumps->val; 1036 #else 1037 mumps->id.a = (mumps_double_complex*)mumps->val; 1038 #endif 1039 #else 1040 mumps->id.a = mumps->val; 1041 #endif 1042 } 1043 } 1044 break; 1045 case 3: /* distributed assembled matrix input (size>1) */ 1046 mumps->id.nz_loc = mumps->nz; 1047 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1048 if (mumps->id.ICNTL(6)>1) { 1049 #if defined(PETSC_USE_COMPLEX) 1050 #if defined(PETSC_USE_REAL_SINGLE) 1051 mumps->id.a_loc = (mumps_complex*)mumps->val; 1052 #else 1053 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1054 #endif 1055 #else 1056 mumps->id.a_loc = mumps->val; 1057 #endif 1058 } 1059 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1060 if (!mumps->myid) { 1061 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1062 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1063 } else { 1064 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1065 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1066 } 1067 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1068 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1069 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1070 ierr = VecDestroy(&b);CHKERRQ(ierr); 1071 break; 1072 } 1073 PetscMUMPS_c(&mumps->id); 1074 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)); 1075 1076 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1077 F->ops->solve = MatSolve_MUMPS; 1078 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1079 PetscFunctionReturn(0); 1080 } 1081 1082 /* Note the Petsc r permutation and factor info are ignored */ 1083 #undef __FUNCT__ 1084 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 1085 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 1086 { 1087 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1088 PetscErrorCode ierr; 1089 Vec b; 1090 IS is_iden; 1091 const PetscInt M = A->rmap->N; 1092 1093 PetscFunctionBegin; 1094 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1095 1096 /* Set MUMPS options from the options database */ 1097 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1098 1099 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1100 1101 /* analysis phase */ 1102 /*----------------*/ 1103 mumps->id.job = JOB_FACTSYMBOLIC; 1104 mumps->id.n = M; 1105 switch (mumps->id.ICNTL(18)) { 1106 case 0: /* centralized assembled matrix input */ 1107 if (!mumps->myid) { 1108 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1109 if (mumps->id.ICNTL(6)>1) { 1110 #if defined(PETSC_USE_COMPLEX) 1111 #if defined(PETSC_USE_REAL_SINGLE) 1112 mumps->id.a = (mumps_complex*)mumps->val; 1113 #else 1114 mumps->id.a = (mumps_double_complex*)mumps->val; 1115 #endif 1116 #else 1117 mumps->id.a = mumps->val; 1118 #endif 1119 } 1120 } 1121 break; 1122 case 3: /* distributed assembled matrix input (size>1) */ 1123 mumps->id.nz_loc = mumps->nz; 1124 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1125 if (mumps->id.ICNTL(6)>1) { 1126 #if defined(PETSC_USE_COMPLEX) 1127 #if defined(PETSC_USE_REAL_SINGLE) 1128 mumps->id.a_loc = (mumps_complex*)mumps->val; 1129 #else 1130 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1131 #endif 1132 #else 1133 mumps->id.a_loc = mumps->val; 1134 #endif 1135 } 1136 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1137 if (!mumps->myid) { 1138 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1139 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1140 } else { 1141 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1142 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1143 } 1144 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1145 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1146 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1147 ierr = VecDestroy(&b);CHKERRQ(ierr); 1148 break; 1149 } 1150 PetscMUMPS_c(&mumps->id); 1151 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)); 1152 1153 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1154 F->ops->solve = MatSolve_MUMPS; 1155 F->ops->solvetranspose = MatSolve_MUMPS; 1156 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 1157 #if !defined(PETSC_USE_COMPLEX) 1158 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1159 #else 1160 F->ops->getinertia = NULL; 1161 #endif 1162 PetscFunctionReturn(0); 1163 } 1164 1165 #undef __FUNCT__ 1166 #define __FUNCT__ "MatView_MUMPS" 1167 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1168 { 1169 PetscErrorCode ierr; 1170 PetscBool iascii; 1171 PetscViewerFormat format; 1172 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1173 1174 PetscFunctionBegin; 1175 /* check if matrix is mumps type */ 1176 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1177 1178 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1179 if (iascii) { 1180 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1181 if (format == PETSC_VIEWER_ASCII_INFO) { 1182 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1183 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1184 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1185 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));CHKERRQ(ierr); 1186 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));CHKERRQ(ierr); 1187 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));CHKERRQ(ierr); 1188 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));CHKERRQ(ierr); 1189 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));CHKERRQ(ierr); 1190 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));CHKERRQ(ierr); 1191 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));CHKERRQ(ierr); 1192 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));CHKERRQ(ierr); 1193 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));CHKERRQ(ierr); 1194 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));CHKERRQ(ierr); 1195 if (mumps->id.ICNTL(11)>0) { 1196 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));CHKERRQ(ierr); 1197 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));CHKERRQ(ierr); 1198 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));CHKERRQ(ierr); 1199 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));CHKERRQ(ierr); 1200 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));CHKERRQ(ierr); 1201 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));CHKERRQ(ierr); 1202 } 1203 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));CHKERRQ(ierr); 1204 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));CHKERRQ(ierr); 1205 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));CHKERRQ(ierr); 1206 /* ICNTL(15-17) not used */ 1207 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));CHKERRQ(ierr); 1208 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));CHKERRQ(ierr); 1209 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));CHKERRQ(ierr); 1210 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));CHKERRQ(ierr); 1211 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));CHKERRQ(ierr); 1212 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));CHKERRQ(ierr); 1213 1214 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));CHKERRQ(ierr); 1215 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));CHKERRQ(ierr); 1216 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));CHKERRQ(ierr); 1217 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));CHKERRQ(ierr); 1218 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));CHKERRQ(ierr); 1219 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));CHKERRQ(ierr); 1220 1221 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));CHKERRQ(ierr); 1222 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));CHKERRQ(ierr); 1223 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));CHKERRQ(ierr); 1224 1225 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));CHKERRQ(ierr); 1226 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));CHKERRQ(ierr); 1227 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));CHKERRQ(ierr); 1228 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));CHKERRQ(ierr); 1229 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));CHKERRQ(ierr); 1230 1231 /* infomation local to each processor */ 1232 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1233 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1234 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));CHKERRQ(ierr); 1235 ierr = PetscViewerFlush(viewer); 1236 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1237 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));CHKERRQ(ierr); 1238 ierr = PetscViewerFlush(viewer); 1239 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1240 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));CHKERRQ(ierr); 1241 ierr = PetscViewerFlush(viewer); 1242 1243 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1244 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));CHKERRQ(ierr); 1245 ierr = PetscViewerFlush(viewer); 1246 1247 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1248 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));CHKERRQ(ierr); 1249 ierr = PetscViewerFlush(viewer); 1250 1251 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1252 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));CHKERRQ(ierr); 1253 ierr = PetscViewerFlush(viewer); 1254 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1255 1256 if (!mumps->myid) { /* information from the host */ 1257 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));CHKERRQ(ierr); 1258 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));CHKERRQ(ierr); 1259 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));CHKERRQ(ierr); 1260 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); 1261 1262 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));CHKERRQ(ierr); 1263 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));CHKERRQ(ierr); 1264 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));CHKERRQ(ierr); 1265 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));CHKERRQ(ierr); 1266 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));CHKERRQ(ierr); 1267 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));CHKERRQ(ierr); 1268 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));CHKERRQ(ierr); 1269 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));CHKERRQ(ierr); 1270 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));CHKERRQ(ierr); 1271 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));CHKERRQ(ierr); 1272 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));CHKERRQ(ierr); 1273 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));CHKERRQ(ierr); 1274 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));CHKERRQ(ierr); 1275 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); 1276 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); 1277 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); 1278 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); 1279 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));CHKERRQ(ierr); 1280 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); 1281 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); 1282 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));CHKERRQ(ierr); 1283 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));CHKERRQ(ierr); 1284 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));CHKERRQ(ierr); 1285 ierr = PetscViewerASCIIPrintf(viewer," INFOG(28) (after factorization: number of null pivots encountered): %d\n",mumps->id.INFOG(28));CHKERRQ(ierr); 1286 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); 1287 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); 1288 ierr = PetscViewerASCIIPrintf(viewer," INFOG(32) (after analysis: type of analysis done): %d\n",mumps->id.INFOG(32));CHKERRQ(ierr); 1289 ierr = PetscViewerASCIIPrintf(viewer," INFOG(33) (value used for ICNTL(8)): %d\n",mumps->id.INFOG(33));CHKERRQ(ierr); 1290 ierr = PetscViewerASCIIPrintf(viewer," INFOG(34) (exponent of the determinant if determinant is requested): %d\n",mumps->id.INFOG(34));CHKERRQ(ierr); 1291 } 1292 } 1293 } 1294 PetscFunctionReturn(0); 1295 } 1296 1297 #undef __FUNCT__ 1298 #define __FUNCT__ "MatGetInfo_MUMPS" 1299 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1300 { 1301 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1302 1303 PetscFunctionBegin; 1304 info->block_size = 1.0; 1305 info->nz_allocated = mumps->id.INFOG(20); 1306 info->nz_used = mumps->id.INFOG(20); 1307 info->nz_unneeded = 0.0; 1308 info->assemblies = 0.0; 1309 info->mallocs = 0.0; 1310 info->memory = 0.0; 1311 info->fill_ratio_given = 0; 1312 info->fill_ratio_needed = 0; 1313 info->factor_mallocs = 0; 1314 PetscFunctionReturn(0); 1315 } 1316 1317 /* -------------------------------------------------------------------------------------------*/ 1318 #undef __FUNCT__ 1319 #define __FUNCT__ "MatMumpsSetSchurIndices_MUMPS" 1320 PetscErrorCode MatMumpsSetSchurIndices_MUMPS(Mat F,PetscInt size,PetscInt idxs[]) 1321 { 1322 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1323 PetscErrorCode ierr; 1324 1325 PetscFunctionBegin; 1326 if (mumps->size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MUMPS parallel computation of Schur complements not yet supported from PETSc\n"); 1327 if (mumps->id.size_schur != size) { 1328 ierr = PetscFree2(mumps->id.listvar_schur,mumps->id.schur);CHKERRQ(ierr); 1329 mumps->id.size_schur = size; 1330 mumps->id.schur_lld = size; 1331 ierr = PetscMalloc2(size,&mumps->id.listvar_schur,size*size,&mumps->id.schur);CHKERRQ(ierr); 1332 } 1333 ierr = PetscMemcpy(mumps->id.listvar_schur,idxs,size*sizeof(PetscInt));CHKERRQ(ierr); 1334 if (F->factortype == MAT_FACTOR_LU) { 1335 mumps->id.ICNTL(19) = 3; /* return full matrix */ 1336 } else { 1337 mumps->id.ICNTL(19) = 2; /* return lower triangular part */ 1338 } 1339 PetscFunctionReturn(0); 1340 } 1341 1342 #undef __FUNCT__ 1343 #define __FUNCT__ "MatMumpsSetSchurIndices" 1344 /*@ 1345 MatMumpsSetSchurIndices - Set indices defining the Schur complement that MUMPS will compute during the factorization steps 1346 1347 Logically Collective on Mat 1348 1349 Input Parameters: 1350 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1351 . size - size of the Schur complement indices 1352 - idxs[] - array of Schur complement indices 1353 1354 Notes: 1355 The user has to free the array idxs[] since it is copied by the routine. 1356 Currently implemented for sequential matrices 1357 1358 Level: advanced 1359 1360 References: MUMPS Users' Guide 1361 1362 .seealso: MatGetFactor() 1363 @*/ 1364 PetscErrorCode MatMumpsSetSchurIndices(Mat F,PetscInt size,PetscInt idxs[]) 1365 { 1366 PetscErrorCode ierr; 1367 1368 PetscFunctionBegin; 1369 PetscValidIntPointer(idxs,3); 1370 ierr = PetscTryMethod(F,"MatMumpsSetSchurIndices_C",(Mat,PetscInt,PetscInt[]),(F,size,idxs));CHKERRQ(ierr); 1371 PetscFunctionReturn(0); 1372 } 1373 1374 /* -------------------------------------------------------------------------------------------*/ 1375 #undef __FUNCT__ 1376 #define __FUNCT__ "MatMumpsGetSchurComplement_MUMPS" 1377 PetscErrorCode MatMumpsGetSchurComplement_MUMPS(Mat F,Mat* S) 1378 { 1379 Mat St; 1380 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1381 PetscScalar *array; 1382 PetscErrorCode ierr; 1383 1384 PetscFunctionBegin; 1385 if (mumps->id.job != JOB_FACTNUMERIC) { 1386 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_WRONG,"Numerical factorization phase not yet performed! You should call MatFactorSymbolic/Numeric before"); 1387 } else if (mumps->id.size_schur == 0) { 1388 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_WRONG,"Schur indices not set! You should call MatMumpsSetSchurIndices before"); 1389 } 1390 1391 ierr = MatCreate(PetscObjectComm((PetscObject)F),&St);CHKERRQ(ierr); 1392 ierr = MatSetSizes(St,PETSC_DECIDE,PETSC_DECIDE,mumps->id.size_schur,mumps->id.size_schur);CHKERRQ(ierr); 1393 ierr = MatSetType(St,MATDENSE);CHKERRQ(ierr); 1394 ierr = MatDenseGetArray(St,&array);CHKERRQ(ierr); 1395 if (mumps->sym == 0) { /* MUMPS returned a full matrix */ 1396 if (mumps->id.ICNTL(19) == 1) { 1397 PetscInt i,j,ncols=F->cmap->n,nrows=F->rmap->n; 1398 for (i=0;i<ncols;i++) { 1399 for (j=0;j<nrows;j++) { 1400 array[i*nrows+j] = mumps->id.schur[j*ncols+i]; 1401 } 1402 } 1403 } else { 1404 ierr = PetscMemcpy(array,mumps->id.schur,mumps->id.size_schur*mumps->id.size_schur*sizeof(PetscScalar));CHKERRQ(ierr); 1405 } 1406 } else { 1407 if (mumps->id.ICNTL(19) == 3) { 1408 ierr = PetscMemcpy(array,mumps->id.schur,mumps->id.size_schur*mumps->id.size_schur*sizeof(PetscScalar));CHKERRQ(ierr); 1409 } else { 1410 PetscInt i,j,ncols=F->cmap->n,nrows=F->rmap->n; 1411 for (i=0;i<ncols;i++) { 1412 array[i*nrows+i] = mumps->id.schur[i*nrows+i]; 1413 for (j=i+1;j<nrows;j++) { 1414 array[i*nrows+j] = mumps->id.schur[i*nrows+j]; 1415 array[j*nrows+i] = mumps->id.schur[i*nrows+j]; 1416 } 1417 } 1418 } 1419 } 1420 ierr = MatDenseRestoreArray(St,&array);CHKERRQ(ierr); 1421 *S = St; 1422 PetscFunctionReturn(0); 1423 } 1424 1425 #undef __FUNCT__ 1426 #define __FUNCT__ "MatMumpsGetSchurComplement" 1427 /*@ 1428 MatMumpsGetSchurComplement - Get Schur complement matrix computed by MUMPS during the factorization step 1429 1430 Logically Collective on Mat 1431 1432 Input Parameters: 1433 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1434 . *S - location where to return the Schur complement (MATDENSE) 1435 1436 Notes: 1437 Currently implemented for sequential matrices 1438 1439 Level: advanced 1440 1441 References: MUMPS Users' Guide 1442 1443 .seealso: MatGetFactor() 1444 @*/ 1445 PetscErrorCode MatMumpsGetSchurComplement(Mat F,Mat* S) 1446 { 1447 PetscErrorCode ierr; 1448 1449 PetscFunctionBegin; 1450 ierr = PetscTryMethod(F,"MatMumpsGetSchurComplement_C",(Mat,Mat*),(F,S));CHKERRQ(ierr); 1451 PetscFunctionReturn(0); 1452 } 1453 1454 /* -------------------------------------------------------------------------------------------*/ 1455 #undef __FUNCT__ 1456 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1457 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1458 { 1459 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1460 1461 PetscFunctionBegin; 1462 mumps->id.ICNTL(icntl) = ival; 1463 PetscFunctionReturn(0); 1464 } 1465 1466 #undef __FUNCT__ 1467 #define __FUNCT__ "MatMumpsGetIcntl_MUMPS" 1468 PetscErrorCode MatMumpsGetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt *ival) 1469 { 1470 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1471 1472 PetscFunctionBegin; 1473 *ival = mumps->id.ICNTL(icntl); 1474 PetscFunctionReturn(0); 1475 } 1476 1477 #undef __FUNCT__ 1478 #define __FUNCT__ "MatMumpsSetIcntl" 1479 /*@ 1480 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1481 1482 Logically Collective on Mat 1483 1484 Input Parameters: 1485 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1486 . icntl - index of MUMPS parameter array ICNTL() 1487 - ival - value of MUMPS ICNTL(icntl) 1488 1489 Options Database: 1490 . -mat_mumps_icntl_<icntl> <ival> 1491 1492 Level: beginner 1493 1494 References: MUMPS Users' Guide 1495 1496 .seealso: MatGetFactor() 1497 @*/ 1498 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1499 { 1500 PetscErrorCode ierr; 1501 1502 PetscFunctionBegin; 1503 PetscValidLogicalCollectiveInt(F,icntl,2); 1504 PetscValidLogicalCollectiveInt(F,ival,3); 1505 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1506 PetscFunctionReturn(0); 1507 } 1508 1509 #undef __FUNCT__ 1510 #define __FUNCT__ "MatMumpsGetIcntl" 1511 /*@ 1512 MatMumpsGetIcntl - Get MUMPS parameter ICNTL() 1513 1514 Logically Collective on Mat 1515 1516 Input Parameters: 1517 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1518 - icntl - index of MUMPS parameter array ICNTL() 1519 1520 Output Parameter: 1521 . ival - value of MUMPS ICNTL(icntl) 1522 1523 Level: beginner 1524 1525 References: MUMPS Users' Guide 1526 1527 .seealso: MatGetFactor() 1528 @*/ 1529 PetscErrorCode MatMumpsGetIcntl(Mat F,PetscInt icntl,PetscInt *ival) 1530 { 1531 PetscErrorCode ierr; 1532 1533 PetscFunctionBegin; 1534 PetscValidLogicalCollectiveInt(F,icntl,2); 1535 PetscValidIntPointer(ival,3); 1536 ierr = PetscTryMethod(F,"MatMumpsGetIcntl_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1537 PetscFunctionReturn(0); 1538 } 1539 1540 /* -------------------------------------------------------------------------------------------*/ 1541 #undef __FUNCT__ 1542 #define __FUNCT__ "MatMumpsSetCntl_MUMPS" 1543 PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val) 1544 { 1545 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1546 1547 PetscFunctionBegin; 1548 mumps->id.CNTL(icntl) = val; 1549 PetscFunctionReturn(0); 1550 } 1551 1552 #undef __FUNCT__ 1553 #define __FUNCT__ "MatMumpsGetCntl_MUMPS" 1554 PetscErrorCode MatMumpsGetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal *val) 1555 { 1556 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1557 1558 PetscFunctionBegin; 1559 *val = mumps->id.CNTL(icntl); 1560 PetscFunctionReturn(0); 1561 } 1562 1563 #undef __FUNCT__ 1564 #define __FUNCT__ "MatMumpsSetCntl" 1565 /*@ 1566 MatMumpsSetCntl - Set MUMPS parameter CNTL() 1567 1568 Logically Collective on Mat 1569 1570 Input Parameters: 1571 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1572 . icntl - index of MUMPS parameter array CNTL() 1573 - val - value of MUMPS CNTL(icntl) 1574 1575 Options Database: 1576 . -mat_mumps_cntl_<icntl> <val> 1577 1578 Level: beginner 1579 1580 References: MUMPS Users' Guide 1581 1582 .seealso: MatGetFactor() 1583 @*/ 1584 PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val) 1585 { 1586 PetscErrorCode ierr; 1587 1588 PetscFunctionBegin; 1589 PetscValidLogicalCollectiveInt(F,icntl,2); 1590 PetscValidLogicalCollectiveReal(F,val,3); 1591 ierr = PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));CHKERRQ(ierr); 1592 PetscFunctionReturn(0); 1593 } 1594 1595 #undef __FUNCT__ 1596 #define __FUNCT__ "MatMumpsGetCntl" 1597 /*@ 1598 MatMumpsGetCntl - Get MUMPS parameter CNTL() 1599 1600 Logically Collective on Mat 1601 1602 Input Parameters: 1603 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1604 - icntl - index of MUMPS parameter array CNTL() 1605 1606 Output Parameter: 1607 . val - value of MUMPS CNTL(icntl) 1608 1609 Level: beginner 1610 1611 References: MUMPS Users' Guide 1612 1613 .seealso: MatGetFactor() 1614 @*/ 1615 PetscErrorCode MatMumpsGetCntl(Mat F,PetscInt icntl,PetscReal *val) 1616 { 1617 PetscErrorCode ierr; 1618 1619 PetscFunctionBegin; 1620 PetscValidLogicalCollectiveInt(F,icntl,2); 1621 PetscValidRealPointer(val,3); 1622 ierr = PetscTryMethod(F,"MatMumpsGetCntl_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1623 PetscFunctionReturn(0); 1624 } 1625 1626 #undef __FUNCT__ 1627 #define __FUNCT__ "MatMumpsGetInfo_MUMPS" 1628 PetscErrorCode MatMumpsGetInfo_MUMPS(Mat F,PetscInt icntl,PetscInt *info) 1629 { 1630 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1631 1632 PetscFunctionBegin; 1633 *info = mumps->id.INFO(icntl); 1634 PetscFunctionReturn(0); 1635 } 1636 1637 #undef __FUNCT__ 1638 #define __FUNCT__ "MatMumpsGetInfog_MUMPS" 1639 PetscErrorCode MatMumpsGetInfog_MUMPS(Mat F,PetscInt icntl,PetscInt *infog) 1640 { 1641 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1642 1643 PetscFunctionBegin; 1644 *infog = mumps->id.INFOG(icntl); 1645 PetscFunctionReturn(0); 1646 } 1647 1648 #undef __FUNCT__ 1649 #define __FUNCT__ "MatMumpsGetRinfo_MUMPS" 1650 PetscErrorCode MatMumpsGetRinfo_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfo) 1651 { 1652 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1653 1654 PetscFunctionBegin; 1655 *rinfo = mumps->id.RINFO(icntl); 1656 PetscFunctionReturn(0); 1657 } 1658 1659 #undef __FUNCT__ 1660 #define __FUNCT__ "MatMumpsGetRinfog_MUMPS" 1661 PetscErrorCode MatMumpsGetRinfog_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfog) 1662 { 1663 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1664 1665 PetscFunctionBegin; 1666 *rinfog = mumps->id.RINFOG(icntl); 1667 PetscFunctionReturn(0); 1668 } 1669 1670 #undef __FUNCT__ 1671 #define __FUNCT__ "MatMumpsGetInfo" 1672 /*@ 1673 MatMumpsGetInfo - Get MUMPS parameter INFO() 1674 1675 Logically Collective on Mat 1676 1677 Input Parameters: 1678 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1679 - icntl - index of MUMPS parameter array INFO() 1680 1681 Output Parameter: 1682 . ival - value of MUMPS INFO(icntl) 1683 1684 Level: beginner 1685 1686 References: MUMPS Users' Guide 1687 1688 .seealso: MatGetFactor() 1689 @*/ 1690 PetscErrorCode MatMumpsGetInfo(Mat F,PetscInt icntl,PetscInt *ival) 1691 { 1692 PetscErrorCode ierr; 1693 1694 PetscFunctionBegin; 1695 PetscValidIntPointer(ival,3); 1696 ierr = PetscTryMethod(F,"MatMumpsGetInfo_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1697 PetscFunctionReturn(0); 1698 } 1699 1700 #undef __FUNCT__ 1701 #define __FUNCT__ "MatMumpsGetInfog" 1702 /*@ 1703 MatMumpsGetInfog - Get MUMPS parameter INFOG() 1704 1705 Logically Collective on Mat 1706 1707 Input Parameters: 1708 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1709 - icntl - index of MUMPS parameter array INFOG() 1710 1711 Output Parameter: 1712 . ival - value of MUMPS INFOG(icntl) 1713 1714 Level: beginner 1715 1716 References: MUMPS Users' Guide 1717 1718 .seealso: MatGetFactor() 1719 @*/ 1720 PetscErrorCode MatMumpsGetInfog(Mat F,PetscInt icntl,PetscInt *ival) 1721 { 1722 PetscErrorCode ierr; 1723 1724 PetscFunctionBegin; 1725 PetscValidIntPointer(ival,3); 1726 ierr = PetscTryMethod(F,"MatMumpsGetInfog_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1727 PetscFunctionReturn(0); 1728 } 1729 1730 #undef __FUNCT__ 1731 #define __FUNCT__ "MatMumpsGetRinfo" 1732 /*@ 1733 MatMumpsGetRinfo - Get MUMPS parameter RINFO() 1734 1735 Logically Collective on Mat 1736 1737 Input Parameters: 1738 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1739 - icntl - index of MUMPS parameter array RINFO() 1740 1741 Output Parameter: 1742 . val - value of MUMPS RINFO(icntl) 1743 1744 Level: beginner 1745 1746 References: MUMPS Users' Guide 1747 1748 .seealso: MatGetFactor() 1749 @*/ 1750 PetscErrorCode MatMumpsGetRinfo(Mat F,PetscInt icntl,PetscReal *val) 1751 { 1752 PetscErrorCode ierr; 1753 1754 PetscFunctionBegin; 1755 PetscValidRealPointer(val,3); 1756 ierr = PetscTryMethod(F,"MatMumpsGetRinfo_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1757 PetscFunctionReturn(0); 1758 } 1759 1760 #undef __FUNCT__ 1761 #define __FUNCT__ "MatMumpsGetRinfog" 1762 /*@ 1763 MatMumpsGetRinfog - Get MUMPS parameter RINFOG() 1764 1765 Logically Collective on Mat 1766 1767 Input Parameters: 1768 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1769 - icntl - index of MUMPS parameter array RINFOG() 1770 1771 Output Parameter: 1772 . val - value of MUMPS RINFOG(icntl) 1773 1774 Level: beginner 1775 1776 References: MUMPS Users' Guide 1777 1778 .seealso: MatGetFactor() 1779 @*/ 1780 PetscErrorCode MatMumpsGetRinfog(Mat F,PetscInt icntl,PetscReal *val) 1781 { 1782 PetscErrorCode ierr; 1783 1784 PetscFunctionBegin; 1785 PetscValidRealPointer(val,3); 1786 ierr = PetscTryMethod(F,"MatMumpsGetRinfog_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1787 PetscFunctionReturn(0); 1788 } 1789 1790 /*MC 1791 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1792 distributed and sequential matrices via the external package MUMPS. 1793 1794 Works with MATAIJ and MATSBAIJ matrices 1795 1796 Options Database Keys: 1797 + -mat_mumps_icntl_4 <0,...,4> - print level 1798 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1799 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec) 1800 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1801 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1802 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1803 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1804 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1805 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1806 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1807 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1808 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1809 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1810 1811 Level: beginner 1812 1813 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1814 1815 M*/ 1816 1817 #undef __FUNCT__ 1818 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1819 static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1820 { 1821 PetscFunctionBegin; 1822 *type = MATSOLVERMUMPS; 1823 PetscFunctionReturn(0); 1824 } 1825 1826 /* MatGetFactor for Seq and MPI AIJ matrices */ 1827 #undef __FUNCT__ 1828 #define __FUNCT__ "MatGetFactor_aij_mumps" 1829 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1830 { 1831 Mat B; 1832 PetscErrorCode ierr; 1833 Mat_MUMPS *mumps; 1834 PetscBool isSeqAIJ; 1835 1836 PetscFunctionBegin; 1837 /* Create the factorization matrix */ 1838 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1839 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1840 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1841 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1842 if (isSeqAIJ) { 1843 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1844 } else { 1845 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1846 } 1847 1848 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1849 1850 B->ops->view = MatView_MUMPS; 1851 B->ops->getinfo = MatGetInfo_MUMPS; 1852 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1853 1854 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1855 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1856 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1857 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1858 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1859 1860 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1861 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1862 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1863 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1864 1865 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 1866 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 1867 if (ftype == MAT_FACTOR_LU) { 1868 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1869 B->factortype = MAT_FACTOR_LU; 1870 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1871 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1872 mumps->sym = 0; 1873 } else { 1874 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1875 B->factortype = MAT_FACTOR_CHOLESKY; 1876 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1877 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1878 if (A->spd_set && A->spd) mumps->sym = 1; 1879 else mumps->sym = 2; 1880 } 1881 1882 mumps->isAIJ = PETSC_TRUE; 1883 mumps->Destroy = B->ops->destroy; 1884 B->ops->destroy = MatDestroy_MUMPS; 1885 B->spptr = (void*)mumps; 1886 1887 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1888 1889 *F = B; 1890 PetscFunctionReturn(0); 1891 } 1892 1893 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1894 #undef __FUNCT__ 1895 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1896 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1897 { 1898 Mat B; 1899 PetscErrorCode ierr; 1900 Mat_MUMPS *mumps; 1901 PetscBool isSeqSBAIJ; 1902 1903 PetscFunctionBegin; 1904 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1905 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"); 1906 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1907 /* Create the factorization matrix */ 1908 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1909 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1910 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1911 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1912 if (isSeqSBAIJ) { 1913 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 1914 1915 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1916 } else { 1917 ierr = MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);CHKERRQ(ierr); 1918 1919 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1920 } 1921 1922 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1923 B->ops->view = MatView_MUMPS; 1924 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1925 1926 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1927 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1928 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1929 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1930 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1931 1932 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1933 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1934 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1935 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1936 1937 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 1938 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 1939 1940 B->factortype = MAT_FACTOR_CHOLESKY; 1941 if (A->spd_set && A->spd) mumps->sym = 1; 1942 else mumps->sym = 2; 1943 1944 mumps->isAIJ = PETSC_FALSE; 1945 mumps->Destroy = B->ops->destroy; 1946 B->ops->destroy = MatDestroy_MUMPS; 1947 B->spptr = (void*)mumps; 1948 1949 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1950 1951 *F = B; 1952 PetscFunctionReturn(0); 1953 } 1954 1955 #undef __FUNCT__ 1956 #define __FUNCT__ "MatGetFactor_baij_mumps" 1957 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 1958 { 1959 Mat B; 1960 PetscErrorCode ierr; 1961 Mat_MUMPS *mumps; 1962 PetscBool isSeqBAIJ; 1963 1964 PetscFunctionBegin; 1965 /* Create the factorization matrix */ 1966 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1967 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1968 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1969 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1970 if (isSeqBAIJ) { 1971 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);CHKERRQ(ierr); 1972 } else { 1973 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 1974 } 1975 1976 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1977 if (ftype == MAT_FACTOR_LU) { 1978 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1979 B->factortype = MAT_FACTOR_LU; 1980 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 1981 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 1982 mumps->sym = 0; 1983 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 1984 1985 B->ops->view = MatView_MUMPS; 1986 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1987 1988 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1989 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1990 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1991 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1992 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1993 1994 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1995 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1996 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1997 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1998 1999 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 2000 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 2001 2002 mumps->isAIJ = PETSC_TRUE; 2003 mumps->Destroy = B->ops->destroy; 2004 B->ops->destroy = MatDestroy_MUMPS; 2005 B->spptr = (void*)mumps; 2006 2007 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 2008 2009 *F = B; 2010 PetscFunctionReturn(0); 2011 } 2012 2013 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 2014 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat,MatFactorType,Mat*); 2015 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 2016 2017 #undef __FUNCT__ 2018 #define __FUNCT__ "MatSolverPackageRegister_MUMPS" 2019 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MUMPS(void) 2020 { 2021 PetscErrorCode ierr; 2022 2023 PetscFunctionBegin; 2024 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2025 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2026 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2027 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2028 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPISBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2029 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2030 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2031 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2032 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2033 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2034 PetscFunctionReturn(0); 2035 } 2036 2037