1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the MUMPS sparse solver 5 */ 6 #include "../src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 7 #include "../src/mat/impls/aij/mpi/mpiaij.h" 8 #include "../src/mat/impls/sbaij/seq/sbaij.h" 9 #include "../src/mat/impls/sbaij/mpi/mpisbaij.h" 10 #include "../src/mat/impls/baij/seq/baij.h" 11 #include "../src/mat/impls/baij/mpi/mpibaij.h" 12 13 EXTERN_C_BEGIN 14 #if defined(PETSC_USE_COMPLEX) 15 #include "zmumps_c.h" 16 #else 17 #include "dmumps_c.h" 18 #endif 19 EXTERN_C_END 20 #define JOB_INIT -1 21 #define JOB_FACTSYMBOLIC 1 22 #define JOB_FACTNUMERIC 2 23 #define JOB_SOLVE 3 24 #define JOB_END -2 25 26 27 /* macros s.t. indices match MUMPS documentation */ 28 #define ICNTL(I) icntl[(I)-1] 29 #define CNTL(I) cntl[(I)-1] 30 #define INFOG(I) infog[(I)-1] 31 #define INFO(I) info[(I)-1] 32 #define RINFOG(I) rinfog[(I)-1] 33 #define RINFO(I) rinfo[(I)-1] 34 35 typedef struct { 36 #if defined(PETSC_USE_COMPLEX) 37 ZMUMPS_STRUC_C id; 38 #else 39 DMUMPS_STRUC_C id; 40 #endif 41 MatStructure matstruc; 42 PetscMPIInt myid,size; 43 PetscInt *irn,*jcn,nz,sym,nSolve; 44 PetscScalar *val; 45 MPI_Comm comm_mumps; 46 VecScatter scat_rhs, scat_sol; 47 PetscBool isAIJ,CleanUpMUMPS,mumpsview; 48 Vec b_seq,x_seq; 49 PetscErrorCode (*MatDestroy)(Mat); 50 PetscErrorCode (*ConvertToTriples)(Mat, int, MatReuse, int*, int**, int**, PetscScalar**); 51 } Mat_MUMPS; 52 53 extern PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*); 54 55 56 /* MatConvertToTriples_A_B */ 57 /*convert Petsc matrix to triples: row[nz], col[nz], val[nz] */ 58 /* 59 input: 60 A - matrix in aij,baij or sbaij (bs=1) format 61 shift - 0: C style output triple; 1: Fortran style output triple. 62 reuse - MAT_INITIAL_MATRIX: spaces are allocated and values are set for the triple 63 MAT_REUSE_MATRIX: only the values in v array are updated 64 output: 65 nnz - dim of r, c, and v (number of local nonzero entries of A) 66 r, c, v - row and col index, matrix values (matrix triples) 67 */ 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "MatConvertToTriples_seqaij_seqaij" 71 PetscErrorCode MatConvertToTriples_seqaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 72 { 73 const PetscInt *ai,*aj,*ajj,M=A->rmap->n; 74 PetscInt nz,rnz,i,j; 75 PetscErrorCode ierr; 76 PetscInt *row,*col; 77 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 78 79 PetscFunctionBegin; 80 *v=aa->a; 81 if (reuse == MAT_INITIAL_MATRIX){ 82 nz = aa->nz; ai = aa->i; aj = aa->j; 83 *nnz = nz; 84 ierr = PetscMalloc(2*nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 85 col = row + nz; 86 87 nz = 0; 88 for(i=0; i<M; i++) { 89 rnz = ai[i+1] - ai[i]; 90 ajj = aj + ai[i]; 91 for(j=0; j<rnz; j++) { 92 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 93 } 94 } 95 *r = row; *c = col; 96 } 97 PetscFunctionReturn(0); 98 } 99 100 #undef __FUNCT__ 101 #define __FUNCT__ "MatConvertToTriples_seqbaij_seqaij" 102 PetscErrorCode MatConvertToTriples_seqbaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 103 { 104 Mat_SeqBAIJ *aa=(Mat_SeqBAIJ*)A->data; 105 const PetscInt *ai,*aj,*ajj,bs=A->rmap->bs,bs2=aa->bs2,M=A->rmap->N/bs; 106 PetscInt nz,idx=0,rnz,i,j,k,m,ii; 107 PetscErrorCode ierr; 108 PetscInt *row,*col; 109 110 PetscFunctionBegin; 111 *v = aa->a; 112 if (reuse == MAT_INITIAL_MATRIX){ 113 ai = aa->i; aj = aa->j; 114 nz = bs2*aa->nz; 115 *nnz = nz; 116 ierr = PetscMalloc(2*nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 117 col = row + nz; 118 119 for(i=0; i<M; i++) { 120 ii = 0; 121 ajj = aj + ai[i]; 122 rnz = ai[i+1] - ai[i]; 123 for(k=0; k<rnz; k++) { 124 for(j=0; j<bs; j++) { 125 for(m=0; m<bs; m++) { 126 row[idx] = i*bs + m + shift; 127 col[idx++] = bs*(ajj[k]) + j + shift; 128 } 129 } 130 } 131 } 132 *r = row; *c = col; 133 } 134 PetscFunctionReturn(0); 135 } 136 137 #undef __FUNCT__ 138 #define __FUNCT__ "MatConvertToTriples_seqsbaij_seqsbaij" 139 PetscErrorCode MatConvertToTriples_seqsbaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 140 { 141 const PetscInt *ai, *aj,*ajj,M=A->rmap->n; 142 PetscInt nz,rnz,i,j; 143 PetscErrorCode ierr; 144 PetscInt *row,*col; 145 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data; 146 147 PetscFunctionBegin; 148 if (reuse == MAT_INITIAL_MATRIX){ 149 nz = aa->nz;ai=aa->i; aj=aa->j;*v=aa->a; 150 *nnz = nz; 151 ierr = PetscMalloc(2*nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 152 col = row + nz; 153 154 nz = 0; 155 for(i=0; i<M; i++) { 156 rnz = ai[i+1] - ai[i]; 157 ajj = aj + ai[i]; 158 for(j=0; j<rnz; j++) { 159 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 160 } 161 } 162 *r = row; *c = col; 163 } 164 PetscFunctionReturn(0); 165 } 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "MatConvertToTriples_seqaij_seqsbaij" 169 PetscErrorCode MatConvertToTriples_seqaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 170 { 171 const PetscInt *ai,*aj,*ajj,*adiag,M=A->rmap->n; 172 PetscInt nz,rnz,i,j; 173 const PetscScalar *av,*v1; 174 PetscScalar *val; 175 PetscErrorCode ierr; 176 PetscInt *row,*col; 177 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data; 178 179 PetscFunctionBegin; 180 ai=aa->i; aj=aa->j;av=aa->a; 181 adiag=aa->diag; 182 if (reuse == MAT_INITIAL_MATRIX){ 183 nz = M + (aa->nz-M)/2; 184 *nnz = nz; 185 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 186 col = row + nz; 187 val = (PetscScalar*)(col + nz); 188 189 nz = 0; 190 for(i=0; i<M; i++) { 191 rnz = ai[i+1] - adiag[i]; 192 ajj = aj + adiag[i]; 193 v1 = av + adiag[i]; 194 for(j=0; j<rnz; j++) { 195 row[nz] = i+shift; col[nz] = ajj[j] + shift; val[nz++] = v1[j]; 196 } 197 } 198 *r = row; *c = col; *v = val; 199 } else { 200 nz = 0; val = *v; 201 for(i=0; i <M; i++) { 202 rnz = ai[i+1] - adiag[i]; 203 ajj = aj + adiag[i]; 204 v1 = av + adiag[i]; 205 for(j=0; j<rnz; j++) { 206 val[nz++] = v1[j]; 207 } 208 } 209 } 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "MatConvertToTriples_mpisbaij_mpisbaij" 215 PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 216 { 217 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 218 PetscErrorCode ierr; 219 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 220 PetscInt *row,*col; 221 const PetscScalar *av, *bv,*v1,*v2; 222 PetscScalar *val; 223 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 224 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)(mat->A)->data; 225 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 226 227 PetscFunctionBegin; 228 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 229 garray = mat->garray; 230 av=aa->a; bv=bb->a; 231 232 if (reuse == MAT_INITIAL_MATRIX){ 233 nz = aa->nz + bb->nz; 234 *nnz = nz; 235 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 236 col = row + nz; 237 val = (PetscScalar*)(col + nz); 238 239 *r = row; *c = col; *v = val; 240 } else { 241 row = *r; col = *c; val = *v; 242 } 243 244 jj = 0; irow = rstart; 245 for ( i=0; i<m; i++ ) { 246 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 247 countA = ai[i+1] - ai[i]; 248 countB = bi[i+1] - bi[i]; 249 bjj = bj + bi[i]; 250 v1 = av + ai[i]; 251 v2 = bv + bi[i]; 252 253 /* A-part */ 254 for (j=0; j<countA; j++){ 255 if (reuse == MAT_INITIAL_MATRIX) { 256 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 257 } 258 val[jj++] = v1[j]; 259 } 260 261 /* B-part */ 262 for(j=0; j < countB; j++){ 263 if (reuse == MAT_INITIAL_MATRIX) { 264 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 265 } 266 val[jj++] = v2[j]; 267 } 268 irow++; 269 } 270 PetscFunctionReturn(0); 271 } 272 273 #undef __FUNCT__ 274 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpiaij" 275 PetscErrorCode MatConvertToTriples_mpiaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 276 { 277 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 278 PetscErrorCode ierr; 279 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 280 PetscInt *row,*col; 281 const PetscScalar *av, *bv,*v1,*v2; 282 PetscScalar *val; 283 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 284 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 285 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 286 287 PetscFunctionBegin; 288 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 289 garray = mat->garray; 290 av=aa->a; bv=bb->a; 291 292 if (reuse == MAT_INITIAL_MATRIX){ 293 nz = aa->nz + bb->nz; 294 *nnz = nz; 295 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 296 col = row + nz; 297 val = (PetscScalar*)(col + nz); 298 299 *r = row; *c = col; *v = val; 300 } else { 301 row = *r; col = *c; val = *v; 302 } 303 304 jj = 0; irow = rstart; 305 for ( i=0; i<m; i++ ) { 306 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 307 countA = ai[i+1] - ai[i]; 308 countB = bi[i+1] - bi[i]; 309 bjj = bj + bi[i]; 310 v1 = av + ai[i]; 311 v2 = bv + bi[i]; 312 313 /* A-part */ 314 for (j=0; j<countA; j++){ 315 if (reuse == MAT_INITIAL_MATRIX){ 316 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 317 } 318 val[jj++] = v1[j]; 319 } 320 321 /* B-part */ 322 for(j=0; j < countB; j++){ 323 if (reuse == MAT_INITIAL_MATRIX){ 324 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 325 } 326 val[jj++] = v2[j]; 327 } 328 irow++; 329 } 330 PetscFunctionReturn(0); 331 } 332 333 #undef __FUNCT__ 334 #define __FUNCT__ "MatConvertToTriples_mpibaij_mpiaij" 335 PetscErrorCode MatConvertToTriples_mpibaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 336 { 337 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data; 338 Mat_SeqBAIJ *aa=(Mat_SeqBAIJ*)(mat->A)->data; 339 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 340 const PetscInt *ai = aa->i, *bi = bb->i, *aj = aa->j, *bj = bb->j,*ajj, *bjj; 341 const PetscInt *garray = mat->garray,mbs=mat->mbs,rstartbs=mat->rstartbs; 342 const PetscInt bs = A->rmap->bs,bs2=mat->bs2; 343 PetscErrorCode ierr; 344 PetscInt nz,i,j,k,n,jj,irow,countA,countB,idx; 345 PetscInt *row,*col; 346 const PetscScalar *av=aa->a, *bv=bb->a,*v1,*v2; 347 PetscScalar *val; 348 349 PetscFunctionBegin; 350 351 if (reuse == MAT_INITIAL_MATRIX) { 352 nz = bs2*(aa->nz + bb->nz); 353 *nnz = nz; 354 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 355 col = row + nz; 356 val = (PetscScalar*)(col + nz); 357 358 *r = row; *c = col; *v = val; 359 } else { 360 row = *r; col = *c; val = *v; 361 } 362 363 jj = 0; irow = rstartbs; 364 for ( i=0; i<mbs; i++ ) { 365 countA = ai[i+1] - ai[i]; 366 countB = bi[i+1] - bi[i]; 367 ajj = aj + ai[i]; 368 bjj = bj + bi[i]; 369 v1 = av + bs2*ai[i]; 370 v2 = bv + bs2*bi[i]; 371 372 idx = 0; 373 /* A-part */ 374 for (k=0; k<countA; k++){ 375 for (j=0; j<bs; j++) { 376 for (n=0; n<bs; n++) { 377 if (reuse == MAT_INITIAL_MATRIX){ 378 row[jj] = bs*irow + n + shift; 379 col[jj] = bs*(rstartbs + ajj[k]) + j + shift; 380 } 381 val[jj++] = v1[idx++]; 382 } 383 } 384 } 385 386 idx = 0; 387 /* B-part */ 388 for(k=0; k<countB; k++){ 389 for (j=0; j<bs; j++) { 390 for (n=0; n<bs; n++) { 391 if (reuse == MAT_INITIAL_MATRIX){ 392 row[jj] = bs*irow + n + shift; 393 col[jj] = bs*(garray[bjj[k]]) + j + shift; 394 } 395 val[jj++] = bv[idx++]; 396 } 397 } 398 } 399 irow++; 400 } 401 PetscFunctionReturn(0); 402 } 403 404 #undef __FUNCT__ 405 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpisbaij" 406 PetscErrorCode MatConvertToTriples_mpiaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 407 { 408 const PetscInt *ai, *aj,*adiag, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 409 PetscErrorCode ierr; 410 PetscInt rstart,nz,nza,nzb_low,i,j,jj,irow,countA,countB; 411 PetscInt *row,*col; 412 const PetscScalar *av, *bv,*v1,*v2; 413 PetscScalar *val; 414 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 415 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 416 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 417 418 PetscFunctionBegin; 419 ai=aa->i; aj=aa->j; adiag=aa->diag; 420 bi=bb->i; bj=bb->j; garray = mat->garray; 421 av=aa->a; bv=bb->a; 422 rstart = A->rmap->rstart; 423 424 if (reuse == MAT_INITIAL_MATRIX) { 425 nza = 0;nzb_low = 0; 426 for(i=0; i<m; i++){ 427 nza = nza + (ai[i+1] - adiag[i]); 428 countB = bi[i+1] - bi[i]; 429 bjj = bj + bi[i]; 430 431 j = 0; 432 while(garray[bjj[j]] < rstart) { 433 if(j == countB) break; 434 j++;nzb_low++; 435 } 436 } 437 /* Total nz = nz for the upper triangular A part + nz for the 2nd B part */ 438 nz = nza + (bb->nz - nzb_low); 439 *nnz = nz; 440 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 441 col = row + nz; 442 val = (PetscScalar*)(col + nz); 443 444 *r = row; *c = col; *v = val; 445 } else { 446 row = *r; col = *c; val = *v; 447 } 448 449 jj = 0; irow = rstart; 450 for ( i=0; i<m; i++ ) { 451 ajj = aj + adiag[i]; /* ptr to the beginning of the diagonal of this row */ 452 v1 = av + adiag[i]; 453 countA = ai[i+1] - adiag[i]; 454 countB = bi[i+1] - bi[i]; 455 bjj = bj + bi[i]; 456 v2 = bv + bi[i]; 457 458 /* A-part */ 459 for (j=0; j<countA; j++){ 460 if (reuse == MAT_INITIAL_MATRIX) { 461 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 462 } 463 val[jj++] = v1[j]; 464 } 465 466 /* B-part */ 467 for(j=0; j < countB; j++){ 468 if (garray[bjj[j]] > rstart) { 469 if (reuse == MAT_INITIAL_MATRIX) { 470 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 471 } 472 val[jj++] = v2[j]; 473 } 474 } 475 irow++; 476 } 477 PetscFunctionReturn(0); 478 } 479 480 #undef __FUNCT__ 481 #define __FUNCT__ "MatDestroy_MUMPS" 482 PetscErrorCode MatDestroy_MUMPS(Mat A) 483 { 484 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 485 PetscErrorCode ierr; 486 487 PetscFunctionBegin; 488 if (lu->CleanUpMUMPS) { 489 /* Terminate instance, deallocate memories */ 490 if (lu->id.sol_loc){ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr);} 491 if (lu->scat_rhs){ierr = VecScatterDestroy(lu->scat_rhs);CHKERRQ(ierr);} 492 if (lu->b_seq) {ierr = VecDestroy(lu->b_seq);CHKERRQ(ierr);} 493 if (lu->nSolve && lu->scat_sol){ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);} 494 if (lu->nSolve && lu->x_seq){ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);} 495 496 ierr = PetscFree(lu->irn);CHKERRQ(ierr); 497 lu->id.job=JOB_END; 498 #if defined(PETSC_USE_COMPLEX) 499 zmumps_c(&lu->id); 500 #else 501 dmumps_c(&lu->id); 502 #endif 503 ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr); 504 } 505 /* clear composed functions */ 506 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatFactorGetSolverPackage_C","",PETSC_NULL);CHKERRQ(ierr); 507 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMumpsSetIcntl_C","",PETSC_NULL);CHKERRQ(ierr); 508 ierr = (lu->MatDestroy)(A);CHKERRQ(ierr); 509 PetscFunctionReturn(0); 510 } 511 512 #undef __FUNCT__ 513 #define __FUNCT__ "MatSolve_MUMPS" 514 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 515 { 516 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 517 PetscScalar *array; 518 Vec b_seq; 519 IS is_iden,is_petsc; 520 PetscErrorCode ierr; 521 PetscInt i; 522 523 PetscFunctionBegin; 524 lu->id.nrhs = 1; 525 b_seq = lu->b_seq; 526 if (lu->size > 1){ 527 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 528 ierr = VecScatterBegin(lu->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 529 ierr = VecScatterEnd(lu->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 530 if (!lu->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 531 } else { /* size == 1 */ 532 ierr = VecCopy(b,x);CHKERRQ(ierr); 533 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 534 } 535 if (!lu->myid) { /* define rhs on the host */ 536 lu->id.nrhs = 1; 537 #if defined(PETSC_USE_COMPLEX) 538 lu->id.rhs = (mumps_double_complex*)array; 539 #else 540 lu->id.rhs = array; 541 #endif 542 } 543 544 /* solve phase */ 545 /*-------------*/ 546 lu->id.job = JOB_SOLVE; 547 #if defined(PETSC_USE_COMPLEX) 548 zmumps_c(&lu->id); 549 #else 550 dmumps_c(&lu->id); 551 #endif 552 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)); 553 554 if (lu->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 555 if (!lu->nSolve){ /* create scatter scat_sol */ 556 ierr = ISCreateStride(PETSC_COMM_SELF,lu->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 557 for (i=0; i<lu->id.lsol_loc; i++){ 558 lu->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 559 } 560 ierr = ISCreateGeneral(PETSC_COMM_SELF,lu->id.lsol_loc,lu->id.isol_loc,PETSC_COPY_VALUES,&is_petsc);CHKERRQ(ierr); /* to */ 561 ierr = VecScatterCreate(lu->x_seq,is_iden,x,is_petsc,&lu->scat_sol);CHKERRQ(ierr); 562 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 563 ierr = ISDestroy(is_petsc);CHKERRQ(ierr); 564 } 565 ierr = VecScatterBegin(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 566 ierr = VecScatterEnd(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 567 } 568 lu->nSolve++; 569 PetscFunctionReturn(0); 570 } 571 572 #if !defined(PETSC_USE_COMPLEX) 573 /* 574 input: 575 F: numeric factor 576 output: 577 nneg: total number of negative pivots 578 nzero: 0 579 npos: (global dimension of F) - nneg 580 */ 581 582 #undef __FUNCT__ 583 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 584 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 585 { 586 Mat_MUMPS *lu =(Mat_MUMPS*)F->spptr; 587 PetscErrorCode ierr; 588 PetscMPIInt size; 589 590 PetscFunctionBegin; 591 ierr = MPI_Comm_size(((PetscObject)F)->comm,&size);CHKERRQ(ierr); 592 /* 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 */ 593 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)); 594 if (nneg){ 595 if (!lu->myid){ 596 *nneg = lu->id.INFOG(12); 597 } 598 ierr = MPI_Bcast(nneg,1,MPI_INT,0,lu->comm_mumps);CHKERRQ(ierr); 599 } 600 if (nzero) *nzero = 0; 601 if (npos) *npos = F->rmap->N - (*nneg); 602 PetscFunctionReturn(0); 603 } 604 #endif /* !defined(PETSC_USE_COMPLEX) */ 605 606 #undef __FUNCT__ 607 #define __FUNCT__ "MatFactorNumeric_MUMPS" 608 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 609 { 610 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 611 PetscErrorCode ierr; 612 MatReuse reuse; 613 Mat F_diag; 614 PetscBool isMPIAIJ; 615 616 PetscFunctionBegin; 617 reuse = MAT_REUSE_MATRIX; 618 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 619 620 /* numerical factorization phase */ 621 /*-------------------------------*/ 622 lu->id.job = JOB_FACTNUMERIC; 623 if(!lu->id.ICNTL(18)) { 624 if (!lu->myid) { 625 #if defined(PETSC_USE_COMPLEX) 626 lu->id.a = (mumps_double_complex*)lu->val; 627 #else 628 lu->id.a = lu->val; 629 #endif 630 } 631 } else { 632 #if defined(PETSC_USE_COMPLEX) 633 lu->id.a_loc = (mumps_double_complex*)lu->val; 634 #else 635 lu->id.a_loc = lu->val; 636 #endif 637 } 638 #if defined(PETSC_USE_COMPLEX) 639 zmumps_c(&lu->id); 640 #else 641 dmumps_c(&lu->id); 642 #endif 643 if (lu->id.INFOG(1) < 0) { 644 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)); 645 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)); 646 } 647 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)); 648 649 if (lu->size > 1){ 650 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 651 if(isMPIAIJ) { 652 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 653 } else { 654 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 655 } 656 F_diag->assembled = PETSC_TRUE; 657 if (lu->nSolve){ 658 ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr); 659 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 660 ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr); 661 } 662 } 663 (F)->assembled = PETSC_TRUE; 664 lu->matstruc = SAME_NONZERO_PATTERN; 665 lu->CleanUpMUMPS = PETSC_TRUE; 666 lu->nSolve = 0; 667 668 if (lu->size > 1){ 669 /* distributed solution */ 670 lu->id.ICNTL(21) = 1; 671 if (!lu->nSolve){ 672 /* Create x_seq=sol_loc for repeated use */ 673 PetscInt lsol_loc; 674 PetscScalar *sol_loc; 675 lsol_loc = lu->id.INFO(23); /* length of sol_loc */ 676 ierr = PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&lu->id.isol_loc);CHKERRQ(ierr); 677 lu->id.lsol_loc = lsol_loc; 678 #if defined(PETSC_USE_COMPLEX) 679 lu->id.sol_loc = (mumps_double_complex*)sol_loc; 680 #else 681 lu->id.sol_loc = sol_loc; 682 #endif 683 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,lsol_loc,sol_loc,&lu->x_seq);CHKERRQ(ierr); 684 } 685 } 686 PetscFunctionReturn(0); 687 } 688 689 #undef __FUNCT__ 690 #define __FUNCT__ "PetscSetMUMPSOptions" 691 PetscErrorCode PetscSetMUMPSOptions(Mat F, Mat A) 692 { 693 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 694 PetscErrorCode ierr; 695 PetscInt icntl; 696 PetscBool flg; 697 698 PetscFunctionBegin; 699 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 700 ierr = PetscOptionsBool("-mat_mumps_view","View MUMPS parameters","None",lu->mumpsview,&lu->mumpsview,PETSC_NULL);CHKERRQ(ierr); 701 if (lu->size == 1){ 702 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 703 } else { 704 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 705 } 706 707 icntl=-1; 708 lu->id.ICNTL(4) = 0; /* level of printing; overwrite mumps default ICNTL(4)=2 */ 709 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 710 if ((flg && icntl > 0) || PetscLogPrintInfo) { 711 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 712 } else { /* no output */ 713 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 714 lu->id.ICNTL(2) = 0; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 715 lu->id.ICNTL(3) = 0; /* output stream for global information, default=6 */ 716 } 717 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); 718 icntl=-1; 719 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 720 if (flg) { 721 if (icntl== 1){ 722 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"); 723 } else { 724 lu->id.ICNTL(7) = icntl; 725 } 726 } 727 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); 728 ierr = PetscOptionsInt("-mat_mumps_icntl_9","ICNTL(9): A or A^T x=b to be solved. 1: A; otherwise: A^T","None",lu->id.ICNTL(9),&lu->id.ICNTL(9),PETSC_NULL);CHKERRQ(ierr); 729 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); 730 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); 731 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); 732 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); 733 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); 734 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",lu->id.ICNTL(19),&lu->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr); 735 736 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); 737 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); 738 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); 739 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); 740 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); 741 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",lu->id.ICNTL(27),&lu->id.ICNTL(27),PETSC_NULL);CHKERRQ(ierr); 742 743 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 744 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); 745 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 746 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); 747 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); 748 PetscOptionsEnd(); 749 PetscFunctionReturn(0); 750 } 751 752 #undef __FUNCT__ 753 #define __FUNCT__ "PetscInitializeMUMPS" 754 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS* mumps) 755 { 756 PetscErrorCode ierr; 757 758 PetscFunctionBegin; 759 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &mumps->myid); 760 ierr = MPI_Comm_size(((PetscObject)A)->comm,&mumps->size);CHKERRQ(ierr); 761 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(mumps->comm_mumps));CHKERRQ(ierr); 762 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 763 764 mumps->id.job = JOB_INIT; 765 mumps->id.par = 1; /* host participates factorizaton and solve */ 766 mumps->id.sym = mumps->sym; 767 #if defined(PETSC_USE_COMPLEX) 768 zmumps_c(&mumps->id); 769 #else 770 dmumps_c(&mumps->id); 771 #endif 772 773 mumps->CleanUpMUMPS = PETSC_FALSE; 774 mumps->scat_rhs = PETSC_NULL; 775 mumps->scat_sol = PETSC_NULL; 776 mumps->nSolve = 0; 777 PetscFunctionReturn(0); 778 } 779 780 /* Note the Petsc r and c permutations are ignored */ 781 #undef __FUNCT__ 782 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 783 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 784 { 785 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 786 PetscErrorCode ierr; 787 MatReuse reuse; 788 Vec b; 789 IS is_iden; 790 const PetscInt M = A->rmap->N; 791 792 PetscFunctionBegin; 793 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 794 795 /* Set MUMPS options */ 796 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 797 798 reuse = MAT_INITIAL_MATRIX; 799 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 800 801 /* analysis phase */ 802 /*----------------*/ 803 lu->id.job = JOB_FACTSYMBOLIC; 804 lu->id.n = M; 805 switch (lu->id.ICNTL(18)){ 806 case 0: /* centralized assembled matrix input */ 807 if (!lu->myid) { 808 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 809 if (lu->id.ICNTL(6)>1){ 810 #if defined(PETSC_USE_COMPLEX) 811 lu->id.a = (mumps_double_complex*)lu->val; 812 #else 813 lu->id.a = lu->val; 814 #endif 815 } 816 } 817 break; 818 case 3: /* distributed assembled matrix input (size>1) */ 819 lu->id.nz_loc = lu->nz; 820 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 821 if (lu->id.ICNTL(6)>1) { 822 #if defined(PETSC_USE_COMPLEX) 823 lu->id.a_loc = (mumps_double_complex*)lu->val; 824 #else 825 lu->id.a_loc = lu->val; 826 #endif 827 } 828 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 829 if (!lu->myid){ 830 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 831 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 832 } else { 833 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 834 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 835 } 836 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 837 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 838 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 839 840 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 841 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 842 ierr = VecDestroy(b);CHKERRQ(ierr); 843 break; 844 } 845 #if defined(PETSC_USE_COMPLEX) 846 zmumps_c(&lu->id); 847 #else 848 dmumps_c(&lu->id); 849 #endif 850 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)); 851 852 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 853 F->ops->solve = MatSolve_MUMPS; 854 PetscFunctionReturn(0); 855 } 856 857 /* Note the Petsc r and c permutations are ignored */ 858 #undef __FUNCT__ 859 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 860 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 861 { 862 863 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 864 PetscErrorCode ierr; 865 MatReuse reuse; 866 Vec b; 867 IS is_iden; 868 const PetscInt M = A->rmap->N; 869 870 PetscFunctionBegin; 871 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 872 873 /* Set MUMPS options */ 874 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 875 876 reuse = MAT_INITIAL_MATRIX; 877 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 878 879 /* analysis phase */ 880 /*----------------*/ 881 lu->id.job = JOB_FACTSYMBOLIC; 882 lu->id.n = M; 883 switch (lu->id.ICNTL(18)){ 884 case 0: /* centralized assembled matrix input */ 885 if (!lu->myid) { 886 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 887 if (lu->id.ICNTL(6)>1){ 888 #if defined(PETSC_USE_COMPLEX) 889 lu->id.a = (mumps_double_complex*)lu->val; 890 #else 891 lu->id.a = lu->val; 892 #endif 893 } 894 } 895 break; 896 case 3: /* distributed assembled matrix input (size>1) */ 897 lu->id.nz_loc = lu->nz; 898 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 899 if (lu->id.ICNTL(6)>1) { 900 #if defined(PETSC_USE_COMPLEX) 901 lu->id.a_loc = (mumps_double_complex*)lu->val; 902 #else 903 lu->id.a_loc = lu->val; 904 #endif 905 } 906 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 907 if (!lu->myid){ 908 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 909 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 910 } else { 911 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 912 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 913 } 914 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 915 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 916 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 917 918 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 919 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 920 ierr = VecDestroy(b);CHKERRQ(ierr); 921 break; 922 } 923 #if defined(PETSC_USE_COMPLEX) 924 zmumps_c(&lu->id); 925 #else 926 dmumps_c(&lu->id); 927 #endif 928 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)); 929 930 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 931 F->ops->solve = MatSolve_MUMPS; 932 PetscFunctionReturn(0); 933 } 934 935 /* Note the Petsc r permutation is ignored */ 936 #undef __FUNCT__ 937 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 938 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 939 { 940 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 941 PetscErrorCode ierr; 942 MatReuse reuse; 943 Vec b; 944 IS is_iden; 945 const PetscInt M = A->rmap->N; 946 947 PetscFunctionBegin; 948 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 949 950 /* Set MUMPS options */ 951 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 952 953 reuse = MAT_INITIAL_MATRIX; 954 ierr = (*lu->ConvertToTriples)(A, 1 , reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 955 956 /* analysis phase */ 957 /*----------------*/ 958 lu->id.job = JOB_FACTSYMBOLIC; 959 lu->id.n = M; 960 switch (lu->id.ICNTL(18)){ 961 case 0: /* centralized assembled matrix input */ 962 if (!lu->myid) { 963 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 964 if (lu->id.ICNTL(6)>1){ 965 #if defined(PETSC_USE_COMPLEX) 966 lu->id.a = (mumps_double_complex*)lu->val; 967 #else 968 lu->id.a = lu->val; 969 #endif 970 } 971 } 972 break; 973 case 3: /* distributed assembled matrix input (size>1) */ 974 lu->id.nz_loc = lu->nz; 975 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 976 if (lu->id.ICNTL(6)>1) { 977 #if defined(PETSC_USE_COMPLEX) 978 lu->id.a_loc = (mumps_double_complex*)lu->val; 979 #else 980 lu->id.a_loc = lu->val; 981 #endif 982 } 983 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 984 if (!lu->myid){ 985 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 986 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 987 } else { 988 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 989 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 990 } 991 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 992 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 993 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 994 995 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 996 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 997 ierr = VecDestroy(b);CHKERRQ(ierr); 998 break; 999 } 1000 #if defined(PETSC_USE_COMPLEX) 1001 zmumps_c(&lu->id); 1002 #else 1003 dmumps_c(&lu->id); 1004 #endif 1005 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)); 1006 1007 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1008 F->ops->solve = MatSolve_MUMPS; 1009 #if !defined(PETSC_USE_COMPLEX) 1010 (F)->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1011 #endif 1012 PetscFunctionReturn(0); 1013 } 1014 1015 #undef __FUNCT__ 1016 #define __FUNCT__ "MatFactorInfo_MUMPS" 1017 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 1018 { 1019 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 1020 PetscErrorCode ierr; 1021 1022 PetscFunctionBegin; 1023 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 1024 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 1025 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 1026 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 1027 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 1028 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 1029 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 1030 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 1031 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 1032 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 1033 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 1034 if (lu->id.ICNTL(11)>0) { 1035 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 1036 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 1037 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 1038 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",lu->id.RINFOG(7),lu->id.RINFOG(8));CHKERRQ(ierr); 1039 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 1040 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 1041 1042 } 1043 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 1044 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 1045 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 1046 /* ICNTL(15-17) not used */ 1047 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 1048 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 1049 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 1050 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 1051 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 1052 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 1053 1054 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 1055 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 1056 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 1057 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 1058 1059 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 1060 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 1061 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 1062 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 1063 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 1064 1065 /* infomation local to each processor */ 1066 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1067 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 1068 ierr = PetscViewerFlush(viewer); 1069 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1070 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 1071 ierr = PetscViewerFlush(viewer); 1072 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1073 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 1074 ierr = PetscViewerFlush(viewer); 1075 1076 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1077 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 1078 ierr = PetscViewerFlush(viewer); 1079 1080 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1081 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 1082 ierr = PetscViewerFlush(viewer); 1083 1084 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1085 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 1086 ierr = PetscViewerFlush(viewer); 1087 1088 if (!lu->myid){ /* information from the host */ 1089 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 1090 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 1091 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 1092 1093 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 1094 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 1095 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 1096 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 1097 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 1098 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 1099 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 1100 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 1101 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 1102 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 1103 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 1104 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 1105 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 1106 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); 1107 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); 1108 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); 1109 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); 1110 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 1111 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); 1112 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); 1113 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 1114 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 1115 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 1116 } 1117 PetscFunctionReturn(0); 1118 } 1119 1120 #undef __FUNCT__ 1121 #define __FUNCT__ "MatView_MUMPS" 1122 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1123 { 1124 PetscErrorCode ierr; 1125 PetscBool iascii; 1126 PetscViewerFormat format; 1127 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1128 1129 PetscFunctionBegin; 1130 /* check if matrix is mumps type */ 1131 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1132 1133 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1134 if (iascii) { 1135 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1136 if (format == PETSC_VIEWER_ASCII_INFO || mumps->mumpsview){ 1137 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1138 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1139 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1140 if (mumps->mumpsview){ /* View all MUMPS parameters */ 1141 ierr = MatFactorInfo_MUMPS(A,viewer);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 Guide) 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