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