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