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