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