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