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