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