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