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; /* used by MatSolve() */ 77 Vec bb_mpi,bb_seq,xx_mpi; /* used by MatMatSolve() */ 78 PetscInt ICNTL9_pre; /* check if ICNTL(9) is changed from previous MatSolve */ 79 80 PetscErrorCode (*Destroy)(Mat); 81 PetscErrorCode (*ConvertToTriples)(Mat, int, MatReuse, int*, int**, int**, PetscScalar**); 82 } Mat_MUMPS; 83 84 extern PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*); 85 86 /* 87 MatConvertToTriples_A_B - 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_SeqAIJ *aa=(Mat_SeqAIJ*)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 /* count nz in the uppper triangular part of A */ 226 nz = 0; 227 for (i=0; i<M; i++) nz += ai[i+1] - adiag[i]; 228 *nnz = nz; 229 230 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 231 col = row + nz; 232 val = (PetscScalar*)(col + nz); 233 234 nz = 0; 235 for (i=0; i<M; i++) { 236 rnz = ai[i+1] - adiag[i]; 237 ajj = aj + adiag[i]; 238 v1 = av + adiag[i]; 239 for (j=0; j<rnz; j++) { 240 row[nz] = i+shift; col[nz] = ajj[j] + shift; val[nz++] = v1[j]; 241 } 242 } 243 *r = row; *c = col; *v = val; 244 } else { 245 nz = 0; val = *v; 246 for (i=0; i <M; i++) { 247 rnz = ai[i+1] - adiag[i]; 248 ajj = aj + adiag[i]; 249 v1 = av + adiag[i]; 250 for (j=0; j<rnz; j++) { 251 val[nz++] = v1[j]; 252 } 253 } 254 } 255 PetscFunctionReturn(0); 256 } 257 258 #undef __FUNCT__ 259 #define __FUNCT__ "MatConvertToTriples_mpisbaij_mpisbaij" 260 PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 261 { 262 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 263 PetscErrorCode ierr; 264 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 265 PetscInt *row,*col; 266 const PetscScalar *av, *bv,*v1,*v2; 267 PetscScalar *val; 268 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 269 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)(mat->A)->data; 270 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 271 272 PetscFunctionBegin; 273 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 274 av=aa->a; bv=bb->a; 275 276 garray = mat->garray; 277 278 if (reuse == MAT_INITIAL_MATRIX) { 279 nz = aa->nz + bb->nz; 280 *nnz = nz; 281 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 282 col = row + nz; 283 val = (PetscScalar*)(col + nz); 284 285 *r = row; *c = col; *v = val; 286 } else { 287 row = *r; col = *c; val = *v; 288 } 289 290 jj = 0; irow = rstart; 291 for (i=0; i<m; i++) { 292 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 293 countA = ai[i+1] - ai[i]; 294 countB = bi[i+1] - bi[i]; 295 bjj = bj + bi[i]; 296 v1 = av + ai[i]; 297 v2 = bv + bi[i]; 298 299 /* A-part */ 300 for (j=0; j<countA; j++) { 301 if (reuse == MAT_INITIAL_MATRIX) { 302 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 303 } 304 val[jj++] = v1[j]; 305 } 306 307 /* B-part */ 308 for (j=0; j < countB; j++) { 309 if (reuse == MAT_INITIAL_MATRIX) { 310 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 311 } 312 val[jj++] = v2[j]; 313 } 314 irow++; 315 } 316 PetscFunctionReturn(0); 317 } 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpiaij" 321 PetscErrorCode MatConvertToTriples_mpiaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 322 { 323 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 324 PetscErrorCode ierr; 325 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 326 PetscInt *row,*col; 327 const PetscScalar *av, *bv,*v1,*v2; 328 PetscScalar *val; 329 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 330 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(mat->A)->data; 331 Mat_SeqAIJ *bb = (Mat_SeqAIJ*)(mat->B)->data; 332 333 PetscFunctionBegin; 334 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 335 av=aa->a; bv=bb->a; 336 337 garray = mat->garray; 338 339 if (reuse == MAT_INITIAL_MATRIX) { 340 nz = aa->nz + bb->nz; 341 *nnz = nz; 342 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 343 col = row + nz; 344 val = (PetscScalar*)(col + nz); 345 346 *r = row; *c = col; *v = val; 347 } else { 348 row = *r; col = *c; val = *v; 349 } 350 351 jj = 0; irow = rstart; 352 for (i=0; i<m; i++) { 353 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 354 countA = ai[i+1] - ai[i]; 355 countB = bi[i+1] - bi[i]; 356 bjj = bj + bi[i]; 357 v1 = av + ai[i]; 358 v2 = bv + bi[i]; 359 360 /* A-part */ 361 for (j=0; j<countA; j++) { 362 if (reuse == MAT_INITIAL_MATRIX) { 363 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 364 } 365 val[jj++] = v1[j]; 366 } 367 368 /* B-part */ 369 for (j=0; j < countB; j++) { 370 if (reuse == MAT_INITIAL_MATRIX) { 371 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 372 } 373 val[jj++] = v2[j]; 374 } 375 irow++; 376 } 377 PetscFunctionReturn(0); 378 } 379 380 #undef __FUNCT__ 381 #define __FUNCT__ "MatConvertToTriples_mpibaij_mpiaij" 382 PetscErrorCode MatConvertToTriples_mpibaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 383 { 384 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data; 385 Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)(mat->A)->data; 386 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 387 const PetscInt *ai = aa->i, *bi = bb->i, *aj = aa->j, *bj = bb->j,*ajj, *bjj; 388 const PetscInt *garray = mat->garray,mbs=mat->mbs,rstart=A->rmap->rstart; 389 const PetscInt bs2=mat->bs2; 390 PetscErrorCode ierr; 391 PetscInt bs,nz,i,j,k,n,jj,irow,countA,countB,idx; 392 PetscInt *row,*col; 393 const PetscScalar *av=aa->a, *bv=bb->a,*v1,*v2; 394 PetscScalar *val; 395 396 PetscFunctionBegin; 397 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 398 if (reuse == MAT_INITIAL_MATRIX) { 399 nz = bs2*(aa->nz + bb->nz); 400 *nnz = nz; 401 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 402 col = row + nz; 403 val = (PetscScalar*)(col + nz); 404 405 *r = row; *c = col; *v = val; 406 } else { 407 row = *r; col = *c; val = *v; 408 } 409 410 jj = 0; irow = rstart; 411 for (i=0; i<mbs; i++) { 412 countA = ai[i+1] - ai[i]; 413 countB = bi[i+1] - bi[i]; 414 ajj = aj + ai[i]; 415 bjj = bj + bi[i]; 416 v1 = av + bs2*ai[i]; 417 v2 = bv + bs2*bi[i]; 418 419 idx = 0; 420 /* A-part */ 421 for (k=0; k<countA; k++) { 422 for (j=0; j<bs; j++) { 423 for (n=0; n<bs; n++) { 424 if (reuse == MAT_INITIAL_MATRIX) { 425 row[jj] = irow + n + shift; 426 col[jj] = rstart + bs*ajj[k] + j + shift; 427 } 428 val[jj++] = v1[idx++]; 429 } 430 } 431 } 432 433 idx = 0; 434 /* B-part */ 435 for (k=0; k<countB; k++) { 436 for (j=0; j<bs; j++) { 437 for (n=0; n<bs; n++) { 438 if (reuse == MAT_INITIAL_MATRIX) { 439 row[jj] = irow + n + shift; 440 col[jj] = bs*garray[bjj[k]] + j + shift; 441 } 442 val[jj++] = v2[idx++]; 443 } 444 } 445 } 446 irow += bs; 447 } 448 PetscFunctionReturn(0); 449 } 450 451 #undef __FUNCT__ 452 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpisbaij" 453 PetscErrorCode MatConvertToTriples_mpiaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 454 { 455 const PetscInt *ai, *aj,*adiag, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 456 PetscErrorCode ierr; 457 PetscInt rstart,nz,nza,nzb,i,j,jj,irow,countA,countB; 458 PetscInt *row,*col; 459 const PetscScalar *av, *bv,*v1,*v2; 460 PetscScalar *val; 461 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 462 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(mat->A)->data; 463 Mat_SeqAIJ *bb =(Mat_SeqAIJ*)(mat->B)->data; 464 465 PetscFunctionBegin; 466 ai=aa->i; aj=aa->j; adiag=aa->diag; 467 bi=bb->i; bj=bb->j; garray = mat->garray; 468 av=aa->a; bv=bb->a; 469 470 rstart = A->rmap->rstart; 471 472 if (reuse == MAT_INITIAL_MATRIX) { 473 nza = 0; /* num of upper triangular entries in mat->A, including diagonals */ 474 nzb = 0; /* num of upper triangular entries in mat->B */ 475 for (i=0; i<m; i++) { 476 nza += (ai[i+1] - adiag[i]); 477 countB = bi[i+1] - bi[i]; 478 bjj = bj + bi[i]; 479 for (j=0; j<countB; j++) { 480 if (garray[bjj[j]] > rstart) nzb++; 481 } 482 } 483 484 nz = nza + nzb; /* total nz of upper triangular part of mat */ 485 *nnz = nz; 486 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 487 col = row + nz; 488 val = (PetscScalar*)(col + nz); 489 490 *r = row; *c = col; *v = val; 491 } else { 492 row = *r; col = *c; val = *v; 493 } 494 495 jj = 0; irow = rstart; 496 for (i=0; i<m; i++) { 497 ajj = aj + adiag[i]; /* ptr to the beginning of the diagonal of this row */ 498 v1 = av + adiag[i]; 499 countA = ai[i+1] - adiag[i]; 500 countB = bi[i+1] - bi[i]; 501 bjj = bj + bi[i]; 502 v2 = bv + bi[i]; 503 504 /* A-part */ 505 for (j=0; j<countA; j++) { 506 if (reuse == MAT_INITIAL_MATRIX) { 507 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 508 } 509 val[jj++] = v1[j]; 510 } 511 512 /* B-part */ 513 for (j=0; j < countB; j++) { 514 if (garray[bjj[j]] > rstart) { 515 if (reuse == MAT_INITIAL_MATRIX) { 516 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 517 } 518 val[jj++] = v2[j]; 519 } 520 } 521 irow++; 522 } 523 PetscFunctionReturn(0); 524 } 525 526 #undef __FUNCT__ 527 #define __FUNCT__ "MatGetDiagonal_MUMPS" 528 PetscErrorCode MatGetDiagonal_MUMPS(Mat A,Vec v) 529 { 530 PetscFunctionBegin; 531 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: MUMPS factor"); 532 PetscFunctionReturn(0); 533 } 534 535 #undef __FUNCT__ 536 #define __FUNCT__ "MatDestroy_MUMPS" 537 PetscErrorCode MatDestroy_MUMPS(Mat A) 538 { 539 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 540 PetscErrorCode ierr; 541 542 PetscFunctionBegin; 543 if (mumps->CleanUpMUMPS) { 544 /* Terminate instance, deallocate memories */ 545 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 546 ierr = VecScatterDestroy(&mumps->scat_rhs);CHKERRQ(ierr); 547 ierr = VecDestroy(&mumps->b_seq);CHKERRQ(ierr); 548 ierr = VecDestroy(&mumps->bb_seq);CHKERRQ(ierr); 549 ierr = VecDestroy(&mumps->bb_mpi);CHKERRQ(ierr); 550 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 551 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 552 ierr = PetscFree(mumps->id.perm_in);CHKERRQ(ierr); 553 ierr = PetscFree(mumps->irn);CHKERRQ(ierr); 554 555 mumps->id.job = JOB_END; 556 PetscMUMPS_c(&mumps->id); 557 ierr = MPI_Comm_free(&(mumps->comm_mumps));CHKERRQ(ierr); 558 } 559 if (mumps->Destroy) { 560 ierr = (mumps->Destroy)(A);CHKERRQ(ierr); 561 } 562 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 563 564 /* clear composed functions */ 565 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);CHKERRQ(ierr); 566 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetIcntl_C",NULL);CHKERRQ(ierr); 567 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetIcntl_C",NULL);CHKERRQ(ierr); 568 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetCntl_C",NULL);CHKERRQ(ierr); 569 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetCntl_C",NULL);CHKERRQ(ierr); 570 571 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfo_C",NULL);CHKERRQ(ierr); 572 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfog_C",NULL);CHKERRQ(ierr); 573 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfo_C",NULL);CHKERRQ(ierr); 574 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfog_C",NULL);CHKERRQ(ierr); 575 PetscFunctionReturn(0); 576 } 577 578 #undef __FUNCT__ 579 #define __FUNCT__ "MatSolve_MUMPS" 580 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 581 { 582 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 583 PetscScalar *array; 584 Vec b_seq; 585 IS is_iden,is_petsc; 586 PetscErrorCode ierr; 587 PetscInt i; 588 static PetscBool cite1 = PETSC_FALSE,cite2 = PETSC_FALSE; 589 590 PetscFunctionBegin; 591 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); 592 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); 593 mumps->id.nrhs = 1; 594 b_seq = mumps->b_seq; 595 if (mumps->size > 1) { 596 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 597 ierr = VecScatterBegin(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 598 ierr = VecScatterEnd(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 599 if (!mumps->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 600 } else { /* size == 1 */ 601 ierr = VecCopy(b,x);CHKERRQ(ierr); 602 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 603 } 604 if (!mumps->myid) { /* define rhs on the host */ 605 mumps->id.nrhs = 1; 606 #if defined(PETSC_USE_COMPLEX) 607 #if defined(PETSC_USE_REAL_SINGLE) 608 mumps->id.rhs = (mumps_complex*)array; 609 #else 610 mumps->id.rhs = (mumps_double_complex*)array; 611 #endif 612 #else 613 mumps->id.rhs = array; 614 #endif 615 } 616 617 /* solve phase */ 618 /*-------------*/ 619 mumps->id.job = JOB_SOLVE; 620 PetscMUMPS_c(&mumps->id); 621 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)); 622 623 if (mumps->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 624 if (mumps->scat_sol && mumps->ICNTL9_pre != mumps->id.ICNTL(9)) { 625 /* when id.ICNTL(9) changes, the contents of lsol_loc may change (not its size, lsol_loc), recreates scat_sol */ 626 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 627 } 628 if (!mumps->scat_sol) { /* create scatter scat_sol */ 629 ierr = ISCreateStride(PETSC_COMM_SELF,mumps->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 630 for (i=0; i<mumps->id.lsol_loc; i++) { 631 mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 632 } 633 ierr = ISCreateGeneral(PETSC_COMM_SELF,mumps->id.lsol_loc,mumps->id.isol_loc,PETSC_COPY_VALUES,&is_petsc);CHKERRQ(ierr); /* to */ 634 ierr = VecScatterCreate(mumps->x_seq,is_iden,x,is_petsc,&mumps->scat_sol);CHKERRQ(ierr); 635 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 636 ierr = ISDestroy(&is_petsc);CHKERRQ(ierr); 637 638 mumps->ICNTL9_pre = mumps->id.ICNTL(9); /* save current value of id.ICNTL(9) */ 639 } 640 641 ierr = VecScatterBegin(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 642 ierr = VecScatterEnd(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 643 } 644 PetscFunctionReturn(0); 645 } 646 647 #undef __FUNCT__ 648 #define __FUNCT__ "MatSolveTranspose_MUMPS" 649 PetscErrorCode MatSolveTranspose_MUMPS(Mat A,Vec b,Vec x) 650 { 651 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 652 PetscErrorCode ierr; 653 654 PetscFunctionBegin; 655 mumps->id.ICNTL(9) = 0; 656 657 ierr = MatSolve_MUMPS(A,b,x);CHKERRQ(ierr); 658 659 mumps->id.ICNTL(9) = 1; 660 PetscFunctionReturn(0); 661 } 662 663 #undef __FUNCT__ 664 #define __FUNCT__ "MatMatSolve_MUMPS" 665 PetscErrorCode MatMatSolve_MUMPS(Mat A,Mat B,Mat X) 666 { 667 PetscErrorCode ierr; 668 PetscBool flg; 669 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 670 PetscInt i,nrhs,m,M; 671 PetscScalar *array,*bray; 672 673 PetscFunctionBegin; 674 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 675 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 676 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 677 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 678 679 ierr = MatGetLocalSize(B,&m,NULL);CHKERRQ(ierr); 680 ierr = MatGetSize(B,&M,&nrhs);CHKERRQ(ierr); 681 682 if (mumps->size == 1) { 683 /* copy B to X */ 684 ierr = MatDenseGetArray(B,&bray);CHKERRQ(ierr); 685 ierr = MatDenseGetArray(X,&array);CHKERRQ(ierr); 686 for (i=0; i<M*nrhs; i++) array[i] = bray[i]; 687 ierr = MatDenseRestoreArray(B,&bray);CHKERRQ(ierr); 688 689 //if (!mumps->myid) { /* define rhs on the host */ 690 mumps->id.nrhs = nrhs; 691 mumps->id.lrhs = M; 692 #if defined(PETSC_USE_COMPLEX) 693 #if defined(PETSC_USE_REAL_SINGLE) 694 mumps->id.rhs = (mumps_complex*)array; 695 #else 696 mumps->id.rhs = (mumps_double_complex*)array; 697 #endif 698 #else 699 mumps->id.rhs = array; 700 #endif 701 //} 702 /* solve phase */ 703 /*-------------*/ 704 mumps->id.job = JOB_SOLVE; 705 PetscMUMPS_c(&mumps->id); 706 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)); 707 708 ierr = MatDenseRestoreArray(X,&array);CHKERRQ(ierr); 709 PetscFunctionReturn(0); 710 } 711 712 /* copy rhs matrix B into vector bb_mpi */ 713 if (mumps->bb_mpi) { 714 ierr = VecDestroy(&mumps->bb_mpi);CHKERRQ(ierr); 715 ierr = VecDestroy(&mumps->bb_seq);CHKERRQ(ierr); 716 } 717 ierr = MatDenseGetArray(B,&array);CHKERRQ(ierr); 718 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)B),1,nrhs*m,nrhs*M,(const PetscScalar*)array,&mumps->bb_mpi);CHKERRQ(ierr); 719 //ierr = VecView(mumps->bb_mpi,PETSC_VIEWER_STDOUT_WORLD); 720 ierr = MatDenseRestoreArray(B,&array);CHKERRQ(ierr); 721 722 /* scatter bb_mpi to bb_seq because MUMPS only supports centralized rhs */ 723 if (!mumps->myid) { 724 ierr = VecCreateSeq(PETSC_COMM_SELF,nrhs*M,&mumps->bb_seq);CHKERRQ(ierr); 725 //ierr = ISCreateStride(PETSC_COMM_SELF,A->rmap->N,0,1,&is_iden);CHKERRQ(ierr); 726 } else { 727 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->bb_seq);CHKERRQ(ierr); 728 //ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 729 } 730 731 //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_MUMPS() is not implemented yet"); 732 PetscFunctionReturn(0); 733 } 734 735 #if !defined(PETSC_USE_COMPLEX) 736 /* 737 input: 738 F: numeric factor 739 output: 740 nneg: total number of negative pivots 741 nzero: 0 742 npos: (global dimension of F) - nneg 743 */ 744 745 #undef __FUNCT__ 746 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 747 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 748 { 749 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 750 PetscErrorCode ierr; 751 PetscMPIInt size; 752 753 PetscFunctionBegin; 754 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);CHKERRQ(ierr); 755 /* 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 */ 756 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)); 757 758 if (nneg) *nneg = mumps->id.INFOG(12); 759 if (nzero || npos) { 760 if (mumps->id.ICNTL(24) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"-mat_mumps_icntl_24 must be set as 1 for null pivot row detection"); 761 if (nzero) *nzero = mumps->id.INFOG(28); 762 if (npos) *npos = F->rmap->N - (mumps->id.INFOG(12) + mumps->id.INFOG(28)); 763 } 764 PetscFunctionReturn(0); 765 } 766 #endif /* !defined(PETSC_USE_COMPLEX) */ 767 768 #undef __FUNCT__ 769 #define __FUNCT__ "MatFactorNumeric_MUMPS" 770 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 771 { 772 Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr; 773 PetscErrorCode ierr; 774 Mat F_diag; 775 PetscBool isMPIAIJ; 776 777 PetscFunctionBegin; 778 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 779 780 /* numerical factorization phase */ 781 /*-------------------------------*/ 782 mumps->id.job = JOB_FACTNUMERIC; 783 if (!mumps->id.ICNTL(18)) { /* A is centralized */ 784 if (!mumps->myid) { 785 #if defined(PETSC_USE_COMPLEX) 786 #if defined(PETSC_USE_REAL_SINGLE) 787 mumps->id.a = (mumps_complex*)mumps->val; 788 #else 789 mumps->id.a = (mumps_double_complex*)mumps->val; 790 #endif 791 #else 792 mumps->id.a = mumps->val; 793 #endif 794 } 795 } else { 796 #if defined(PETSC_USE_COMPLEX) 797 #if defined(PETSC_USE_REAL_SINGLE) 798 mumps->id.a_loc = (mumps_complex*)mumps->val; 799 #else 800 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 801 #endif 802 #else 803 mumps->id.a_loc = mumps->val; 804 #endif 805 } 806 PetscMUMPS_c(&mumps->id); 807 if (mumps->id.INFOG(1) < 0) { 808 if (mumps->id.INFO(1) == -13) { 809 if (mumps->id.INFO(2) < 0) { 810 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)); 811 } else { 812 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d bytes\n",mumps->id.INFO(2)); 813 } 814 } 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)); 815 } 816 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)); 817 818 (F)->assembled = PETSC_TRUE; 819 mumps->matstruc = SAME_NONZERO_PATTERN; 820 mumps->CleanUpMUMPS = PETSC_TRUE; 821 822 if (mumps->size > 1) { 823 PetscInt lsol_loc; 824 PetscScalar *sol_loc; 825 826 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 827 if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A; 828 else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A; 829 F_diag->assembled = PETSC_TRUE; 830 831 /* distributed solution; Create x_seq=sol_loc for repeated use */ 832 if (mumps->x_seq) { 833 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 834 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 835 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 836 } 837 lsol_loc = mumps->id.INFO(23); /* length of sol_loc */ 838 ierr = PetscMalloc2(lsol_loc,&sol_loc,lsol_loc,&mumps->id.isol_loc);CHKERRQ(ierr); 839 mumps->id.lsol_loc = lsol_loc; 840 #if defined(PETSC_USE_COMPLEX) 841 #if defined(PETSC_USE_REAL_SINGLE) 842 mumps->id.sol_loc = (mumps_complex*)sol_loc; 843 #else 844 mumps->id.sol_loc = (mumps_double_complex*)sol_loc; 845 #endif 846 #else 847 mumps->id.sol_loc = sol_loc; 848 #endif 849 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 850 } 851 PetscFunctionReturn(0); 852 } 853 854 /* Sets MUMPS options from the options database */ 855 #undef __FUNCT__ 856 #define __FUNCT__ "PetscSetMUMPSFromOptions" 857 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 858 { 859 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 860 PetscErrorCode ierr; 861 PetscInt icntl; 862 PetscBool flg; 863 864 PetscFunctionBegin; 865 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 866 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 867 if (flg) mumps->id.ICNTL(1) = icntl; 868 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); 869 if (flg) mumps->id.ICNTL(2) = icntl; 870 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); 871 if (flg) mumps->id.ICNTL(3) = icntl; 872 873 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 874 if (flg) mumps->id.ICNTL(4) = icntl; 875 if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 876 877 ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): permutes to a zero-free diagonal and/or scale the matrix (0 to 7)","None",mumps->id.ICNTL(6),&icntl,&flg);CHKERRQ(ierr); 878 if (flg) mumps->id.ICNTL(6) = icntl; 879 880 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): computes a symmetric permutation in sequential analysis (0 to 7). 3=Scotch, 4=PORD, 5=Metis","None",mumps->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 881 if (flg) { 882 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"); 883 else mumps->id.ICNTL(7) = icntl; 884 } 885 886 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); 887 /* ierr = PetscOptionsInt("-mat_mumps_icntl_9","ICNTL(9): computes the solution using A or A^T","None",mumps->id.ICNTL(9),&mumps->id.ICNTL(9),NULL);CHKERRQ(ierr); handled by MatSolveTranspose_MUMPS() */ 888 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);CHKERRQ(ierr); 889 ierr = PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): statistics related to an error analysis (via -ksp_view)","None",mumps->id.ICNTL(11),&mumps->id.ICNTL(11),NULL);CHKERRQ(ierr); 890 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): an ordering strategy for symmetric matrices (0 to 3)","None",mumps->id.ICNTL(12),&mumps->id.ICNTL(12),NULL);CHKERRQ(ierr); 891 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): parallelism of the root node (enable ScaLAPACK) and its splitting","None",mumps->id.ICNTL(13),&mumps->id.ICNTL(13),NULL);CHKERRQ(ierr); 892 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): percentage increase in the estimated working space","None",mumps->id.ICNTL(14),&mumps->id.ICNTL(14),NULL);CHKERRQ(ierr); 893 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): computes the Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);CHKERRQ(ierr); 894 /* ierr = PetscOptionsInt("-mat_mumps_icntl_20","ICNTL(20): the format (dense or sparse) of the right-hand sides","None",mumps->id.ICNTL(20),&mumps->id.ICNTL(20),NULL);CHKERRQ(ierr); -- sparse rhs is not supported in PETSc API */ 895 /* ierr = PetscOptionsInt("-mat_mumps_icntl_21","ICNTL(21): the distribution (centralized or distributed) of the solution vectors","None",mumps->id.ICNTL(21),&mumps->id.ICNTL(21),NULL);CHKERRQ(ierr); we only use distributed solution vector */ 896 897 ierr = PetscOptionsInt("-mat_mumps_icntl_22","ICNTL(22): in-core/out-of-core factorization and solve (0 or 1)","None",mumps->id.ICNTL(22),&mumps->id.ICNTL(22),NULL);CHKERRQ(ierr); 898 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); 899 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); 900 if (mumps->id.ICNTL(24)) { 901 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 902 } 903 904 ierr = PetscOptionsInt("-mat_mumps_icntl_25","ICNTL(25): compute a solution of a deficient matrix and a null space basis","None",mumps->id.ICNTL(25),&mumps->id.ICNTL(25),NULL);CHKERRQ(ierr); 905 ierr = PetscOptionsInt("-mat_mumps_icntl_26","ICNTL(26): drives the solution phase if a Schur complement matrix","None",mumps->id.ICNTL(26),&mumps->id.ICNTL(26),NULL);CHKERRQ(ierr); 906 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): the blocking size for multiple right-hand sides","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);CHKERRQ(ierr); 907 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); 908 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); 909 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); 910 ierr = PetscOptionsInt("-mat_mumps_icntl_31","ICNTL(31): indicates which factors may be discarded during factorization","None",mumps->id.ICNTL(31),&mumps->id.ICNTL(31),NULL);CHKERRQ(ierr); 911 /* ierr = PetscOptionsInt("-mat_mumps_icntl_32","ICNTL(32): performs the forward elemination of the right-hand sides during factorization","None",mumps->id.ICNTL(32),&mumps->id.ICNTL(32),NULL);CHKERRQ(ierr); -- not supported by PETSc API */ 912 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);CHKERRQ(ierr); 913 914 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);CHKERRQ(ierr); 915 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);CHKERRQ(ierr); 916 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);CHKERRQ(ierr); 917 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);CHKERRQ(ierr); 918 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);CHKERRQ(ierr); 919 920 ierr = PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL); 921 PetscOptionsEnd(); 922 PetscFunctionReturn(0); 923 } 924 925 #undef __FUNCT__ 926 #define __FUNCT__ "PetscInitializeMUMPS" 927 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps) 928 { 929 PetscErrorCode ierr; 930 931 PetscFunctionBegin; 932 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid); 933 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);CHKERRQ(ierr); 934 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));CHKERRQ(ierr); 935 936 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 937 938 mumps->id.job = JOB_INIT; 939 mumps->id.par = 1; /* host participates factorizaton and solve */ 940 mumps->id.sym = mumps->sym; 941 PetscMUMPS_c(&mumps->id); 942 943 mumps->CleanUpMUMPS = PETSC_FALSE; 944 mumps->scat_rhs = NULL; 945 mumps->scat_sol = NULL; 946 947 /* set PETSc-MUMPS default options - override MUMPS default */ 948 mumps->id.ICNTL(3) = 0; 949 mumps->id.ICNTL(4) = 0; 950 if (mumps->size == 1) { 951 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 952 } else { 953 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 954 mumps->id.ICNTL(20) = 0; /* rhs is in dense format */ 955 mumps->id.ICNTL(21) = 1; /* distributed solution */ 956 } 957 PetscFunctionReturn(0); 958 } 959 960 /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */ 961 #undef __FUNCT__ 962 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 963 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 964 { 965 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 966 PetscErrorCode ierr; 967 Vec b; 968 IS is_iden; 969 const PetscInt M = A->rmap->N; 970 971 PetscFunctionBegin; 972 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 973 974 /* Set MUMPS options from the options database */ 975 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 976 977 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 978 979 /* analysis phase */ 980 /*----------------*/ 981 mumps->id.job = JOB_FACTSYMBOLIC; 982 mumps->id.n = M; 983 switch (mumps->id.ICNTL(18)) { 984 case 0: /* centralized assembled matrix input */ 985 if (!mumps->myid) { 986 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 987 if (mumps->id.ICNTL(6)>1) { 988 #if defined(PETSC_USE_COMPLEX) 989 #if defined(PETSC_USE_REAL_SINGLE) 990 mumps->id.a = (mumps_complex*)mumps->val; 991 #else 992 mumps->id.a = (mumps_double_complex*)mumps->val; 993 #endif 994 #else 995 mumps->id.a = mumps->val; 996 #endif 997 } 998 if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */ 999 /* 1000 PetscBool flag; 1001 ierr = ISEqual(r,c,&flag);CHKERRQ(ierr); 1002 if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm"); 1003 ierr = ISView(r,PETSC_VIEWER_STDOUT_SELF); 1004 */ 1005 if (!mumps->myid) { 1006 const PetscInt *idx; 1007 PetscInt i,*perm_in; 1008 1009 ierr = PetscMalloc1(M,&perm_in);CHKERRQ(ierr); 1010 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 1011 1012 mumps->id.perm_in = perm_in; 1013 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 1014 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 1015 } 1016 } 1017 } 1018 break; 1019 case 3: /* distributed assembled matrix input (size>1) */ 1020 mumps->id.nz_loc = mumps->nz; 1021 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1022 if (mumps->id.ICNTL(6)>1) { 1023 #if defined(PETSC_USE_COMPLEX) 1024 #if defined(PETSC_USE_REAL_SINGLE) 1025 mumps->id.a_loc = (mumps_complex*)mumps->val; 1026 #else 1027 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1028 #endif 1029 #else 1030 mumps->id.a_loc = mumps->val; 1031 #endif 1032 } 1033 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1034 if (!mumps->myid) { 1035 ierr = VecCreateSeq(PETSC_COMM_SELF,A->rmap->N,&mumps->b_seq);CHKERRQ(ierr); 1036 ierr = ISCreateStride(PETSC_COMM_SELF,A->rmap->N,0,1,&is_iden);CHKERRQ(ierr); 1037 } else { 1038 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1039 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1040 } 1041 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1042 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1043 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1044 ierr = VecDestroy(&b);CHKERRQ(ierr); 1045 break; 1046 } 1047 PetscMUMPS_c(&mumps->id); 1048 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)); 1049 1050 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1051 F->ops->solve = MatSolve_MUMPS; 1052 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1053 F->ops->matsolve = MatMatSolve_MUMPS; 1054 PetscFunctionReturn(0); 1055 } 1056 1057 /* Note the Petsc r and c permutations are ignored */ 1058 #undef __FUNCT__ 1059 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 1060 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1061 { 1062 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1063 PetscErrorCode ierr; 1064 Vec b; 1065 IS is_iden; 1066 const PetscInt M = A->rmap->N; 1067 1068 PetscFunctionBegin; 1069 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1070 1071 /* Set MUMPS options from the options database */ 1072 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1073 1074 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1075 1076 /* analysis phase */ 1077 /*----------------*/ 1078 mumps->id.job = JOB_FACTSYMBOLIC; 1079 mumps->id.n = M; 1080 switch (mumps->id.ICNTL(18)) { 1081 case 0: /* centralized assembled matrix input */ 1082 if (!mumps->myid) { 1083 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1084 if (mumps->id.ICNTL(6)>1) { 1085 #if defined(PETSC_USE_COMPLEX) 1086 #if defined(PETSC_USE_REAL_SINGLE) 1087 mumps->id.a = (mumps_complex*)mumps->val; 1088 #else 1089 mumps->id.a = (mumps_double_complex*)mumps->val; 1090 #endif 1091 #else 1092 mumps->id.a = mumps->val; 1093 #endif 1094 } 1095 } 1096 break; 1097 case 3: /* distributed assembled matrix input (size>1) */ 1098 mumps->id.nz_loc = mumps->nz; 1099 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1100 if (mumps->id.ICNTL(6)>1) { 1101 #if defined(PETSC_USE_COMPLEX) 1102 #if defined(PETSC_USE_REAL_SINGLE) 1103 mumps->id.a_loc = (mumps_complex*)mumps->val; 1104 #else 1105 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1106 #endif 1107 #else 1108 mumps->id.a_loc = mumps->val; 1109 #endif 1110 } 1111 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1112 if (!mumps->myid) { 1113 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1114 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1115 } else { 1116 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1117 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1118 } 1119 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1120 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1121 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1122 ierr = VecDestroy(&b);CHKERRQ(ierr); 1123 break; 1124 } 1125 PetscMUMPS_c(&mumps->id); 1126 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)); 1127 1128 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1129 F->ops->solve = MatSolve_MUMPS; 1130 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1131 PetscFunctionReturn(0); 1132 } 1133 1134 /* Note the Petsc r permutation and factor info are ignored */ 1135 #undef __FUNCT__ 1136 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 1137 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 1138 { 1139 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1140 PetscErrorCode ierr; 1141 Vec b; 1142 IS is_iden; 1143 const PetscInt M = A->rmap->N; 1144 1145 PetscFunctionBegin; 1146 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1147 1148 /* Set MUMPS options from the options database */ 1149 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1150 1151 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1152 1153 /* analysis phase */ 1154 /*----------------*/ 1155 mumps->id.job = JOB_FACTSYMBOLIC; 1156 mumps->id.n = M; 1157 switch (mumps->id.ICNTL(18)) { 1158 case 0: /* centralized assembled matrix input */ 1159 if (!mumps->myid) { 1160 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1161 if (mumps->id.ICNTL(6)>1) { 1162 #if defined(PETSC_USE_COMPLEX) 1163 #if defined(PETSC_USE_REAL_SINGLE) 1164 mumps->id.a = (mumps_complex*)mumps->val; 1165 #else 1166 mumps->id.a = (mumps_double_complex*)mumps->val; 1167 #endif 1168 #else 1169 mumps->id.a = mumps->val; 1170 #endif 1171 } 1172 } 1173 break; 1174 case 3: /* distributed assembled matrix input (size>1) */ 1175 mumps->id.nz_loc = mumps->nz; 1176 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1177 if (mumps->id.ICNTL(6)>1) { 1178 #if defined(PETSC_USE_COMPLEX) 1179 #if defined(PETSC_USE_REAL_SINGLE) 1180 mumps->id.a_loc = (mumps_complex*)mumps->val; 1181 #else 1182 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1183 #endif 1184 #else 1185 mumps->id.a_loc = mumps->val; 1186 #endif 1187 } 1188 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1189 if (!mumps->myid) { 1190 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1191 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1192 } else { 1193 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1194 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1195 } 1196 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1197 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1198 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1199 ierr = VecDestroy(&b);CHKERRQ(ierr); 1200 break; 1201 } 1202 PetscMUMPS_c(&mumps->id); 1203 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)); 1204 1205 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1206 F->ops->solve = MatSolve_MUMPS; 1207 F->ops->solvetranspose = MatSolve_MUMPS; 1208 F->ops->matsolve = MatMatSolve_MUMPS; 1209 #if defined(PETSC_USE_COMPLEX) 1210 F->ops->getinertia = NULL; 1211 #else 1212 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1213 #endif 1214 PetscFunctionReturn(0); 1215 } 1216 1217 //update!!! 1218 #undef __FUNCT__ 1219 #define __FUNCT__ "MatView_MUMPS" 1220 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1221 { 1222 PetscErrorCode ierr; 1223 PetscBool iascii; 1224 PetscViewerFormat format; 1225 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1226 1227 PetscFunctionBegin; 1228 /* check if matrix is mumps type */ 1229 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1230 1231 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1232 if (iascii) { 1233 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1234 if (format == PETSC_VIEWER_ASCII_INFO) { 1235 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1236 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1237 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1238 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));CHKERRQ(ierr); 1239 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));CHKERRQ(ierr); 1240 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));CHKERRQ(ierr); 1241 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));CHKERRQ(ierr); 1242 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));CHKERRQ(ierr); 1243 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));CHKERRQ(ierr); 1244 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));CHKERRQ(ierr); 1245 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));CHKERRQ(ierr); 1246 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));CHKERRQ(ierr); 1247 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));CHKERRQ(ierr); 1248 if (mumps->id.ICNTL(11)>0) { 1249 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));CHKERRQ(ierr); 1250 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));CHKERRQ(ierr); 1251 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));CHKERRQ(ierr); 1252 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));CHKERRQ(ierr); 1253 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));CHKERRQ(ierr); 1254 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));CHKERRQ(ierr); 1255 } 1256 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));CHKERRQ(ierr); 1257 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));CHKERRQ(ierr); 1258 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));CHKERRQ(ierr); 1259 /* ICNTL(15-17) not used */ 1260 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));CHKERRQ(ierr); 1261 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));CHKERRQ(ierr); 1262 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));CHKERRQ(ierr); 1263 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));CHKERRQ(ierr); 1264 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));CHKERRQ(ierr); 1265 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));CHKERRQ(ierr); 1266 1267 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));CHKERRQ(ierr); 1268 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));CHKERRQ(ierr); 1269 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));CHKERRQ(ierr); 1270 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));CHKERRQ(ierr); 1271 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));CHKERRQ(ierr); 1272 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));CHKERRQ(ierr); 1273 1274 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));CHKERRQ(ierr); 1275 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));CHKERRQ(ierr); 1276 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));CHKERRQ(ierr); 1277 1278 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));CHKERRQ(ierr); 1279 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));CHKERRQ(ierr); 1280 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));CHKERRQ(ierr); 1281 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));CHKERRQ(ierr); 1282 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));CHKERRQ(ierr); 1283 1284 /* infomation local to each processor */ 1285 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1286 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1287 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));CHKERRQ(ierr); 1288 ierr = PetscViewerFlush(viewer); 1289 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1290 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));CHKERRQ(ierr); 1291 ierr = PetscViewerFlush(viewer); 1292 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1293 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));CHKERRQ(ierr); 1294 ierr = PetscViewerFlush(viewer); 1295 1296 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1297 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));CHKERRQ(ierr); 1298 ierr = PetscViewerFlush(viewer); 1299 1300 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1301 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));CHKERRQ(ierr); 1302 ierr = PetscViewerFlush(viewer); 1303 1304 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1305 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));CHKERRQ(ierr); 1306 ierr = PetscViewerFlush(viewer); 1307 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1308 1309 if (!mumps->myid) { /* information from the host */ 1310 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));CHKERRQ(ierr); 1311 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));CHKERRQ(ierr); 1312 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));CHKERRQ(ierr); 1313 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); 1314 1315 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));CHKERRQ(ierr); 1316 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));CHKERRQ(ierr); 1317 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));CHKERRQ(ierr); 1318 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));CHKERRQ(ierr); 1319 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));CHKERRQ(ierr); 1320 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));CHKERRQ(ierr); 1321 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));CHKERRQ(ierr); 1322 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));CHKERRQ(ierr); 1323 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));CHKERRQ(ierr); 1324 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));CHKERRQ(ierr); 1325 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));CHKERRQ(ierr); 1326 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));CHKERRQ(ierr); 1327 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));CHKERRQ(ierr); 1328 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); 1329 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); 1330 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); 1331 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); 1332 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));CHKERRQ(ierr); 1333 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); 1334 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); 1335 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));CHKERRQ(ierr); 1336 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));CHKERRQ(ierr); 1337 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));CHKERRQ(ierr); 1338 ierr = PetscViewerASCIIPrintf(viewer," INFOG(28) (after factorization: number of null pivots encountered): %d\n",mumps->id.INFOG(28));CHKERRQ(ierr); 1339 ierr = PetscViewerASCIIPrintf(viewer," INFOG(29) (after factorization: effective number of entries in the factors (sum over all processors)): %d\n",mumps->id.INFOG(29));CHKERRQ(ierr); 1340 ierr = PetscViewerASCIIPrintf(viewer," INFOG(30, 31) (after solution: size in Mbytes of memory used during solution phase): %d, %d\n",mumps->id.INFOG(30),mumps->id.INFOG(31));CHKERRQ(ierr); 1341 ierr = PetscViewerASCIIPrintf(viewer," INFOG(32) (after analysis: type of analysis done): %d\n",mumps->id.INFOG(32));CHKERRQ(ierr); 1342 ierr = PetscViewerASCIIPrintf(viewer," INFOG(33) (value used for ICNTL(8)): %d\n",mumps->id.INFOG(33));CHKERRQ(ierr); 1343 ierr = PetscViewerASCIIPrintf(viewer," INFOG(34) (exponent of the determinant if determinant is requested): %d\n",mumps->id.INFOG(34));CHKERRQ(ierr); 1344 } 1345 } 1346 } 1347 PetscFunctionReturn(0); 1348 } 1349 1350 #undef __FUNCT__ 1351 #define __FUNCT__ "MatGetInfo_MUMPS" 1352 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1353 { 1354 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1355 1356 PetscFunctionBegin; 1357 info->block_size = 1.0; 1358 info->nz_allocated = mumps->id.INFOG(20); 1359 info->nz_used = mumps->id.INFOG(20); 1360 info->nz_unneeded = 0.0; 1361 info->assemblies = 0.0; 1362 info->mallocs = 0.0; 1363 info->memory = 0.0; 1364 info->fill_ratio_given = 0; 1365 info->fill_ratio_needed = 0; 1366 info->factor_mallocs = 0; 1367 PetscFunctionReturn(0); 1368 } 1369 1370 /* -------------------------------------------------------------------------------------------*/ 1371 #undef __FUNCT__ 1372 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1373 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1374 { 1375 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1376 1377 PetscFunctionBegin; 1378 mumps->id.ICNTL(icntl) = ival; 1379 PetscFunctionReturn(0); 1380 } 1381 1382 #undef __FUNCT__ 1383 #define __FUNCT__ "MatMumpsGetIcntl_MUMPS" 1384 PetscErrorCode MatMumpsGetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt *ival) 1385 { 1386 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1387 1388 PetscFunctionBegin; 1389 *ival = mumps->id.ICNTL(icntl); 1390 PetscFunctionReturn(0); 1391 } 1392 1393 #undef __FUNCT__ 1394 #define __FUNCT__ "MatMumpsSetIcntl" 1395 /*@ 1396 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1397 1398 Logically Collective on Mat 1399 1400 Input Parameters: 1401 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1402 . icntl - index of MUMPS parameter array ICNTL() 1403 - ival - value of MUMPS ICNTL(icntl) 1404 1405 Options Database: 1406 . -mat_mumps_icntl_<icntl> <ival> 1407 1408 Level: beginner 1409 1410 References: MUMPS Users' Guide 1411 1412 .seealso: MatGetFactor() 1413 @*/ 1414 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1415 { 1416 PetscErrorCode ierr; 1417 1418 PetscFunctionBegin; 1419 PetscValidLogicalCollectiveInt(F,icntl,2); 1420 PetscValidLogicalCollectiveInt(F,ival,3); 1421 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1422 PetscFunctionReturn(0); 1423 } 1424 1425 #undef __FUNCT__ 1426 #define __FUNCT__ "MatMumpsGetIcntl" 1427 /*@ 1428 MatMumpsGetIcntl - Get MUMPS parameter ICNTL() 1429 1430 Logically Collective on Mat 1431 1432 Input Parameters: 1433 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1434 - icntl - index of MUMPS parameter array ICNTL() 1435 1436 Output Parameter: 1437 . ival - value of MUMPS ICNTL(icntl) 1438 1439 Level: beginner 1440 1441 References: MUMPS Users' Guide 1442 1443 .seealso: MatGetFactor() 1444 @*/ 1445 PetscErrorCode MatMumpsGetIcntl(Mat F,PetscInt icntl,PetscInt *ival) 1446 { 1447 PetscErrorCode ierr; 1448 1449 PetscFunctionBegin; 1450 PetscValidLogicalCollectiveInt(F,icntl,2); 1451 PetscValidIntPointer(ival,3); 1452 ierr = PetscTryMethod(F,"MatMumpsGetIcntl_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1453 PetscFunctionReturn(0); 1454 } 1455 1456 /* -------------------------------------------------------------------------------------------*/ 1457 #undef __FUNCT__ 1458 #define __FUNCT__ "MatMumpsSetCntl_MUMPS" 1459 PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val) 1460 { 1461 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1462 1463 PetscFunctionBegin; 1464 mumps->id.CNTL(icntl) = val; 1465 PetscFunctionReturn(0); 1466 } 1467 1468 #undef __FUNCT__ 1469 #define __FUNCT__ "MatMumpsGetCntl_MUMPS" 1470 PetscErrorCode MatMumpsGetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal *val) 1471 { 1472 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1473 1474 PetscFunctionBegin; 1475 *val = mumps->id.CNTL(icntl); 1476 PetscFunctionReturn(0); 1477 } 1478 1479 #undef __FUNCT__ 1480 #define __FUNCT__ "MatMumpsSetCntl" 1481 /*@ 1482 MatMumpsSetCntl - Set MUMPS parameter CNTL() 1483 1484 Logically Collective on Mat 1485 1486 Input Parameters: 1487 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1488 . icntl - index of MUMPS parameter array CNTL() 1489 - val - value of MUMPS CNTL(icntl) 1490 1491 Options Database: 1492 . -mat_mumps_cntl_<icntl> <val> 1493 1494 Level: beginner 1495 1496 References: MUMPS Users' Guide 1497 1498 .seealso: MatGetFactor() 1499 @*/ 1500 PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val) 1501 { 1502 PetscErrorCode ierr; 1503 1504 PetscFunctionBegin; 1505 PetscValidLogicalCollectiveInt(F,icntl,2); 1506 PetscValidLogicalCollectiveReal(F,val,3); 1507 ierr = PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));CHKERRQ(ierr); 1508 PetscFunctionReturn(0); 1509 } 1510 1511 #undef __FUNCT__ 1512 #define __FUNCT__ "MatMumpsGetCntl" 1513 /*@ 1514 MatMumpsGetCntl - Get MUMPS parameter CNTL() 1515 1516 Logically Collective on Mat 1517 1518 Input Parameters: 1519 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1520 - icntl - index of MUMPS parameter array CNTL() 1521 1522 Output Parameter: 1523 . val - value of MUMPS CNTL(icntl) 1524 1525 Level: beginner 1526 1527 References: MUMPS Users' Guide 1528 1529 .seealso: MatGetFactor() 1530 @*/ 1531 PetscErrorCode MatMumpsGetCntl(Mat F,PetscInt icntl,PetscReal *val) 1532 { 1533 PetscErrorCode ierr; 1534 1535 PetscFunctionBegin; 1536 PetscValidLogicalCollectiveInt(F,icntl,2); 1537 PetscValidRealPointer(val,3); 1538 ierr = PetscTryMethod(F,"MatMumpsGetCntl_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1539 PetscFunctionReturn(0); 1540 } 1541 1542 #undef __FUNCT__ 1543 #define __FUNCT__ "MatMumpsGetInfo_MUMPS" 1544 PetscErrorCode MatMumpsGetInfo_MUMPS(Mat F,PetscInt icntl,PetscInt *info) 1545 { 1546 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1547 1548 PetscFunctionBegin; 1549 *info = mumps->id.INFO(icntl); 1550 PetscFunctionReturn(0); 1551 } 1552 1553 #undef __FUNCT__ 1554 #define __FUNCT__ "MatMumpsGetInfog_MUMPS" 1555 PetscErrorCode MatMumpsGetInfog_MUMPS(Mat F,PetscInt icntl,PetscInt *infog) 1556 { 1557 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1558 1559 PetscFunctionBegin; 1560 *infog = mumps->id.INFOG(icntl); 1561 PetscFunctionReturn(0); 1562 } 1563 1564 #undef __FUNCT__ 1565 #define __FUNCT__ "MatMumpsGetRinfo_MUMPS" 1566 PetscErrorCode MatMumpsGetRinfo_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfo) 1567 { 1568 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1569 1570 PetscFunctionBegin; 1571 *rinfo = mumps->id.RINFO(icntl); 1572 PetscFunctionReturn(0); 1573 } 1574 1575 #undef __FUNCT__ 1576 #define __FUNCT__ "MatMumpsGetRinfog_MUMPS" 1577 PetscErrorCode MatMumpsGetRinfog_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfog) 1578 { 1579 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1580 1581 PetscFunctionBegin; 1582 *rinfog = mumps->id.RINFOG(icntl); 1583 PetscFunctionReturn(0); 1584 } 1585 1586 #undef __FUNCT__ 1587 #define __FUNCT__ "MatMumpsGetInfo" 1588 /*@ 1589 MatMumpsGetInfo - Get MUMPS parameter INFO() 1590 1591 Logically Collective on Mat 1592 1593 Input Parameters: 1594 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1595 - icntl - index of MUMPS parameter array INFO() 1596 1597 Output Parameter: 1598 . ival - value of MUMPS INFO(icntl) 1599 1600 Level: beginner 1601 1602 References: MUMPS Users' Guide 1603 1604 .seealso: MatGetFactor() 1605 @*/ 1606 PetscErrorCode MatMumpsGetInfo(Mat F,PetscInt icntl,PetscInt *ival) 1607 { 1608 PetscErrorCode ierr; 1609 1610 PetscFunctionBegin; 1611 PetscValidIntPointer(ival,3); 1612 ierr = PetscTryMethod(F,"MatMumpsGetInfo_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1613 PetscFunctionReturn(0); 1614 } 1615 1616 #undef __FUNCT__ 1617 #define __FUNCT__ "MatMumpsGetInfog" 1618 /*@ 1619 MatMumpsGetInfog - Get MUMPS parameter INFOG() 1620 1621 Logically Collective on Mat 1622 1623 Input Parameters: 1624 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1625 - icntl - index of MUMPS parameter array INFOG() 1626 1627 Output Parameter: 1628 . ival - value of MUMPS INFOG(icntl) 1629 1630 Level: beginner 1631 1632 References: MUMPS Users' Guide 1633 1634 .seealso: MatGetFactor() 1635 @*/ 1636 PetscErrorCode MatMumpsGetInfog(Mat F,PetscInt icntl,PetscInt *ival) 1637 { 1638 PetscErrorCode ierr; 1639 1640 PetscFunctionBegin; 1641 PetscValidIntPointer(ival,3); 1642 ierr = PetscTryMethod(F,"MatMumpsGetInfog_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1643 PetscFunctionReturn(0); 1644 } 1645 1646 #undef __FUNCT__ 1647 #define __FUNCT__ "MatMumpsGetRinfo" 1648 /*@ 1649 MatMumpsGetRinfo - Get MUMPS parameter RINFO() 1650 1651 Logically Collective on Mat 1652 1653 Input Parameters: 1654 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1655 - icntl - index of MUMPS parameter array RINFO() 1656 1657 Output Parameter: 1658 . val - value of MUMPS RINFO(icntl) 1659 1660 Level: beginner 1661 1662 References: MUMPS Users' Guide 1663 1664 .seealso: MatGetFactor() 1665 @*/ 1666 PetscErrorCode MatMumpsGetRinfo(Mat F,PetscInt icntl,PetscReal *val) 1667 { 1668 PetscErrorCode ierr; 1669 1670 PetscFunctionBegin; 1671 PetscValidRealPointer(val,3); 1672 ierr = PetscTryMethod(F,"MatMumpsGetRinfo_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1673 PetscFunctionReturn(0); 1674 } 1675 1676 #undef __FUNCT__ 1677 #define __FUNCT__ "MatMumpsGetRinfog" 1678 /*@ 1679 MatMumpsGetRinfog - Get MUMPS parameter RINFOG() 1680 1681 Logically Collective on Mat 1682 1683 Input Parameters: 1684 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1685 - icntl - index of MUMPS parameter array RINFOG() 1686 1687 Output Parameter: 1688 . val - value of MUMPS RINFOG(icntl) 1689 1690 Level: beginner 1691 1692 References: MUMPS Users' Guide 1693 1694 .seealso: MatGetFactor() 1695 @*/ 1696 PetscErrorCode MatMumpsGetRinfog(Mat F,PetscInt icntl,PetscReal *val) 1697 { 1698 PetscErrorCode ierr; 1699 1700 PetscFunctionBegin; 1701 PetscValidRealPointer(val,3); 1702 ierr = PetscTryMethod(F,"MatMumpsGetRinfog_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1703 PetscFunctionReturn(0); 1704 } 1705 1706 /*MC 1707 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1708 distributed and sequential matrices via the external package MUMPS. 1709 1710 Works with MATAIJ and MATSBAIJ matrices 1711 1712 Options Database Keys: 1713 + -mat_mumps_icntl_1 <6>: ICNTL(1): output stream for error messages (None) 1714 . -mat_mumps_icntl_2 <0>: ICNTL(2): output stream for diagnostic printing, statistics, and warning (None) 1715 . -mat_mumps_icntl_3 <0>: ICNTL(3): output stream for global information, collected on the host (None) 1716 . -mat_mumps_icntl_4 <0>: ICNTL(4): level of printing (0 to 4) (None) 1717 . -mat_mumps_icntl_6 <7>: ICNTL(6): permutes to a zero-free diagonal and/or scale the matrix (0 to 7) (None) 1718 . -mat_mumps_icntl_7 <7>: ICNTL(7): computes a symmetric permutation in sequential analysis (0 to 7). 3=Scotch, 4=PORD, 5=Metis (None) 1719 . -mat_mumps_icntl_8 <77>: ICNTL(8): scaling strategy (-2 to 8 or 77) (None) 1720 . -mat_mumps_icntl_10 <0>: ICNTL(10): max num of refinements (None) 1721 . -mat_mumps_icntl_11 <0>: ICNTL(11): statistics related to an error analysis (via -ksp_view) (None) 1722 . -mat_mumps_icntl_12 <1>: ICNTL(12): an ordering strategy for symmetric matrices (0 to 3) (None) 1723 . -mat_mumps_icntl_13 <0>: ICNTL(13): parallelism of the root node (enable ScaLAPACK) and its splitting (None) 1724 . -mat_mumps_icntl_14 <20>: ICNTL(14): percentage increase in the estimated working space (None) 1725 . -mat_mumps_icntl_19 <0>: ICNTL(19): computes the Schur complement (None) 1726 . -mat_mumps_icntl_22 <0>: ICNTL(22): in-core/out-of-core factorization and solve (0 or 1) (None) 1727 . -mat_mumps_icntl_23 <0>: ICNTL(23): max size of the working memory (MB) that can allocate per processor (None) 1728 . -mat_mumps_icntl_24 <0>: ICNTL(24): detection of null pivot rows (0 or 1) (None) 1729 . -mat_mumps_icntl_25 <0>: ICNTL(25): compute a solution of a deficient matrix and a null space basis (None) 1730 . -mat_mumps_icntl_26 <0>: ICNTL(26): drives the solution phase if a Schur complement matrix (None) 1731 . -mat_mumps_icntl_28 <1>: ICNTL(28): use 1 for sequential analysis and ictnl(7) ordering, or 2 for parallel analysis and ictnl(29) ordering (None) 1732 . -mat_mumps_icntl_29 <0>: ICNTL(29): parallel ordering 1 = ptscotch, 2 = parmetis (None) 1733 . -mat_mumps_icntl_30 <0>: ICNTL(30): compute user-specified set of entries in inv(A) (None) 1734 . -mat_mumps_icntl_31 <0>: ICNTL(31): indicates which factors may be discarded during factorization (None) 1735 . -mat_mumps_icntl_33 <0>: ICNTL(33): compute determinant (None) 1736 . -mat_mumps_cntl_1 <0.01>: CNTL(1): relative pivoting threshold (None) 1737 . -mat_mumps_cntl_2 <1.49012e-08>: CNTL(2): stopping criterion of refinement (None) 1738 . -mat_mumps_cntl_3 <0>: CNTL(3): absolute pivoting threshold (None) 1739 . -mat_mumps_cntl_4 <-1>: CNTL(4): value for static pivoting (None) 1740 - -mat_mumps_cntl_5 <0>: CNTL(5): fixation for null pivots (None) 1741 1742 Level: beginner 1743 1744 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1745 1746 M*/ 1747 1748 #undef __FUNCT__ 1749 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1750 static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1751 { 1752 PetscFunctionBegin; 1753 *type = MATSOLVERMUMPS; 1754 PetscFunctionReturn(0); 1755 } 1756 1757 /* MatGetFactor for Seq and MPI AIJ matrices */ 1758 #undef __FUNCT__ 1759 #define __FUNCT__ "MatGetFactor_aij_mumps" 1760 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1761 { 1762 Mat B; 1763 PetscErrorCode ierr; 1764 Mat_MUMPS *mumps; 1765 PetscBool isSeqAIJ; 1766 1767 PetscFunctionBegin; 1768 /* Create the factorization matrix */ 1769 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1770 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1771 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1772 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1773 if (isSeqAIJ) { 1774 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1775 } else { 1776 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1777 } 1778 1779 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1780 1781 B->ops->view = MatView_MUMPS; 1782 B->ops->getinfo = MatGetInfo_MUMPS; 1783 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1784 1785 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1786 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1787 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1788 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1789 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1790 1791 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1792 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1793 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1794 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1795 if (ftype == MAT_FACTOR_LU) { 1796 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1797 B->factortype = MAT_FACTOR_LU; 1798 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1799 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1800 mumps->sym = 0; 1801 } else { 1802 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1803 B->factortype = MAT_FACTOR_CHOLESKY; 1804 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1805 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1806 if (A->spd_set && A->spd) mumps->sym = 1; 1807 else mumps->sym = 2; 1808 } 1809 1810 mumps->isAIJ = PETSC_TRUE; 1811 mumps->Destroy = B->ops->destroy; 1812 B->ops->destroy = MatDestroy_MUMPS; 1813 B->spptr = (void*)mumps; 1814 1815 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1816 1817 *F = B; 1818 PetscFunctionReturn(0); 1819 } 1820 1821 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1822 #undef __FUNCT__ 1823 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1824 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1825 { 1826 Mat B; 1827 PetscErrorCode ierr; 1828 Mat_MUMPS *mumps; 1829 PetscBool isSeqSBAIJ; 1830 1831 PetscFunctionBegin; 1832 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1833 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"); 1834 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1835 /* Create the factorization matrix */ 1836 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1837 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1838 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1839 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1840 if (isSeqSBAIJ) { 1841 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 1842 1843 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1844 } else { 1845 ierr = MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);CHKERRQ(ierr); 1846 1847 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1848 } 1849 1850 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1851 B->ops->view = MatView_MUMPS; 1852 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1853 1854 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1855 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1856 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1857 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1858 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1859 1860 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1861 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1862 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1863 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1864 1865 B->factortype = MAT_FACTOR_CHOLESKY; 1866 if (A->spd_set && A->spd) mumps->sym = 1; 1867 else mumps->sym = 2; 1868 1869 mumps->isAIJ = PETSC_FALSE; 1870 mumps->Destroy = B->ops->destroy; 1871 B->ops->destroy = MatDestroy_MUMPS; 1872 B->spptr = (void*)mumps; 1873 1874 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1875 1876 *F = B; 1877 PetscFunctionReturn(0); 1878 } 1879 1880 #undef __FUNCT__ 1881 #define __FUNCT__ "MatGetFactor_baij_mumps" 1882 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 1883 { 1884 Mat B; 1885 PetscErrorCode ierr; 1886 Mat_MUMPS *mumps; 1887 PetscBool isSeqBAIJ; 1888 1889 PetscFunctionBegin; 1890 /* Create the factorization matrix */ 1891 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1892 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1893 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1894 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1895 if (isSeqBAIJ) { 1896 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);CHKERRQ(ierr); 1897 } else { 1898 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 1899 } 1900 1901 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1902 if (ftype == MAT_FACTOR_LU) { 1903 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1904 B->factortype = MAT_FACTOR_LU; 1905 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 1906 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 1907 mumps->sym = 0; 1908 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 1909 1910 B->ops->view = MatView_MUMPS; 1911 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1912 1913 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1914 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1915 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1916 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1917 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1918 1919 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1920 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1921 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1922 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1923 1924 mumps->isAIJ = PETSC_TRUE; 1925 mumps->Destroy = B->ops->destroy; 1926 B->ops->destroy = MatDestroy_MUMPS; 1927 B->spptr = (void*)mumps; 1928 1929 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1930 1931 *F = B; 1932 PetscFunctionReturn(0); 1933 } 1934 1935 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 1936 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat,MatFactorType,Mat*); 1937 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 1938 1939 #undef __FUNCT__ 1940 #define __FUNCT__ "MatSolverPackageRegister_MUMPS" 1941 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MUMPS(void) 1942 { 1943 PetscErrorCode ierr; 1944 1945 PetscFunctionBegin; 1946 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 1947 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 1948 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 1949 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 1950 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPISBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 1951 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 1952 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 1953 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 1954 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 1955 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 1956 PetscFunctionReturn(0); 1957 } 1958 1959