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