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