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