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 99 #undef __FUNCT__ 100 #define __FUNCT__ "MatConvertToTriples_seqaij_seqaij" 101 PetscErrorCode MatConvertToTriples_seqaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 102 { 103 const PetscInt *ai,*aj,*ajj,M=A->rmap->n; 104 PetscInt nz,rnz,i,j; 105 PetscErrorCode ierr; 106 PetscInt *row,*col; 107 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 108 109 PetscFunctionBegin; 110 *v=aa->a; 111 if (reuse == MAT_INITIAL_MATRIX) { 112 nz = aa->nz; 113 ai = aa->i; 114 aj = aa->j; 115 *nnz = nz; 116 ierr = PetscMalloc(2*nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 117 col = row + nz; 118 119 nz = 0; 120 for (i=0; i<M; i++) { 121 rnz = ai[i+1] - ai[i]; 122 ajj = aj + ai[i]; 123 for (j=0; j<rnz; j++) { 124 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 125 } 126 } 127 *r = row; *c = col; 128 } 129 PetscFunctionReturn(0); 130 } 131 132 #undef __FUNCT__ 133 #define __FUNCT__ "MatConvertToTriples_seqbaij_seqaij" 134 PetscErrorCode MatConvertToTriples_seqbaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 135 { 136 Mat_SeqBAIJ *aa=(Mat_SeqBAIJ*)A->data; 137 const PetscInt *ai,*aj,*ajj,bs=A->rmap->bs,bs2=aa->bs2,M=A->rmap->N/bs; 138 PetscInt nz,idx=0,rnz,i,j,k,m; 139 PetscErrorCode ierr; 140 PetscInt *row,*col; 141 142 PetscFunctionBegin; 143 *v = aa->a; 144 if (reuse == MAT_INITIAL_MATRIX) { 145 ai = aa->i; aj = aa->j; 146 nz = bs2*aa->nz; 147 *nnz = nz; 148 ierr = PetscMalloc(2*nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 149 col = row + nz; 150 151 for (i=0; i<M; i++) { 152 ajj = aj + ai[i]; 153 rnz = ai[i+1] - ai[i]; 154 for (k=0; k<rnz; k++) { 155 for (j=0; j<bs; j++) { 156 for (m=0; m<bs; m++) { 157 row[idx] = i*bs + m + shift; 158 col[idx++] = bs*(ajj[k]) + j + shift; 159 } 160 } 161 } 162 } 163 *r = row; *c = col; 164 } 165 PetscFunctionReturn(0); 166 } 167 168 #undef __FUNCT__ 169 #define __FUNCT__ "MatConvertToTriples_seqsbaij_seqsbaij" 170 PetscErrorCode MatConvertToTriples_seqsbaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 171 { 172 const PetscInt *ai, *aj,*ajj,M=A->rmap->n; 173 PetscInt nz,rnz,i,j; 174 PetscErrorCode ierr; 175 PetscInt *row,*col; 176 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data; 177 178 PetscFunctionBegin; 179 *v = aa->a; 180 if (reuse == MAT_INITIAL_MATRIX) { 181 nz = aa->nz; 182 ai = aa->i; 183 aj = aa->j; 184 *v = aa->a; 185 *nnz = nz; 186 ierr = PetscMalloc(2*nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 187 col = row + nz; 188 189 nz = 0; 190 for (i=0; i<M; i++) { 191 rnz = ai[i+1] - ai[i]; 192 ajj = aj + ai[i]; 193 for (j=0; j<rnz; j++) { 194 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 195 } 196 } 197 *r = row; *c = col; 198 } 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "MatConvertToTriples_seqaij_seqsbaij" 204 PetscErrorCode MatConvertToTriples_seqaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 205 { 206 const PetscInt *ai,*aj,*ajj,*adiag,M=A->rmap->n; 207 PetscInt nz,rnz,i,j; 208 const PetscScalar *av,*v1; 209 PetscScalar *val; 210 PetscErrorCode ierr; 211 PetscInt *row,*col; 212 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data; 213 214 PetscFunctionBegin; 215 ai =aa->i; aj=aa->j;av=aa->a; 216 adiag=aa->diag; 217 if (reuse == MAT_INITIAL_MATRIX) { 218 nz = M + (aa->nz-M)/2; 219 *nnz = nz; 220 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 221 col = row + nz; 222 val = (PetscScalar*)(col + nz); 223 224 nz = 0; 225 for (i=0; i<M; i++) { 226 rnz = ai[i+1] - adiag[i]; 227 ajj = aj + adiag[i]; 228 v1 = av + adiag[i]; 229 for (j=0; j<rnz; j++) { 230 row[nz] = i+shift; col[nz] = ajj[j] + shift; val[nz++] = v1[j]; 231 } 232 } 233 *r = row; *c = col; *v = val; 234 } else { 235 nz = 0; val = *v; 236 for (i=0; i <M; i++) { 237 rnz = ai[i+1] - adiag[i]; 238 ajj = aj + adiag[i]; 239 v1 = av + adiag[i]; 240 for (j=0; j<rnz; j++) { 241 val[nz++] = v1[j]; 242 } 243 } 244 } 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "MatConvertToTriples_mpisbaij_mpisbaij" 250 PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 251 { 252 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 253 PetscErrorCode ierr; 254 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 255 PetscInt *row,*col; 256 const PetscScalar *av, *bv,*v1,*v2; 257 PetscScalar *val; 258 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 259 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)(mat->A)->data; 260 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 261 262 PetscFunctionBegin; 263 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 264 av=aa->a; bv=bb->a; 265 266 garray = mat->garray; 267 268 if (reuse == MAT_INITIAL_MATRIX) { 269 nz = aa->nz + bb->nz; 270 *nnz = nz; 271 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 272 col = row + nz; 273 val = (PetscScalar*)(col + nz); 274 275 *r = row; *c = col; *v = val; 276 } else { 277 row = *r; col = *c; val = *v; 278 } 279 280 jj = 0; irow = rstart; 281 for (i=0; i<m; i++) { 282 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 283 countA = ai[i+1] - ai[i]; 284 countB = bi[i+1] - bi[i]; 285 bjj = bj + bi[i]; 286 v1 = av + ai[i]; 287 v2 = bv + bi[i]; 288 289 /* A-part */ 290 for (j=0; j<countA; j++) { 291 if (reuse == MAT_INITIAL_MATRIX) { 292 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 293 } 294 val[jj++] = v1[j]; 295 } 296 297 /* B-part */ 298 for (j=0; j < countB; j++) { 299 if (reuse == MAT_INITIAL_MATRIX) { 300 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 301 } 302 val[jj++] = v2[j]; 303 } 304 irow++; 305 } 306 PetscFunctionReturn(0); 307 } 308 309 #undef __FUNCT__ 310 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpiaij" 311 PetscErrorCode MatConvertToTriples_mpiaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 312 { 313 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 314 PetscErrorCode ierr; 315 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 316 PetscInt *row,*col; 317 const PetscScalar *av, *bv,*v1,*v2; 318 PetscScalar *val; 319 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 320 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(mat->A)->data; 321 Mat_SeqAIJ *bb = (Mat_SeqAIJ*)(mat->B)->data; 322 323 PetscFunctionBegin; 324 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 325 av=aa->a; bv=bb->a; 326 327 garray = mat->garray; 328 329 if (reuse == MAT_INITIAL_MATRIX) { 330 nz = aa->nz + bb->nz; 331 *nnz = nz; 332 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 333 col = row + nz; 334 val = (PetscScalar*)(col + nz); 335 336 *r = row; *c = col; *v = val; 337 } else { 338 row = *r; col = *c; val = *v; 339 } 340 341 jj = 0; irow = rstart; 342 for (i=0; i<m; i++) { 343 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 344 countA = ai[i+1] - ai[i]; 345 countB = bi[i+1] - bi[i]; 346 bjj = bj + bi[i]; 347 v1 = av + ai[i]; 348 v2 = bv + bi[i]; 349 350 /* A-part */ 351 for (j=0; j<countA; j++) { 352 if (reuse == MAT_INITIAL_MATRIX) { 353 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 354 } 355 val[jj++] = v1[j]; 356 } 357 358 /* B-part */ 359 for (j=0; j < countB; j++) { 360 if (reuse == MAT_INITIAL_MATRIX) { 361 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 362 } 363 val[jj++] = v2[j]; 364 } 365 irow++; 366 } 367 PetscFunctionReturn(0); 368 } 369 370 #undef __FUNCT__ 371 #define __FUNCT__ "MatConvertToTriples_mpibaij_mpiaij" 372 PetscErrorCode MatConvertToTriples_mpibaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 373 { 374 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data; 375 Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)(mat->A)->data; 376 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 377 const PetscInt *ai = aa->i, *bi = bb->i, *aj = aa->j, *bj = bb->j,*ajj, *bjj; 378 const PetscInt *garray = mat->garray,mbs=mat->mbs,rstart=A->rmap->rstart; 379 const PetscInt bs = A->rmap->bs,bs2=mat->bs2; 380 PetscErrorCode ierr; 381 PetscInt nz,i,j,k,n,jj,irow,countA,countB,idx; 382 PetscInt *row,*col; 383 const PetscScalar *av=aa->a, *bv=bb->a,*v1,*v2; 384 PetscScalar *val; 385 386 PetscFunctionBegin; 387 if (reuse == MAT_INITIAL_MATRIX) { 388 nz = bs2*(aa->nz + bb->nz); 389 *nnz = nz; 390 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 391 col = row + nz; 392 val = (PetscScalar*)(col + nz); 393 394 *r = row; *c = col; *v = val; 395 } else { 396 row = *r; col = *c; val = *v; 397 } 398 399 jj = 0; irow = rstart; 400 for (i=0; i<mbs; i++) { 401 countA = ai[i+1] - ai[i]; 402 countB = bi[i+1] - bi[i]; 403 ajj = aj + ai[i]; 404 bjj = bj + bi[i]; 405 v1 = av + bs2*ai[i]; 406 v2 = bv + bs2*bi[i]; 407 408 idx = 0; 409 /* A-part */ 410 for (k=0; k<countA; k++) { 411 for (j=0; j<bs; j++) { 412 for (n=0; n<bs; n++) { 413 if (reuse == MAT_INITIAL_MATRIX) { 414 row[jj] = irow + n + shift; 415 col[jj] = rstart + bs*ajj[k] + j + shift; 416 } 417 val[jj++] = v1[idx++]; 418 } 419 } 420 } 421 422 idx = 0; 423 /* B-part */ 424 for (k=0; k<countB; k++) { 425 for (j=0; j<bs; j++) { 426 for (n=0; n<bs; n++) { 427 if (reuse == MAT_INITIAL_MATRIX) { 428 row[jj] = irow + n + shift; 429 col[jj] = bs*garray[bjj[k]] + j + shift; 430 } 431 val[jj++] = v2[idx++]; 432 } 433 } 434 } 435 irow += bs; 436 } 437 PetscFunctionReturn(0); 438 } 439 440 #undef __FUNCT__ 441 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpisbaij" 442 PetscErrorCode MatConvertToTriples_mpiaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 443 { 444 const PetscInt *ai, *aj,*adiag, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 445 PetscErrorCode ierr; 446 PetscInt rstart,nz,nza,nzb,i,j,jj,irow,countA,countB; 447 PetscInt *row,*col; 448 const PetscScalar *av, *bv,*v1,*v2; 449 PetscScalar *val; 450 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 451 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(mat->A)->data; 452 Mat_SeqAIJ *bb =(Mat_SeqAIJ*)(mat->B)->data; 453 454 PetscFunctionBegin; 455 ai=aa->i; aj=aa->j; adiag=aa->diag; 456 bi=bb->i; bj=bb->j; garray = mat->garray; 457 av=aa->a; bv=bb->a; 458 459 rstart = A->rmap->rstart; 460 461 if (reuse == MAT_INITIAL_MATRIX) { 462 nza = 0; /* num of upper triangular entries in mat->A, including diagonals */ 463 nzb = 0; /* num of upper triangular entries in mat->B */ 464 for (i=0; i<m; i++) { 465 nza += (ai[i+1] - adiag[i]); 466 countB = bi[i+1] - bi[i]; 467 bjj = bj + bi[i]; 468 for (j=0; j<countB; j++) { 469 if (garray[bjj[j]] > rstart) nzb++; 470 } 471 } 472 473 nz = nza + nzb; /* total nz of upper triangular part of mat */ 474 *nnz = nz; 475 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 476 col = row + nz; 477 val = (PetscScalar*)(col + nz); 478 479 *r = row; *c = col; *v = val; 480 } else { 481 row = *r; col = *c; val = *v; 482 } 483 484 jj = 0; irow = rstart; 485 for (i=0; i<m; i++) { 486 ajj = aj + adiag[i]; /* ptr to the beginning of the diagonal of this row */ 487 v1 = av + adiag[i]; 488 countA = ai[i+1] - adiag[i]; 489 countB = bi[i+1] - bi[i]; 490 bjj = bj + bi[i]; 491 v2 = bv + bi[i]; 492 493 /* A-part */ 494 for (j=0; j<countA; j++) { 495 if (reuse == MAT_INITIAL_MATRIX) { 496 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 497 } 498 val[jj++] = v1[j]; 499 } 500 501 /* B-part */ 502 for (j=0; j < countB; j++) { 503 if (garray[bjj[j]] > rstart) { 504 if (reuse == MAT_INITIAL_MATRIX) { 505 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 506 } 507 val[jj++] = v2[j]; 508 } 509 } 510 irow++; 511 } 512 PetscFunctionReturn(0); 513 } 514 515 #undef __FUNCT__ 516 #define __FUNCT__ "MatDestroy_MUMPS" 517 PetscErrorCode MatDestroy_MUMPS(Mat A) 518 { 519 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 520 PetscErrorCode ierr; 521 522 PetscFunctionBegin; 523 if (mumps->CleanUpMUMPS) { 524 /* Terminate instance, deallocate memories */ 525 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 526 ierr = VecScatterDestroy(&mumps->scat_rhs);CHKERRQ(ierr); 527 ierr = VecDestroy(&mumps->b_seq);CHKERRQ(ierr); 528 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 529 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 530 ierr = PetscFree(mumps->id.perm_in);CHKERRQ(ierr); 531 ierr = PetscFree(mumps->irn);CHKERRQ(ierr); 532 533 mumps->id.job = JOB_END; 534 PetscMUMPS_c(&mumps->id); 535 ierr = MPI_Comm_free(&(mumps->comm_mumps));CHKERRQ(ierr); 536 } 537 if (mumps->Destroy) { 538 ierr = (mumps->Destroy)(A);CHKERRQ(ierr); 539 } 540 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 541 542 /* clear composed functions */ 543 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);CHKERRQ(ierr); 544 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetIcntl_C",NULL);CHKERRQ(ierr); 545 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetCntl_C",NULL);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 #undef __FUNCT__ 550 #define __FUNCT__ "MatSolve_MUMPS" 551 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 552 { 553 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 554 PetscScalar *array; 555 Vec b_seq; 556 IS is_iden,is_petsc; 557 PetscErrorCode ierr; 558 PetscInt i; 559 static PetscBool cite1 = PETSC_FALSE,cite2 = PETSC_FALSE; 560 561 PetscFunctionBegin; 562 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); 563 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); 564 mumps->id.nrhs = 1; 565 b_seq = mumps->b_seq; 566 if (mumps->size > 1) { 567 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 568 ierr = VecScatterBegin(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 569 ierr = VecScatterEnd(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 570 if (!mumps->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 571 } else { /* size == 1 */ 572 ierr = VecCopy(b,x);CHKERRQ(ierr); 573 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 574 } 575 if (!mumps->myid) { /* define rhs on the host */ 576 mumps->id.nrhs = 1; 577 #if defined(PETSC_USE_COMPLEX) 578 #if defined(PETSC_USE_REAL_SINGLE) 579 mumps->id.rhs = (mumps_complex*)array; 580 #else 581 mumps->id.rhs = (mumps_double_complex*)array; 582 #endif 583 #else 584 mumps->id.rhs = array; 585 #endif 586 } 587 588 /* solve phase */ 589 /*-------------*/ 590 mumps->id.job = JOB_SOLVE; 591 PetscMUMPS_c(&mumps->id); 592 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)); 593 594 if (mumps->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 595 if (mumps->scat_sol && mumps->ICNTL9_pre != mumps->id.ICNTL(9)) { 596 /* when id.ICNTL(9) changes, the contents of lsol_loc may change (not its size, lsol_loc), recreates scat_sol */ 597 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 598 } 599 if (!mumps->scat_sol) { /* create scatter scat_sol */ 600 ierr = ISCreateStride(PETSC_COMM_SELF,mumps->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 601 for (i=0; i<mumps->id.lsol_loc; i++) { 602 mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 603 } 604 ierr = ISCreateGeneral(PETSC_COMM_SELF,mumps->id.lsol_loc,mumps->id.isol_loc,PETSC_COPY_VALUES,&is_petsc);CHKERRQ(ierr); /* to */ 605 ierr = VecScatterCreate(mumps->x_seq,is_iden,x,is_petsc,&mumps->scat_sol);CHKERRQ(ierr); 606 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 607 ierr = ISDestroy(&is_petsc);CHKERRQ(ierr); 608 609 mumps->ICNTL9_pre = mumps->id.ICNTL(9); /* save current value of id.ICNTL(9) */ 610 } 611 612 ierr = VecScatterBegin(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 613 ierr = VecScatterEnd(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 614 } 615 PetscFunctionReturn(0); 616 } 617 618 #undef __FUNCT__ 619 #define __FUNCT__ "MatSolveTranspose_MUMPS" 620 PetscErrorCode MatSolveTranspose_MUMPS(Mat A,Vec b,Vec x) 621 { 622 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 623 PetscErrorCode ierr; 624 625 PetscFunctionBegin; 626 mumps->id.ICNTL(9) = 0; 627 628 ierr = MatSolve_MUMPS(A,b,x);CHKERRQ(ierr); 629 630 mumps->id.ICNTL(9) = 1; 631 PetscFunctionReturn(0); 632 } 633 634 #undef __FUNCT__ 635 #define __FUNCT__ "MatMatSolve_MUMPS" 636 PetscErrorCode MatMatSolve_MUMPS(Mat A,Mat B,Mat X) 637 { 638 PetscErrorCode ierr; 639 PetscBool flg; 640 641 PetscFunctionBegin; 642 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 643 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 644 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 645 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 646 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_MUMPS() is not implemented yet"); 647 PetscFunctionReturn(0); 648 } 649 650 #if !defined(PETSC_USE_COMPLEX) 651 /* 652 input: 653 F: numeric factor 654 output: 655 nneg: total number of negative pivots 656 nzero: 0 657 npos: (global dimension of F) - nneg 658 */ 659 660 #undef __FUNCT__ 661 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 662 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 663 { 664 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 665 PetscErrorCode ierr; 666 PetscMPIInt size; 667 668 PetscFunctionBegin; 669 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);CHKERRQ(ierr); 670 /* 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 */ 671 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)); 672 if (nneg) { 673 if (!mumps->myid) { 674 *nneg = mumps->id.INFOG(12); 675 } 676 ierr = MPI_Bcast(nneg,1,MPI_INT,0,mumps->comm_mumps);CHKERRQ(ierr); 677 } 678 if (nzero) *nzero = 0; 679 if (npos) *npos = F->rmap->N - (*nneg); 680 PetscFunctionReturn(0); 681 } 682 #endif /* !defined(PETSC_USE_COMPLEX) */ 683 684 #undef __FUNCT__ 685 #define __FUNCT__ "MatFactorNumeric_MUMPS" 686 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 687 { 688 Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr; 689 PetscErrorCode ierr; 690 Mat F_diag; 691 PetscBool isMPIAIJ; 692 693 PetscFunctionBegin; 694 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 695 696 /* numerical factorization phase */ 697 /*-------------------------------*/ 698 mumps->id.job = JOB_FACTNUMERIC; 699 if (!mumps->id.ICNTL(18)) { 700 if (!mumps->myid) { 701 #if defined(PETSC_USE_COMPLEX) 702 #if defined(PETSC_USE_REAL_SINGLE) 703 mumps->id.a = (mumps_complex*)mumps->val; 704 #else 705 mumps->id.a = (mumps_double_complex*)mumps->val; 706 #endif 707 #else 708 mumps->id.a = mumps->val; 709 #endif 710 } 711 } else { 712 #if defined(PETSC_USE_COMPLEX) 713 #if defined(PETSC_USE_REAL_SINGLE) 714 mumps->id.a_loc = (mumps_complex*)mumps->val; 715 #else 716 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 717 #endif 718 #else 719 mumps->id.a_loc = mumps->val; 720 #endif 721 } 722 PetscMUMPS_c(&mumps->id); 723 if (mumps->id.INFOG(1) < 0) { 724 if (mumps->id.INFO(1) == -13) 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)); 725 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)); 726 } 727 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)); 728 729 if (mumps->size > 1) { 730 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 731 if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A; 732 else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A; 733 F_diag->assembled = PETSC_TRUE; 734 if (mumps->scat_sol) { 735 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 736 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 737 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 738 } 739 } 740 (F)->assembled = PETSC_TRUE; 741 mumps->matstruc = SAME_NONZERO_PATTERN; 742 mumps->CleanUpMUMPS = PETSC_TRUE; 743 744 if (mumps->size > 1) { 745 /* distributed solution */ 746 if (!mumps->scat_sol) { 747 /* Create x_seq=sol_loc for repeated use */ 748 PetscInt lsol_loc; 749 PetscScalar *sol_loc; 750 751 lsol_loc = mumps->id.INFO(23); /* length of sol_loc */ 752 753 ierr = PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&mumps->id.isol_loc);CHKERRQ(ierr); 754 755 mumps->id.lsol_loc = lsol_loc; 756 #if defined(PETSC_USE_COMPLEX) 757 #if defined(PETSC_USE_REAL_SINGLE) 758 mumps->id.sol_loc = (mumps_complex*)sol_loc; 759 #else 760 mumps->id.sol_loc = (mumps_double_complex*)sol_loc; 761 #endif 762 #else 763 mumps->id.sol_loc = sol_loc; 764 #endif 765 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 766 } 767 } 768 PetscFunctionReturn(0); 769 } 770 771 /* Sets MUMPS options from the options database */ 772 #undef __FUNCT__ 773 #define __FUNCT__ "PetscSetMUMPSFromOptions" 774 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 775 { 776 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 777 PetscErrorCode ierr; 778 PetscInt icntl; 779 PetscBool flg; 780 781 PetscFunctionBegin; 782 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 783 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 784 if (flg) mumps->id.ICNTL(1) = icntl; 785 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); 786 if (flg) mumps->id.ICNTL(2) = icntl; 787 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); 788 if (flg) mumps->id.ICNTL(3) = icntl; 789 790 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 791 if (flg) mumps->id.ICNTL(4) = icntl; 792 if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 793 794 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); 795 if (flg) mumps->id.ICNTL(6) = icntl; 796 797 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); 798 if (flg) { 799 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"); 800 else mumps->id.ICNTL(7) = icntl; 801 } 802 803 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); 804 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);CHKERRQ(ierr); 805 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); 806 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); 807 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); 808 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); 809 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);CHKERRQ(ierr); 810 811 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); 812 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); 813 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); 814 if (mumps->id.ICNTL(24)) { 815 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 816 } 817 818 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); 819 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); 820 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);CHKERRQ(ierr); 821 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); 822 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); 823 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); 824 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); 825 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);CHKERRQ(ierr); 826 827 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);CHKERRQ(ierr); 828 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);CHKERRQ(ierr); 829 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);CHKERRQ(ierr); 830 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);CHKERRQ(ierr); 831 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);CHKERRQ(ierr); 832 833 ierr = PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL); 834 PetscOptionsEnd(); 835 PetscFunctionReturn(0); 836 } 837 838 #undef __FUNCT__ 839 #define __FUNCT__ "PetscInitializeMUMPS" 840 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps) 841 { 842 PetscErrorCode ierr; 843 844 PetscFunctionBegin; 845 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid); 846 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);CHKERRQ(ierr); 847 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));CHKERRQ(ierr); 848 849 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 850 851 mumps->id.job = JOB_INIT; 852 mumps->id.par = 1; /* host participates factorizaton and solve */ 853 mumps->id.sym = mumps->sym; 854 PetscMUMPS_c(&mumps->id); 855 856 mumps->CleanUpMUMPS = PETSC_FALSE; 857 mumps->scat_rhs = NULL; 858 mumps->scat_sol = NULL; 859 860 /* set PETSc-MUMPS default options - override MUMPS default */ 861 mumps->id.ICNTL(3) = 0; 862 mumps->id.ICNTL(4) = 0; 863 if (mumps->size == 1) { 864 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 865 } else { 866 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 867 mumps->id.ICNTL(21) = 1; /* distributed solution */ 868 } 869 PetscFunctionReturn(0); 870 } 871 872 /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */ 873 #undef __FUNCT__ 874 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 875 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 876 { 877 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 878 PetscErrorCode ierr; 879 Vec b; 880 IS is_iden; 881 const PetscInt M = A->rmap->N; 882 883 PetscFunctionBegin; 884 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 885 886 /* Set MUMPS options from the options database */ 887 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 888 889 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 890 891 /* analysis phase */ 892 /*----------------*/ 893 mumps->id.job = JOB_FACTSYMBOLIC; 894 mumps->id.n = M; 895 switch (mumps->id.ICNTL(18)) { 896 case 0: /* centralized assembled matrix input */ 897 if (!mumps->myid) { 898 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 899 if (mumps->id.ICNTL(6)>1) { 900 #if defined(PETSC_USE_COMPLEX) 901 #if defined(PETSC_USE_REAL_SINGLE) 902 mumps->id.a = (mumps_complex*)mumps->val; 903 #else 904 mumps->id.a = (mumps_double_complex*)mumps->val; 905 #endif 906 #else 907 mumps->id.a = mumps->val; 908 #endif 909 } 910 if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */ 911 /* 912 PetscBool flag; 913 ierr = ISEqual(r,c,&flag);CHKERRQ(ierr); 914 if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm"); 915 ierr = ISView(r,PETSC_VIEWER_STDOUT_SELF); 916 */ 917 if (!mumps->myid) { 918 const PetscInt *idx; 919 PetscInt i,*perm_in; 920 921 ierr = PetscMalloc(M*sizeof(PetscInt),&perm_in);CHKERRQ(ierr); 922 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 923 924 mumps->id.perm_in = perm_in; 925 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 926 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 927 } 928 } 929 } 930 break; 931 case 3: /* distributed assembled matrix input (size>1) */ 932 mumps->id.nz_loc = mumps->nz; 933 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 934 if (mumps->id.ICNTL(6)>1) { 935 #if defined(PETSC_USE_COMPLEX) 936 #if defined(PETSC_USE_REAL_SINGLE) 937 mumps->id.a_loc = (mumps_complex*)mumps->val; 938 #else 939 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 940 #endif 941 #else 942 mumps->id.a_loc = mumps->val; 943 #endif 944 } 945 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 946 if (!mumps->myid) { 947 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 948 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 949 } else { 950 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 951 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 952 } 953 ierr = VecCreate(PetscObjectComm((PetscObject)A),&b);CHKERRQ(ierr); 954 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 955 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 956 957 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 958 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 959 ierr = VecDestroy(&b);CHKERRQ(ierr); 960 break; 961 } 962 PetscMUMPS_c(&mumps->id); 963 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)); 964 965 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 966 F->ops->solve = MatSolve_MUMPS; 967 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 968 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 969 PetscFunctionReturn(0); 970 } 971 972 /* Note the Petsc r and c permutations are ignored */ 973 #undef __FUNCT__ 974 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 975 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 976 { 977 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 978 PetscErrorCode ierr; 979 Vec b; 980 IS is_iden; 981 const PetscInt M = A->rmap->N; 982 983 PetscFunctionBegin; 984 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 985 986 /* Set MUMPS options from the options database */ 987 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 988 989 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 990 991 /* analysis phase */ 992 /*----------------*/ 993 mumps->id.job = JOB_FACTSYMBOLIC; 994 mumps->id.n = M; 995 switch (mumps->id.ICNTL(18)) { 996 case 0: /* centralized assembled matrix input */ 997 if (!mumps->myid) { 998 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 999 if (mumps->id.ICNTL(6)>1) { 1000 #if defined(PETSC_USE_COMPLEX) 1001 #if defined(PETSC_USE_REAL_SINGLE) 1002 mumps->id.a = (mumps_complex*)mumps->val; 1003 #else 1004 mumps->id.a = (mumps_double_complex*)mumps->val; 1005 #endif 1006 #else 1007 mumps->id.a = mumps->val; 1008 #endif 1009 } 1010 } 1011 break; 1012 case 3: /* distributed assembled matrix input (size>1) */ 1013 mumps->id.nz_loc = mumps->nz; 1014 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1015 if (mumps->id.ICNTL(6)>1) { 1016 #if defined(PETSC_USE_COMPLEX) 1017 #if defined(PETSC_USE_REAL_SINGLE) 1018 mumps->id.a_loc = (mumps_complex*)mumps->val; 1019 #else 1020 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1021 #endif 1022 #else 1023 mumps->id.a_loc = mumps->val; 1024 #endif 1025 } 1026 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1027 if (!mumps->myid) { 1028 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1029 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1030 } else { 1031 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1032 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1033 } 1034 ierr = VecCreate(PetscObjectComm((PetscObject)A),&b);CHKERRQ(ierr); 1035 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 1036 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 1037 1038 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1039 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1040 ierr = VecDestroy(&b);CHKERRQ(ierr); 1041 break; 1042 } 1043 PetscMUMPS_c(&mumps->id); 1044 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)); 1045 1046 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1047 F->ops->solve = MatSolve_MUMPS; 1048 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1049 PetscFunctionReturn(0); 1050 } 1051 1052 /* Note the Petsc r permutation and factor info are ignored */ 1053 #undef __FUNCT__ 1054 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 1055 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 1056 { 1057 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1058 PetscErrorCode ierr; 1059 Vec b; 1060 IS is_iden; 1061 const PetscInt M = A->rmap->N; 1062 1063 PetscFunctionBegin; 1064 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1065 1066 /* Set MUMPS options from the options database */ 1067 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1068 1069 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1070 1071 /* analysis phase */ 1072 /*----------------*/ 1073 mumps->id.job = JOB_FACTSYMBOLIC; 1074 mumps->id.n = M; 1075 switch (mumps->id.ICNTL(18)) { 1076 case 0: /* centralized assembled matrix input */ 1077 if (!mumps->myid) { 1078 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1079 if (mumps->id.ICNTL(6)>1) { 1080 #if defined(PETSC_USE_COMPLEX) 1081 #if defined(PETSC_USE_REAL_SINGLE) 1082 mumps->id.a = (mumps_complex*)mumps->val; 1083 #else 1084 mumps->id.a = (mumps_double_complex*)mumps->val; 1085 #endif 1086 #else 1087 mumps->id.a = mumps->val; 1088 #endif 1089 } 1090 } 1091 break; 1092 case 3: /* distributed assembled matrix input (size>1) */ 1093 mumps->id.nz_loc = mumps->nz; 1094 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1095 if (mumps->id.ICNTL(6)>1) { 1096 #if defined(PETSC_USE_COMPLEX) 1097 #if defined(PETSC_USE_REAL_SINGLE) 1098 mumps->id.a_loc = (mumps_complex*)mumps->val; 1099 #else 1100 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1101 #endif 1102 #else 1103 mumps->id.a_loc = mumps->val; 1104 #endif 1105 } 1106 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1107 if (!mumps->myid) { 1108 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1109 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1110 } else { 1111 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1112 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1113 } 1114 ierr = VecCreate(PetscObjectComm((PetscObject)A),&b);CHKERRQ(ierr); 1115 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 1116 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 1117 1118 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1119 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1120 ierr = VecDestroy(&b);CHKERRQ(ierr); 1121 break; 1122 } 1123 PetscMUMPS_c(&mumps->id); 1124 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)); 1125 1126 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1127 F->ops->solve = MatSolve_MUMPS; 1128 F->ops->solvetranspose = MatSolve_MUMPS; 1129 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 1130 #if !defined(PETSC_USE_COMPLEX) 1131 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1132 #else 1133 F->ops->getinertia = NULL; 1134 #endif 1135 PetscFunctionReturn(0); 1136 } 1137 1138 #undef __FUNCT__ 1139 #define __FUNCT__ "MatView_MUMPS" 1140 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1141 { 1142 PetscErrorCode ierr; 1143 PetscBool iascii; 1144 PetscViewerFormat format; 1145 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1146 1147 PetscFunctionBegin; 1148 /* check if matrix is mumps type */ 1149 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1150 1151 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1152 if (iascii) { 1153 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1154 if (format == PETSC_VIEWER_ASCII_INFO) { 1155 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1156 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1157 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1158 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));CHKERRQ(ierr); 1159 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));CHKERRQ(ierr); 1160 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));CHKERRQ(ierr); 1161 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));CHKERRQ(ierr); 1162 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));CHKERRQ(ierr); 1163 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));CHKERRQ(ierr); 1164 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));CHKERRQ(ierr); 1165 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));CHKERRQ(ierr); 1166 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));CHKERRQ(ierr); 1167 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));CHKERRQ(ierr); 1168 if (mumps->id.ICNTL(11)>0) { 1169 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));CHKERRQ(ierr); 1170 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));CHKERRQ(ierr); 1171 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));CHKERRQ(ierr); 1172 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));CHKERRQ(ierr); 1173 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));CHKERRQ(ierr); 1174 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));CHKERRQ(ierr); 1175 } 1176 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));CHKERRQ(ierr); 1177 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));CHKERRQ(ierr); 1178 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));CHKERRQ(ierr); 1179 /* ICNTL(15-17) not used */ 1180 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));CHKERRQ(ierr); 1181 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));CHKERRQ(ierr); 1182 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));CHKERRQ(ierr); 1183 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));CHKERRQ(ierr); 1184 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));CHKERRQ(ierr); 1185 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));CHKERRQ(ierr); 1186 1187 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));CHKERRQ(ierr); 1188 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));CHKERRQ(ierr); 1189 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));CHKERRQ(ierr); 1190 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));CHKERRQ(ierr); 1191 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));CHKERRQ(ierr); 1192 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));CHKERRQ(ierr); 1193 1194 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));CHKERRQ(ierr); 1195 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));CHKERRQ(ierr); 1196 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));CHKERRQ(ierr); 1197 1198 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));CHKERRQ(ierr); 1199 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));CHKERRQ(ierr); 1200 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));CHKERRQ(ierr); 1201 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));CHKERRQ(ierr); 1202 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));CHKERRQ(ierr); 1203 1204 /* infomation local to each processor */ 1205 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1206 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1207 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));CHKERRQ(ierr); 1208 ierr = PetscViewerFlush(viewer); 1209 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1210 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));CHKERRQ(ierr); 1211 ierr = PetscViewerFlush(viewer); 1212 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1213 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));CHKERRQ(ierr); 1214 ierr = PetscViewerFlush(viewer); 1215 1216 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1217 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));CHKERRQ(ierr); 1218 ierr = PetscViewerFlush(viewer); 1219 1220 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1221 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));CHKERRQ(ierr); 1222 ierr = PetscViewerFlush(viewer); 1223 1224 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1225 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));CHKERRQ(ierr); 1226 ierr = PetscViewerFlush(viewer); 1227 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1228 1229 if (!mumps->myid) { /* information from the host */ 1230 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));CHKERRQ(ierr); 1231 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));CHKERRQ(ierr); 1232 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));CHKERRQ(ierr); 1233 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); 1234 1235 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));CHKERRQ(ierr); 1236 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));CHKERRQ(ierr); 1237 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));CHKERRQ(ierr); 1238 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));CHKERRQ(ierr); 1239 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));CHKERRQ(ierr); 1240 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));CHKERRQ(ierr); 1241 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));CHKERRQ(ierr); 1242 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));CHKERRQ(ierr); 1243 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));CHKERRQ(ierr); 1244 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));CHKERRQ(ierr); 1245 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));CHKERRQ(ierr); 1246 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));CHKERRQ(ierr); 1247 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));CHKERRQ(ierr); 1248 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); 1249 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); 1250 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); 1251 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); 1252 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));CHKERRQ(ierr); 1253 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); 1254 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); 1255 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));CHKERRQ(ierr); 1256 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));CHKERRQ(ierr); 1257 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));CHKERRQ(ierr); 1258 } 1259 } 1260 } 1261 PetscFunctionReturn(0); 1262 } 1263 1264 #undef __FUNCT__ 1265 #define __FUNCT__ "MatGetInfo_MUMPS" 1266 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1267 { 1268 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1269 1270 PetscFunctionBegin; 1271 info->block_size = 1.0; 1272 info->nz_allocated = mumps->id.INFOG(20); 1273 info->nz_used = mumps->id.INFOG(20); 1274 info->nz_unneeded = 0.0; 1275 info->assemblies = 0.0; 1276 info->mallocs = 0.0; 1277 info->memory = 0.0; 1278 info->fill_ratio_given = 0; 1279 info->fill_ratio_needed = 0; 1280 info->factor_mallocs = 0; 1281 PetscFunctionReturn(0); 1282 } 1283 1284 /* -------------------------------------------------------------------------------------------*/ 1285 #undef __FUNCT__ 1286 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1287 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1288 { 1289 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1290 1291 PetscFunctionBegin; 1292 mumps->id.ICNTL(icntl) = ival; 1293 PetscFunctionReturn(0); 1294 } 1295 1296 #undef __FUNCT__ 1297 #define __FUNCT__ "MatMumpsSetIcntl" 1298 /*@ 1299 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1300 1301 Logically Collective on Mat 1302 1303 Input Parameters: 1304 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1305 . icntl - index of MUMPS parameter array ICNTL() 1306 - ival - value of MUMPS ICNTL(icntl) 1307 1308 Options Database: 1309 . -mat_mumps_icntl_<icntl> <ival> 1310 1311 Level: beginner 1312 1313 References: MUMPS Users' Guide 1314 1315 .seealso: MatGetFactor() 1316 @*/ 1317 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1318 { 1319 PetscErrorCode ierr; 1320 1321 PetscFunctionBegin; 1322 PetscValidLogicalCollectiveInt(F,icntl,2); 1323 PetscValidLogicalCollectiveInt(F,ival,3); 1324 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1325 PetscFunctionReturn(0); 1326 } 1327 1328 /* -------------------------------------------------------------------------------------------*/ 1329 #undef __FUNCT__ 1330 #define __FUNCT__ "MatMumpsSetCntl_MUMPS" 1331 PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val) 1332 { 1333 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1334 1335 PetscFunctionBegin; 1336 mumps->id.CNTL(icntl) = val; 1337 PetscFunctionReturn(0); 1338 } 1339 1340 #undef __FUNCT__ 1341 #define __FUNCT__ "MatMumpsSetCntl" 1342 /*@ 1343 MatMumpsSetCntl - Set MUMPS parameter CNTL() 1344 1345 Logically Collective on Mat 1346 1347 Input Parameters: 1348 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1349 . icntl - index of MUMPS parameter array CNTL() 1350 - val - value of MUMPS CNTL(icntl) 1351 1352 Options Database: 1353 . -mat_mumps_cntl_<icntl> <val> 1354 1355 Level: beginner 1356 1357 References: MUMPS Users' Guide 1358 1359 .seealso: MatGetFactor() 1360 @*/ 1361 PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val) 1362 { 1363 PetscErrorCode ierr; 1364 1365 PetscFunctionBegin; 1366 PetscValidLogicalCollectiveInt(F,icntl,2); 1367 PetscValidLogicalCollectiveInt(F,val,3); 1368 ierr = PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));CHKERRQ(ierr); 1369 PetscFunctionReturn(0); 1370 } 1371 1372 /*MC 1373 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1374 distributed and sequential matrices via the external package MUMPS. 1375 1376 Works with MATAIJ and MATSBAIJ matrices 1377 1378 Options Database Keys: 1379 + -mat_mumps_icntl_4 <0,...,4> - print level 1380 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1381 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec) 1382 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1383 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1384 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1385 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1386 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1387 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1388 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1389 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1390 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1391 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1392 1393 Level: beginner 1394 1395 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1396 1397 M*/ 1398 1399 #undef __FUNCT__ 1400 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1401 static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1402 { 1403 PetscFunctionBegin; 1404 *type = MATSOLVERMUMPS; 1405 PetscFunctionReturn(0); 1406 } 1407 1408 /* MatGetFactor for Seq and MPI AIJ matrices */ 1409 #undef __FUNCT__ 1410 #define __FUNCT__ "MatGetFactor_aij_mumps" 1411 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1412 { 1413 Mat B; 1414 PetscErrorCode ierr; 1415 Mat_MUMPS *mumps; 1416 PetscBool isSeqAIJ; 1417 1418 PetscFunctionBegin; 1419 /* Create the factorization matrix */ 1420 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1421 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1422 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1423 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1424 if (isSeqAIJ) { 1425 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1426 } else { 1427 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1428 } 1429 1430 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1431 1432 B->ops->view = MatView_MUMPS; 1433 B->ops->getinfo = MatGetInfo_MUMPS; 1434 1435 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1436 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1437 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1438 if (ftype == MAT_FACTOR_LU) { 1439 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1440 B->factortype = MAT_FACTOR_LU; 1441 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1442 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1443 mumps->sym = 0; 1444 } else { 1445 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1446 B->factortype = MAT_FACTOR_CHOLESKY; 1447 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1448 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1449 if (A->spd_set && A->spd) mumps->sym = 1; 1450 else mumps->sym = 2; 1451 } 1452 1453 mumps->isAIJ = PETSC_TRUE; 1454 mumps->Destroy = B->ops->destroy; 1455 B->ops->destroy = MatDestroy_MUMPS; 1456 B->spptr = (void*)mumps; 1457 1458 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1459 1460 *F = B; 1461 PetscFunctionReturn(0); 1462 } 1463 1464 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1465 #undef __FUNCT__ 1466 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1467 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1468 { 1469 Mat B; 1470 PetscErrorCode ierr; 1471 Mat_MUMPS *mumps; 1472 PetscBool isSeqSBAIJ; 1473 1474 PetscFunctionBegin; 1475 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1476 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"); 1477 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1478 /* Create the factorization matrix */ 1479 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1480 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1481 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1482 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1483 if (isSeqSBAIJ) { 1484 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 1485 1486 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1487 } else { 1488 ierr = MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);CHKERRQ(ierr); 1489 1490 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1491 } 1492 1493 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1494 B->ops->view = MatView_MUMPS; 1495 1496 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1497 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl);CHKERRQ(ierr); 1498 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl);CHKERRQ(ierr); 1499 1500 B->factortype = MAT_FACTOR_CHOLESKY; 1501 if (A->spd_set && A->spd) mumps->sym = 1; 1502 else mumps->sym = 2; 1503 1504 mumps->isAIJ = PETSC_FALSE; 1505 mumps->Destroy = B->ops->destroy; 1506 B->ops->destroy = MatDestroy_MUMPS; 1507 B->spptr = (void*)mumps; 1508 1509 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1510 1511 *F = B; 1512 PetscFunctionReturn(0); 1513 } 1514 1515 #undef __FUNCT__ 1516 #define __FUNCT__ "MatGetFactor_baij_mumps" 1517 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 1518 { 1519 Mat B; 1520 PetscErrorCode ierr; 1521 Mat_MUMPS *mumps; 1522 PetscBool isSeqBAIJ; 1523 1524 PetscFunctionBegin; 1525 /* Create the factorization matrix */ 1526 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1527 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1528 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1529 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1530 if (isSeqBAIJ) { 1531 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);CHKERRQ(ierr); 1532 } else { 1533 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 1534 } 1535 1536 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1537 if (ftype == MAT_FACTOR_LU) { 1538 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1539 B->factortype = MAT_FACTOR_LU; 1540 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 1541 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 1542 mumps->sym = 0; 1543 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 1544 1545 B->ops->view = MatView_MUMPS; 1546 1547 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1548 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1549 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1550 1551 mumps->isAIJ = PETSC_TRUE; 1552 mumps->Destroy = B->ops->destroy; 1553 B->ops->destroy = MatDestroy_MUMPS; 1554 B->spptr = (void*)mumps; 1555 1556 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1557 1558 *F = B; 1559 PetscFunctionReturn(0); 1560 } 1561 1562