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