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