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 560 PetscFunctionBegin; 561 mumps->id.nrhs = 1; 562 b_seq = mumps->b_seq; 563 if (mumps->size > 1) { 564 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 565 ierr = VecScatterBegin(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 566 ierr = VecScatterEnd(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 567 if (!mumps->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 568 } else { /* size == 1 */ 569 ierr = VecCopy(b,x);CHKERRQ(ierr); 570 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 571 } 572 if (!mumps->myid) { /* define rhs on the host */ 573 mumps->id.nrhs = 1; 574 #if defined(PETSC_USE_COMPLEX) 575 #if defined(PETSC_USE_REAL_SINGLE) 576 mumps->id.rhs = (mumps_complex*)array; 577 #else 578 mumps->id.rhs = (mumps_double_complex*)array; 579 #endif 580 #else 581 mumps->id.rhs = array; 582 #endif 583 } 584 585 /* solve phase */ 586 /*-------------*/ 587 mumps->id.job = JOB_SOLVE; 588 PetscMUMPS_c(&mumps->id); 589 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)); 590 591 if (mumps->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 592 if (mumps->scat_sol && mumps->ICNTL9_pre != mumps->id.ICNTL(9)) { 593 /* when id.ICNTL(9) changes, the contents of lsol_loc may change (not its size, lsol_loc), recreates scat_sol */ 594 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 595 } 596 if (!mumps->scat_sol) { /* create scatter scat_sol */ 597 ierr = ISCreateStride(PETSC_COMM_SELF,mumps->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 598 for (i=0; i<mumps->id.lsol_loc; i++) { 599 mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 600 } 601 ierr = ISCreateGeneral(PETSC_COMM_SELF,mumps->id.lsol_loc,mumps->id.isol_loc,PETSC_COPY_VALUES,&is_petsc);CHKERRQ(ierr); /* to */ 602 ierr = VecScatterCreate(mumps->x_seq,is_iden,x,is_petsc,&mumps->scat_sol);CHKERRQ(ierr); 603 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 604 ierr = ISDestroy(&is_petsc);CHKERRQ(ierr); 605 606 mumps->ICNTL9_pre = mumps->id.ICNTL(9); /* save current value of id.ICNTL(9) */ 607 } 608 609 ierr = VecScatterBegin(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 610 ierr = VecScatterEnd(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 611 } 612 PetscFunctionReturn(0); 613 } 614 615 #undef __FUNCT__ 616 #define __FUNCT__ "MatSolveTranspose_MUMPS" 617 PetscErrorCode MatSolveTranspose_MUMPS(Mat A,Vec b,Vec x) 618 { 619 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 620 PetscErrorCode ierr; 621 622 PetscFunctionBegin; 623 mumps->id.ICNTL(9) = 0; 624 625 ierr = MatSolve_MUMPS(A,b,x);CHKERRQ(ierr); 626 627 mumps->id.ICNTL(9) = 1; 628 PetscFunctionReturn(0); 629 } 630 631 #undef __FUNCT__ 632 #define __FUNCT__ "MatMatSolve_MUMPS" 633 PetscErrorCode MatMatSolve_MUMPS(Mat A,Mat B,Mat X) 634 { 635 PetscErrorCode ierr; 636 PetscBool flg; 637 638 PetscFunctionBegin; 639 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 640 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 641 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 642 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 643 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_MUMPS() is not implemented yet"); 644 PetscFunctionReturn(0); 645 } 646 647 #if !defined(PETSC_USE_COMPLEX) 648 /* 649 input: 650 F: numeric factor 651 output: 652 nneg: total number of negative pivots 653 nzero: 0 654 npos: (global dimension of F) - nneg 655 */ 656 657 #undef __FUNCT__ 658 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 659 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 660 { 661 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 662 PetscErrorCode ierr; 663 PetscMPIInt size; 664 665 PetscFunctionBegin; 666 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);CHKERRQ(ierr); 667 /* 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 */ 668 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)); 669 if (nneg) { 670 if (!mumps->myid) { 671 *nneg = mumps->id.INFOG(12); 672 } 673 ierr = MPI_Bcast(nneg,1,MPI_INT,0,mumps->comm_mumps);CHKERRQ(ierr); 674 } 675 if (nzero) *nzero = 0; 676 if (npos) *npos = F->rmap->N - (*nneg); 677 PetscFunctionReturn(0); 678 } 679 #endif /* !defined(PETSC_USE_COMPLEX) */ 680 681 #undef __FUNCT__ 682 #define __FUNCT__ "MatFactorNumeric_MUMPS" 683 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 684 { 685 Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr; 686 PetscErrorCode ierr; 687 Mat F_diag; 688 PetscBool isMPIAIJ; 689 690 PetscFunctionBegin; 691 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 692 693 /* numerical factorization phase */ 694 /*-------------------------------*/ 695 mumps->id.job = JOB_FACTNUMERIC; 696 if (!mumps->id.ICNTL(18)) { 697 if (!mumps->myid) { 698 #if defined(PETSC_USE_COMPLEX) 699 #if defined(PETSC_USE_REAL_SINGLE) 700 mumps->id.a = (mumps_complex*)mumps->val; 701 #else 702 mumps->id.a = (mumps_double_complex*)mumps->val; 703 #endif 704 #else 705 mumps->id.a = mumps->val; 706 #endif 707 } 708 } else { 709 #if defined(PETSC_USE_COMPLEX) 710 #if defined(PETSC_USE_REAL_SINGLE) 711 mumps->id.a_loc = (mumps_complex*)mumps->val; 712 #else 713 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 714 #endif 715 #else 716 mumps->id.a_loc = mumps->val; 717 #endif 718 } 719 PetscMUMPS_c(&mumps->id); 720 if (mumps->id.INFOG(1) < 0) { 721 if (mumps->id.INFO(1) == -13) { 722 if (mumps->id.INFO(2) < 0) { 723 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)); 724 } else { 725 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)); 726 } 727 } 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)); 728 } 729 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)); 730 731 if (mumps->size > 1) { 732 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 733 if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A; 734 else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A; 735 F_diag->assembled = PETSC_TRUE; 736 if (mumps->scat_sol) { 737 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 738 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 739 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 740 } 741 } 742 (F)->assembled = PETSC_TRUE; 743 mumps->matstruc = SAME_NONZERO_PATTERN; 744 mumps->CleanUpMUMPS = PETSC_TRUE; 745 746 if (mumps->size > 1) { 747 /* distributed solution */ 748 if (!mumps->scat_sol) { 749 /* Create x_seq=sol_loc for repeated use */ 750 PetscInt lsol_loc; 751 PetscScalar *sol_loc; 752 753 lsol_loc = mumps->id.INFO(23); /* length of sol_loc */ 754 755 ierr = PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&mumps->id.isol_loc);CHKERRQ(ierr); 756 757 mumps->id.lsol_loc = lsol_loc; 758 #if defined(PETSC_USE_COMPLEX) 759 #if defined(PETSC_USE_REAL_SINGLE) 760 mumps->id.sol_loc = (mumps_complex*)sol_loc; 761 #else 762 mumps->id.sol_loc = (mumps_double_complex*)sol_loc; 763 #endif 764 #else 765 mumps->id.sol_loc = sol_loc; 766 #endif 767 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 768 } 769 } 770 PetscFunctionReturn(0); 771 } 772 773 /* Sets MUMPS options from the options database */ 774 #undef __FUNCT__ 775 #define __FUNCT__ "PetscSetMUMPSFromOptions" 776 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 777 { 778 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 779 PetscErrorCode ierr; 780 PetscInt icntl; 781 PetscBool flg; 782 783 PetscFunctionBegin; 784 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 785 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 786 if (flg) mumps->id.ICNTL(1) = icntl; 787 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); 788 if (flg) mumps->id.ICNTL(2) = icntl; 789 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); 790 if (flg) mumps->id.ICNTL(3) = icntl; 791 792 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 793 if (flg) mumps->id.ICNTL(4) = icntl; 794 if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 795 796 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); 797 if (flg) mumps->id.ICNTL(6) = icntl; 798 799 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); 800 if (flg) { 801 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"); 802 else mumps->id.ICNTL(7) = icntl; 803 } 804 805 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); 806 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);CHKERRQ(ierr); 807 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); 808 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); 809 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); 810 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); 811 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);CHKERRQ(ierr); 812 813 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); 814 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); 815 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); 816 if (mumps->id.ICNTL(24)) { 817 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 818 } 819 820 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); 821 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); 822 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);CHKERRQ(ierr); 823 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); 824 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); 825 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); 826 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); 827 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);CHKERRQ(ierr); 828 829 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);CHKERRQ(ierr); 830 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);CHKERRQ(ierr); 831 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);CHKERRQ(ierr); 832 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);CHKERRQ(ierr); 833 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);CHKERRQ(ierr); 834 835 ierr = PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL); 836 PetscOptionsEnd(); 837 PetscFunctionReturn(0); 838 } 839 840 #undef __FUNCT__ 841 #define __FUNCT__ "PetscInitializeMUMPS" 842 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps) 843 { 844 PetscErrorCode ierr; 845 846 PetscFunctionBegin; 847 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid); 848 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);CHKERRQ(ierr); 849 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));CHKERRQ(ierr); 850 851 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 852 853 mumps->id.job = JOB_INIT; 854 mumps->id.par = 1; /* host participates factorizaton and solve */ 855 mumps->id.sym = mumps->sym; 856 PetscMUMPS_c(&mumps->id); 857 858 mumps->CleanUpMUMPS = PETSC_FALSE; 859 mumps->scat_rhs = NULL; 860 mumps->scat_sol = NULL; 861 862 /* set PETSc-MUMPS default options - override MUMPS default */ 863 mumps->id.ICNTL(3) = 0; 864 mumps->id.ICNTL(4) = 0; 865 if (mumps->size == 1) { 866 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 867 } else { 868 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 869 mumps->id.ICNTL(21) = 1; /* distributed solution */ 870 } 871 PetscFunctionReturn(0); 872 } 873 874 /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */ 875 #undef __FUNCT__ 876 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 877 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 878 { 879 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 880 PetscErrorCode ierr; 881 Vec b; 882 IS is_iden; 883 const PetscInt M = A->rmap->N; 884 885 PetscFunctionBegin; 886 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 887 888 /* Set MUMPS options from the options database */ 889 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 890 891 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 892 893 /* analysis phase */ 894 /*----------------*/ 895 mumps->id.job = JOB_FACTSYMBOLIC; 896 mumps->id.n = M; 897 switch (mumps->id.ICNTL(18)) { 898 case 0: /* centralized assembled matrix input */ 899 if (!mumps->myid) { 900 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 901 if (mumps->id.ICNTL(6)>1) { 902 #if defined(PETSC_USE_COMPLEX) 903 #if defined(PETSC_USE_REAL_SINGLE) 904 mumps->id.a = (mumps_complex*)mumps->val; 905 #else 906 mumps->id.a = (mumps_double_complex*)mumps->val; 907 #endif 908 #else 909 mumps->id.a = mumps->val; 910 #endif 911 } 912 if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */ 913 /* 914 PetscBool flag; 915 ierr = ISEqual(r,c,&flag);CHKERRQ(ierr); 916 if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm"); 917 ierr = ISView(r,PETSC_VIEWER_STDOUT_SELF); 918 */ 919 if (!mumps->myid) { 920 const PetscInt *idx; 921 PetscInt i,*perm_in; 922 923 ierr = PetscMalloc(M*sizeof(PetscInt),&perm_in);CHKERRQ(ierr); 924 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 925 926 mumps->id.perm_in = perm_in; 927 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 928 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 929 } 930 } 931 } 932 break; 933 case 3: /* distributed assembled matrix input (size>1) */ 934 mumps->id.nz_loc = mumps->nz; 935 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 936 if (mumps->id.ICNTL(6)>1) { 937 #if defined(PETSC_USE_COMPLEX) 938 #if defined(PETSC_USE_REAL_SINGLE) 939 mumps->id.a_loc = (mumps_complex*)mumps->val; 940 #else 941 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 942 #endif 943 #else 944 mumps->id.a_loc = mumps->val; 945 #endif 946 } 947 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 948 if (!mumps->myid) { 949 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 950 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 951 } else { 952 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 953 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 954 } 955 ierr = VecCreate(PetscObjectComm((PetscObject)A),&b);CHKERRQ(ierr); 956 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 957 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 958 959 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 960 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 961 ierr = VecDestroy(&b);CHKERRQ(ierr); 962 break; 963 } 964 PetscMUMPS_c(&mumps->id); 965 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)); 966 967 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 968 F->ops->solve = MatSolve_MUMPS; 969 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 970 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 971 PetscFunctionReturn(0); 972 } 973 974 /* Note the Petsc r and c permutations are ignored */ 975 #undef __FUNCT__ 976 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 977 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 978 { 979 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 980 PetscErrorCode ierr; 981 Vec b; 982 IS is_iden; 983 const PetscInt M = A->rmap->N; 984 985 PetscFunctionBegin; 986 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 987 988 /* Set MUMPS options from the options database */ 989 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 990 991 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 992 993 /* analysis phase */ 994 /*----------------*/ 995 mumps->id.job = JOB_FACTSYMBOLIC; 996 mumps->id.n = M; 997 switch (mumps->id.ICNTL(18)) { 998 case 0: /* centralized assembled matrix input */ 999 if (!mumps->myid) { 1000 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=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 = (mumps_complex*)mumps->val; 1005 #else 1006 mumps->id.a = (mumps_double_complex*)mumps->val; 1007 #endif 1008 #else 1009 mumps->id.a = mumps->val; 1010 #endif 1011 } 1012 } 1013 break; 1014 case 3: /* distributed assembled matrix input (size>1) */ 1015 mumps->id.nz_loc = mumps->nz; 1016 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1017 if (mumps->id.ICNTL(6)>1) { 1018 #if defined(PETSC_USE_COMPLEX) 1019 #if defined(PETSC_USE_REAL_SINGLE) 1020 mumps->id.a_loc = (mumps_complex*)mumps->val; 1021 #else 1022 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1023 #endif 1024 #else 1025 mumps->id.a_loc = mumps->val; 1026 #endif 1027 } 1028 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1029 if (!mumps->myid) { 1030 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1031 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1032 } else { 1033 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1034 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1035 } 1036 ierr = VecCreate(PetscObjectComm((PetscObject)A),&b);CHKERRQ(ierr); 1037 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 1038 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 1039 1040 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1041 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1042 ierr = VecDestroy(&b);CHKERRQ(ierr); 1043 break; 1044 } 1045 PetscMUMPS_c(&mumps->id); 1046 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)); 1047 1048 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1049 F->ops->solve = MatSolve_MUMPS; 1050 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1051 PetscFunctionReturn(0); 1052 } 1053 1054 /* Note the Petsc r permutation and factor info are ignored */ 1055 #undef __FUNCT__ 1056 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 1057 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 1058 { 1059 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1060 PetscErrorCode ierr; 1061 Vec b; 1062 IS is_iden; 1063 const PetscInt M = A->rmap->N; 1064 1065 PetscFunctionBegin; 1066 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1067 1068 /* Set MUMPS options from the options database */ 1069 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1070 1071 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1072 1073 /* analysis phase */ 1074 /*----------------*/ 1075 mumps->id.job = JOB_FACTSYMBOLIC; 1076 mumps->id.n = M; 1077 switch (mumps->id.ICNTL(18)) { 1078 case 0: /* centralized assembled matrix input */ 1079 if (!mumps->myid) { 1080 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1081 if (mumps->id.ICNTL(6)>1) { 1082 #if defined(PETSC_USE_COMPLEX) 1083 #if defined(PETSC_USE_REAL_SINGLE) 1084 mumps->id.a = (mumps_complex*)mumps->val; 1085 #else 1086 mumps->id.a = (mumps_double_complex*)mumps->val; 1087 #endif 1088 #else 1089 mumps->id.a = mumps->val; 1090 #endif 1091 } 1092 } 1093 break; 1094 case 3: /* distributed assembled matrix input (size>1) */ 1095 mumps->id.nz_loc = mumps->nz; 1096 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1097 if (mumps->id.ICNTL(6)>1) { 1098 #if defined(PETSC_USE_COMPLEX) 1099 #if defined(PETSC_USE_REAL_SINGLE) 1100 mumps->id.a_loc = (mumps_complex*)mumps->val; 1101 #else 1102 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1103 #endif 1104 #else 1105 mumps->id.a_loc = mumps->val; 1106 #endif 1107 } 1108 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1109 if (!mumps->myid) { 1110 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1111 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1112 } else { 1113 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1114 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1115 } 1116 ierr = VecCreate(PetscObjectComm((PetscObject)A),&b);CHKERRQ(ierr); 1117 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 1118 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 1119 1120 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1121 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1122 ierr = VecDestroy(&b);CHKERRQ(ierr); 1123 break; 1124 } 1125 PetscMUMPS_c(&mumps->id); 1126 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)); 1127 1128 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1129 F->ops->solve = MatSolve_MUMPS; 1130 F->ops->solvetranspose = MatSolve_MUMPS; 1131 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 1132 #if !defined(PETSC_USE_COMPLEX) 1133 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1134 #else 1135 F->ops->getinertia = NULL; 1136 #endif 1137 PetscFunctionReturn(0); 1138 } 1139 1140 #undef __FUNCT__ 1141 #define __FUNCT__ "MatView_MUMPS" 1142 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1143 { 1144 PetscErrorCode ierr; 1145 PetscBool iascii; 1146 PetscViewerFormat format; 1147 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1148 1149 PetscFunctionBegin; 1150 /* check if matrix is mumps type */ 1151 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1152 1153 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1154 if (iascii) { 1155 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1156 if (format == PETSC_VIEWER_ASCII_INFO) { 1157 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1158 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1159 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1160 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));CHKERRQ(ierr); 1161 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));CHKERRQ(ierr); 1162 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));CHKERRQ(ierr); 1163 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));CHKERRQ(ierr); 1164 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));CHKERRQ(ierr); 1165 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));CHKERRQ(ierr); 1166 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));CHKERRQ(ierr); 1167 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));CHKERRQ(ierr); 1168 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));CHKERRQ(ierr); 1169 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));CHKERRQ(ierr); 1170 if (mumps->id.ICNTL(11)>0) { 1171 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));CHKERRQ(ierr); 1172 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));CHKERRQ(ierr); 1173 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));CHKERRQ(ierr); 1174 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));CHKERRQ(ierr); 1175 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));CHKERRQ(ierr); 1176 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));CHKERRQ(ierr); 1177 } 1178 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));CHKERRQ(ierr); 1179 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));CHKERRQ(ierr); 1180 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));CHKERRQ(ierr); 1181 /* ICNTL(15-17) not used */ 1182 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));CHKERRQ(ierr); 1183 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));CHKERRQ(ierr); 1184 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));CHKERRQ(ierr); 1185 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));CHKERRQ(ierr); 1186 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));CHKERRQ(ierr); 1187 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));CHKERRQ(ierr); 1188 1189 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));CHKERRQ(ierr); 1190 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));CHKERRQ(ierr); 1191 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));CHKERRQ(ierr); 1192 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));CHKERRQ(ierr); 1193 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));CHKERRQ(ierr); 1194 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));CHKERRQ(ierr); 1195 1196 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));CHKERRQ(ierr); 1197 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));CHKERRQ(ierr); 1198 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));CHKERRQ(ierr); 1199 1200 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));CHKERRQ(ierr); 1201 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));CHKERRQ(ierr); 1202 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));CHKERRQ(ierr); 1203 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));CHKERRQ(ierr); 1204 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));CHKERRQ(ierr); 1205 1206 /* infomation local to each processor */ 1207 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1208 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1209 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));CHKERRQ(ierr); 1210 ierr = PetscViewerFlush(viewer); 1211 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1212 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));CHKERRQ(ierr); 1213 ierr = PetscViewerFlush(viewer); 1214 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1215 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));CHKERRQ(ierr); 1216 ierr = PetscViewerFlush(viewer); 1217 1218 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1219 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));CHKERRQ(ierr); 1220 ierr = PetscViewerFlush(viewer); 1221 1222 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1223 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));CHKERRQ(ierr); 1224 ierr = PetscViewerFlush(viewer); 1225 1226 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1227 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));CHKERRQ(ierr); 1228 ierr = PetscViewerFlush(viewer); 1229 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1230 1231 if (!mumps->myid) { /* information from the host */ 1232 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));CHKERRQ(ierr); 1233 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));CHKERRQ(ierr); 1234 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));CHKERRQ(ierr); 1235 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); 1236 1237 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));CHKERRQ(ierr); 1238 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));CHKERRQ(ierr); 1239 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));CHKERRQ(ierr); 1240 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));CHKERRQ(ierr); 1241 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));CHKERRQ(ierr); 1242 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));CHKERRQ(ierr); 1243 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));CHKERRQ(ierr); 1244 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));CHKERRQ(ierr); 1245 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));CHKERRQ(ierr); 1246 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));CHKERRQ(ierr); 1247 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));CHKERRQ(ierr); 1248 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));CHKERRQ(ierr); 1249 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));CHKERRQ(ierr); 1250 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); 1251 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); 1252 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); 1253 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); 1254 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));CHKERRQ(ierr); 1255 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); 1256 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); 1257 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));CHKERRQ(ierr); 1258 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));CHKERRQ(ierr); 1259 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));CHKERRQ(ierr); 1260 } 1261 } 1262 } 1263 PetscFunctionReturn(0); 1264 } 1265 1266 #undef __FUNCT__ 1267 #define __FUNCT__ "MatGetInfo_MUMPS" 1268 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1269 { 1270 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1271 1272 PetscFunctionBegin; 1273 info->block_size = 1.0; 1274 info->nz_allocated = mumps->id.INFOG(20); 1275 info->nz_used = mumps->id.INFOG(20); 1276 info->nz_unneeded = 0.0; 1277 info->assemblies = 0.0; 1278 info->mallocs = 0.0; 1279 info->memory = 0.0; 1280 info->fill_ratio_given = 0; 1281 info->fill_ratio_needed = 0; 1282 info->factor_mallocs = 0; 1283 PetscFunctionReturn(0); 1284 } 1285 1286 /* -------------------------------------------------------------------------------------------*/ 1287 #undef __FUNCT__ 1288 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1289 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1290 { 1291 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1292 1293 PetscFunctionBegin; 1294 mumps->id.ICNTL(icntl) = ival; 1295 PetscFunctionReturn(0); 1296 } 1297 1298 #undef __FUNCT__ 1299 #define __FUNCT__ "MatMumpsSetIcntl" 1300 /*@ 1301 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1302 1303 Logically Collective on Mat 1304 1305 Input Parameters: 1306 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1307 . icntl - index of MUMPS parameter array ICNTL() 1308 - ival - value of MUMPS ICNTL(icntl) 1309 1310 Options Database: 1311 . -mat_mumps_icntl_<icntl> <ival> 1312 1313 Level: beginner 1314 1315 References: MUMPS Users' Guide 1316 1317 .seealso: MatGetFactor() 1318 @*/ 1319 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1320 { 1321 PetscErrorCode ierr; 1322 1323 PetscFunctionBegin; 1324 PetscValidLogicalCollectiveInt(F,icntl,2); 1325 PetscValidLogicalCollectiveInt(F,ival,3); 1326 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1327 PetscFunctionReturn(0); 1328 } 1329 1330 /* -------------------------------------------------------------------------------------------*/ 1331 #undef __FUNCT__ 1332 #define __FUNCT__ "MatMumpsSetCntl_MUMPS" 1333 PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val) 1334 { 1335 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1336 1337 PetscFunctionBegin; 1338 mumps->id.CNTL(icntl) = val; 1339 PetscFunctionReturn(0); 1340 } 1341 1342 #undef __FUNCT__ 1343 #define __FUNCT__ "MatMumpsSetCntl" 1344 /*@ 1345 MatMumpsSetCntl - Set MUMPS parameter CNTL() 1346 1347 Logically Collective on Mat 1348 1349 Input Parameters: 1350 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1351 . icntl - index of MUMPS parameter array CNTL() 1352 - val - value of MUMPS CNTL(icntl) 1353 1354 Options Database: 1355 . -mat_mumps_cntl_<icntl> <val> 1356 1357 Level: beginner 1358 1359 References: MUMPS Users' Guide 1360 1361 .seealso: MatGetFactor() 1362 @*/ 1363 PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val) 1364 { 1365 PetscErrorCode ierr; 1366 1367 PetscFunctionBegin; 1368 PetscValidLogicalCollectiveInt(F,icntl,2); 1369 PetscValidLogicalCollectiveInt(F,val,3); 1370 ierr = PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));CHKERRQ(ierr); 1371 PetscFunctionReturn(0); 1372 } 1373 1374 /*MC 1375 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1376 distributed and sequential matrices via the external package MUMPS. 1377 1378 Works with MATAIJ and MATSBAIJ matrices 1379 1380 Options Database Keys: 1381 + -mat_mumps_icntl_4 <0,...,4> - print level 1382 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1383 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec) 1384 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1385 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1386 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1387 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1388 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1389 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1390 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1391 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1392 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1393 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1394 1395 Level: beginner 1396 1397 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1398 1399 M*/ 1400 1401 #undef __FUNCT__ 1402 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1403 static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1404 { 1405 PetscFunctionBegin; 1406 *type = MATSOLVERMUMPS; 1407 PetscFunctionReturn(0); 1408 } 1409 1410 /* MatGetFactor for Seq and MPI AIJ matrices */ 1411 #undef __FUNCT__ 1412 #define __FUNCT__ "MatGetFactor_aij_mumps" 1413 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1414 { 1415 Mat B; 1416 PetscErrorCode ierr; 1417 Mat_MUMPS *mumps; 1418 PetscBool isSeqAIJ; 1419 1420 PetscFunctionBegin; 1421 /* Create the factorization matrix */ 1422 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1423 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1424 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1425 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1426 if (isSeqAIJ) { 1427 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1428 } else { 1429 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1430 } 1431 1432 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1433 1434 B->ops->view = MatView_MUMPS; 1435 B->ops->getinfo = MatGetInfo_MUMPS; 1436 1437 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1438 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1439 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1440 if (ftype == MAT_FACTOR_LU) { 1441 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1442 B->factortype = MAT_FACTOR_LU; 1443 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1444 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1445 mumps->sym = 0; 1446 } else { 1447 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1448 B->factortype = MAT_FACTOR_CHOLESKY; 1449 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1450 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1451 if (A->spd_set && A->spd) mumps->sym = 1; 1452 else mumps->sym = 2; 1453 } 1454 1455 mumps->isAIJ = PETSC_TRUE; 1456 mumps->Destroy = B->ops->destroy; 1457 B->ops->destroy = MatDestroy_MUMPS; 1458 B->spptr = (void*)mumps; 1459 1460 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1461 1462 *F = B; 1463 PetscFunctionReturn(0); 1464 } 1465 1466 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1467 #undef __FUNCT__ 1468 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1469 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1470 { 1471 Mat B; 1472 PetscErrorCode ierr; 1473 Mat_MUMPS *mumps; 1474 PetscBool isSeqSBAIJ; 1475 1476 PetscFunctionBegin; 1477 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1478 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"); 1479 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1480 /* Create the factorization matrix */ 1481 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1482 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1483 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1484 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1485 if (isSeqSBAIJ) { 1486 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 1487 1488 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1489 } else { 1490 ierr = MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);CHKERRQ(ierr); 1491 1492 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1493 } 1494 1495 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1496 B->ops->view = MatView_MUMPS; 1497 1498 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1499 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl);CHKERRQ(ierr); 1500 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl);CHKERRQ(ierr); 1501 1502 B->factortype = MAT_FACTOR_CHOLESKY; 1503 if (A->spd_set && A->spd) mumps->sym = 1; 1504 else mumps->sym = 2; 1505 1506 mumps->isAIJ = PETSC_FALSE; 1507 mumps->Destroy = B->ops->destroy; 1508 B->ops->destroy = MatDestroy_MUMPS; 1509 B->spptr = (void*)mumps; 1510 1511 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1512 1513 *F = B; 1514 PetscFunctionReturn(0); 1515 } 1516 1517 #undef __FUNCT__ 1518 #define __FUNCT__ "MatGetFactor_baij_mumps" 1519 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 1520 { 1521 Mat B; 1522 PetscErrorCode ierr; 1523 Mat_MUMPS *mumps; 1524 PetscBool isSeqBAIJ; 1525 1526 PetscFunctionBegin; 1527 /* Create the factorization matrix */ 1528 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1529 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1530 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1531 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1532 if (isSeqBAIJ) { 1533 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);CHKERRQ(ierr); 1534 } else { 1535 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 1536 } 1537 1538 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1539 if (ftype == MAT_FACTOR_LU) { 1540 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1541 B->factortype = MAT_FACTOR_LU; 1542 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 1543 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 1544 mumps->sym = 0; 1545 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 1546 1547 B->ops->view = MatView_MUMPS; 1548 1549 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1550 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1551 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1552 1553 mumps->isAIJ = PETSC_TRUE; 1554 mumps->Destroy = B->ops->destroy; 1555 B->ops->destroy = MatDestroy_MUMPS; 1556 B->spptr = (void*)mumps; 1557 1558 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1559 1560 *F = B; 1561 PetscFunctionReturn(0); 1562 } 1563 1564