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; 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 if (lu->size == 1){ 701 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 702 } else { 703 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 704 } 705 706 icntl=-1; 707 lu->id.ICNTL(4) = 0; /* level of printing; overwrite mumps default ICNTL(4)=2 */ 708 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 709 if ((flg && icntl > 0) || PetscLogPrintInfo) { 710 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 711 } else { /* no output */ 712 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 713 lu->id.ICNTL(2) = 0; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 714 lu->id.ICNTL(3) = 0; /* output stream for global information, default=6 */ 715 } 716 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); 717 icntl=-1; 718 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): sequential matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 719 if (flg) { 720 if (icntl== 1){ 721 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"); 722 } else { 723 lu->id.ICNTL(7) = icntl; 724 } 725 } 726 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); 727 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); 728 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); 729 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); 730 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); 731 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); 732 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); 733 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",lu->id.ICNTL(19),&lu->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr); 734 735 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); 736 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); 737 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); 738 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); 739 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); 740 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",lu->id.ICNTL(27),&lu->id.ICNTL(27),PETSC_NULL);CHKERRQ(ierr); 741 ierr = PetscOptionsInt("-mat_mumps_icntl_28","ICNTL(28): use 1 sequential or 2 parallel ordering","None",lu->id.ICNTL(28),&lu->id.ICNTL(28),PETSC_NULL);CHKERRQ(ierr); 742 ierr = PetscOptionsInt("-mat_mumps_icntl_29","ICNTL(29): parallel ordering 2 for parmetis","None",lu->id.ICNTL(29),&lu->id.ICNTL(29),PETSC_NULL);CHKERRQ(ierr); 743 744 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 745 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); 746 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 747 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); 748 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); 749 PetscOptionsEnd(); 750 PetscFunctionReturn(0); 751 } 752 753 #undef __FUNCT__ 754 #define __FUNCT__ "PetscInitializeMUMPS" 755 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS* mumps) 756 { 757 PetscErrorCode ierr; 758 759 PetscFunctionBegin; 760 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &mumps->myid); 761 ierr = MPI_Comm_size(((PetscObject)A)->comm,&mumps->size);CHKERRQ(ierr); 762 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(mumps->comm_mumps));CHKERRQ(ierr); 763 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 764 765 mumps->id.job = JOB_INIT; 766 mumps->id.par = 1; /* host participates factorizaton and solve */ 767 mumps->id.sym = mumps->sym; 768 #if defined(PETSC_USE_COMPLEX) 769 zmumps_c(&mumps->id); 770 #else 771 dmumps_c(&mumps->id); 772 #endif 773 774 mumps->CleanUpMUMPS = PETSC_FALSE; 775 mumps->scat_rhs = PETSC_NULL; 776 mumps->scat_sol = PETSC_NULL; 777 mumps->nSolve = 0; 778 PetscFunctionReturn(0); 779 } 780 781 /* Note the Petsc r and c permutations are ignored */ 782 #undef __FUNCT__ 783 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 784 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 785 { 786 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 787 PetscErrorCode ierr; 788 MatReuse reuse; 789 Vec b; 790 IS is_iden; 791 const PetscInt M = A->rmap->N; 792 793 PetscFunctionBegin; 794 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 795 796 /* Set MUMPS options */ 797 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 798 799 reuse = MAT_INITIAL_MATRIX; 800 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 801 802 /* analysis phase */ 803 /*----------------*/ 804 lu->id.job = JOB_FACTSYMBOLIC; 805 lu->id.n = M; 806 switch (lu->id.ICNTL(18)){ 807 case 0: /* centralized assembled matrix input */ 808 if (!lu->myid) { 809 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 810 if (lu->id.ICNTL(6)>1){ 811 #if defined(PETSC_USE_COMPLEX) 812 lu->id.a = (mumps_double_complex*)lu->val; 813 #else 814 lu->id.a = lu->val; 815 #endif 816 } 817 } 818 break; 819 case 3: /* distributed assembled matrix input (size>1) */ 820 lu->id.nz_loc = lu->nz; 821 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 822 if (lu->id.ICNTL(6)>1) { 823 #if defined(PETSC_USE_COMPLEX) 824 lu->id.a_loc = (mumps_double_complex*)lu->val; 825 #else 826 lu->id.a_loc = lu->val; 827 #endif 828 } 829 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 830 if (!lu->myid){ 831 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 832 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 833 } else { 834 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 835 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 836 } 837 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 838 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 839 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 840 841 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 842 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 843 ierr = VecDestroy(b);CHKERRQ(ierr); 844 break; 845 } 846 #if defined(PETSC_USE_COMPLEX) 847 zmumps_c(&lu->id); 848 #else 849 dmumps_c(&lu->id); 850 #endif 851 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)); 852 853 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 854 F->ops->solve = MatSolve_MUMPS; 855 PetscFunctionReturn(0); 856 } 857 858 /* Note the Petsc r and c permutations are ignored */ 859 #undef __FUNCT__ 860 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 861 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 862 { 863 864 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 865 PetscErrorCode ierr; 866 MatReuse reuse; 867 Vec b; 868 IS is_iden; 869 const PetscInt M = A->rmap->N; 870 871 PetscFunctionBegin; 872 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 873 874 /* Set MUMPS options */ 875 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 876 877 reuse = MAT_INITIAL_MATRIX; 878 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 879 880 /* analysis phase */ 881 /*----------------*/ 882 lu->id.job = JOB_FACTSYMBOLIC; 883 lu->id.n = M; 884 switch (lu->id.ICNTL(18)){ 885 case 0: /* centralized assembled matrix input */ 886 if (!lu->myid) { 887 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 888 if (lu->id.ICNTL(6)>1){ 889 #if defined(PETSC_USE_COMPLEX) 890 lu->id.a = (mumps_double_complex*)lu->val; 891 #else 892 lu->id.a = lu->val; 893 #endif 894 } 895 } 896 break; 897 case 3: /* distributed assembled matrix input (size>1) */ 898 lu->id.nz_loc = lu->nz; 899 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 900 if (lu->id.ICNTL(6)>1) { 901 #if defined(PETSC_USE_COMPLEX) 902 lu->id.a_loc = (mumps_double_complex*)lu->val; 903 #else 904 lu->id.a_loc = lu->val; 905 #endif 906 } 907 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 908 if (!lu->myid){ 909 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 910 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 911 } else { 912 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 913 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 914 } 915 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 916 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 917 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 918 919 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 920 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 921 ierr = VecDestroy(b);CHKERRQ(ierr); 922 break; 923 } 924 #if defined(PETSC_USE_COMPLEX) 925 zmumps_c(&lu->id); 926 #else 927 dmumps_c(&lu->id); 928 #endif 929 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)); 930 931 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 932 F->ops->solve = MatSolve_MUMPS; 933 PetscFunctionReturn(0); 934 } 935 936 /* Note the Petsc r permutation is ignored */ 937 #undef __FUNCT__ 938 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 939 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 940 { 941 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 942 PetscErrorCode ierr; 943 MatReuse reuse; 944 Vec b; 945 IS is_iden; 946 const PetscInt M = A->rmap->N; 947 948 PetscFunctionBegin; 949 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 950 951 /* Set MUMPS options */ 952 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 953 954 reuse = MAT_INITIAL_MATRIX; 955 ierr = (*lu->ConvertToTriples)(A, 1 , reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 956 957 /* analysis phase */ 958 /*----------------*/ 959 lu->id.job = JOB_FACTSYMBOLIC; 960 lu->id.n = M; 961 switch (lu->id.ICNTL(18)){ 962 case 0: /* centralized assembled matrix input */ 963 if (!lu->myid) { 964 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 965 if (lu->id.ICNTL(6)>1){ 966 #if defined(PETSC_USE_COMPLEX) 967 lu->id.a = (mumps_double_complex*)lu->val; 968 #else 969 lu->id.a = lu->val; 970 #endif 971 } 972 } 973 break; 974 case 3: /* distributed assembled matrix input (size>1) */ 975 lu->id.nz_loc = lu->nz; 976 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 977 if (lu->id.ICNTL(6)>1) { 978 #if defined(PETSC_USE_COMPLEX) 979 lu->id.a_loc = (mumps_double_complex*)lu->val; 980 #else 981 lu->id.a_loc = lu->val; 982 #endif 983 } 984 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 985 if (!lu->myid){ 986 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 987 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 988 } else { 989 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 990 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 991 } 992 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 993 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 994 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 995 996 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 997 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 998 ierr = VecDestroy(b);CHKERRQ(ierr); 999 break; 1000 } 1001 #if defined(PETSC_USE_COMPLEX) 1002 zmumps_c(&lu->id); 1003 #else 1004 dmumps_c(&lu->id); 1005 #endif 1006 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)); 1007 1008 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1009 F->ops->solve = MatSolve_MUMPS; 1010 #if !defined(PETSC_USE_COMPLEX) 1011 (F)->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1012 #endif 1013 PetscFunctionReturn(0); 1014 } 1015 1016 #undef __FUNCT__ 1017 #define __FUNCT__ "MatView_MUMPS" 1018 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1019 { 1020 PetscErrorCode ierr; 1021 PetscBool iascii; 1022 PetscViewerFormat format; 1023 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 1024 1025 PetscFunctionBegin; 1026 /* check if matrix is mumps type */ 1027 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1028 1029 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1030 if (iascii) { 1031 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1032 if (format == PETSC_VIEWER_ASCII_INFO){ 1033 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1034 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 1035 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 1036 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 1037 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 1038 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 1039 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 1040 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 1041 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 1042 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 1043 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 1044 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 1045 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 1046 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 1047 if (lu->id.ICNTL(11)>0) { 1048 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 1049 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 1050 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 1051 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",lu->id.RINFOG(7),lu->id.RINFOG(8));CHKERRQ(ierr); 1052 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 1053 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 1054 } 1055 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 1056 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 1057 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 1058 /* ICNTL(15-17) not used */ 1059 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 1060 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 1061 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 1062 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 1063 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 1064 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 1065 1066 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 1067 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 1068 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 1069 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 1070 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (Use parallel or sequential ordering): %d \n",lu->id.ICNTL(28));CHKERRQ(ierr); 1071 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (Parallel ordering): %d \n",lu->id.ICNTL(29));CHKERRQ(ierr); 1072 1073 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 1074 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 1075 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 1076 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 1077 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 1078 1079 /* infomation local to each processor */ 1080 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1081 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 1082 ierr = PetscViewerFlush(viewer); 1083 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1084 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 1085 ierr = PetscViewerFlush(viewer); 1086 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1087 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 1088 ierr = PetscViewerFlush(viewer); 1089 1090 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1091 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 1092 ierr = PetscViewerFlush(viewer); 1093 1094 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1095 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 1096 ierr = PetscViewerFlush(viewer); 1097 1098 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1099 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 1100 ierr = PetscViewerFlush(viewer); 1101 1102 if (!lu->myid){ /* information from the host */ 1103 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 1104 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 1105 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 1106 1107 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 1108 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 1109 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 1110 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 1111 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 1112 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 1113 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 1114 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 1115 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 1116 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 1117 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 1118 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 1119 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 1120 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); 1121 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); 1122 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); 1123 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); 1124 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 1125 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); 1126 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); 1127 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 1128 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 1129 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 1130 } 1131 } 1132 } 1133 PetscFunctionReturn(0); 1134 } 1135 1136 #undef __FUNCT__ 1137 #define __FUNCT__ "MatGetInfo_MUMPS" 1138 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1139 { 1140 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1141 1142 PetscFunctionBegin; 1143 info->block_size = 1.0; 1144 info->nz_allocated = mumps->id.INFOG(20); 1145 info->nz_used = mumps->id.INFOG(20); 1146 info->nz_unneeded = 0.0; 1147 info->assemblies = 0.0; 1148 info->mallocs = 0.0; 1149 info->memory = 0.0; 1150 info->fill_ratio_given = 0; 1151 info->fill_ratio_needed = 0; 1152 info->factor_mallocs = 0; 1153 PetscFunctionReturn(0); 1154 } 1155 1156 /* -------------------------------------------------------------------------------------------*/ 1157 #undef __FUNCT__ 1158 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1159 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1160 { 1161 Mat_MUMPS *lu =(Mat_MUMPS*)F->spptr; 1162 1163 PetscFunctionBegin; 1164 lu->id.ICNTL(icntl) = ival; 1165 PetscFunctionReturn(0); 1166 } 1167 1168 #undef __FUNCT__ 1169 #define __FUNCT__ "MatMumpsSetIcntl" 1170 /*@ 1171 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1172 1173 Logically Collective on Mat 1174 1175 Input Parameters: 1176 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1177 . icntl - index of MUMPS parameter array ICNTL() 1178 - ival - value of MUMPS ICNTL(icntl) 1179 1180 Options Database: 1181 . -mat_mumps_icntl_<icntl> <ival> 1182 1183 Level: beginner 1184 1185 References: MUMPS Users' Guide 1186 1187 .seealso: MatGetFactor() 1188 @*/ 1189 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1190 { 1191 PetscErrorCode ierr; 1192 1193 PetscFunctionBegin; 1194 PetscValidLogicalCollectiveInt(F,icntl,2); 1195 PetscValidLogicalCollectiveInt(F,ival,3); 1196 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1197 PetscFunctionReturn(0); 1198 } 1199 1200 /*MC 1201 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1202 distributed and sequential matrices via the external package MUMPS. 1203 1204 Works with MATAIJ and MATSBAIJ matrices 1205 1206 Options Database Keys: 1207 + -mat_mumps_icntl_4 <0,...,4> - print level 1208 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1209 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec) 1210 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1211 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1212 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1213 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1214 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1215 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1216 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1217 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1218 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1219 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1220 1221 Level: beginner 1222 1223 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1224 1225 M*/ 1226 1227 EXTERN_C_BEGIN 1228 #undef __FUNCT__ 1229 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1230 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1231 { 1232 PetscFunctionBegin; 1233 *type = MATSOLVERMUMPS; 1234 PetscFunctionReturn(0); 1235 } 1236 EXTERN_C_END 1237 1238 EXTERN_C_BEGIN 1239 /* MatGetFactor for Seq and MPI AIJ matrices */ 1240 #undef __FUNCT__ 1241 #define __FUNCT__ "MatGetFactor_aij_mumps" 1242 PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1243 { 1244 Mat B; 1245 PetscErrorCode ierr; 1246 Mat_MUMPS *mumps; 1247 PetscBool isSeqAIJ; 1248 1249 PetscFunctionBegin; 1250 /* Create the factorization matrix */ 1251 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1252 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1253 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1254 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1255 if (isSeqAIJ) { 1256 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1257 } else { 1258 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1259 } 1260 1261 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1262 B->ops->view = MatView_MUMPS; 1263 B->ops->getinfo = MatGetInfo_MUMPS; 1264 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1265 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl_MUMPS",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1266 if (ftype == MAT_FACTOR_LU) { 1267 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1268 B->factortype = MAT_FACTOR_LU; 1269 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1270 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1271 mumps->sym = 0; 1272 } else { 1273 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1274 B->factortype = MAT_FACTOR_CHOLESKY; 1275 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1276 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1277 if (A->spd_set && A->spd) mumps->sym = 1; 1278 else mumps->sym = 2; 1279 } 1280 1281 mumps->isAIJ = PETSC_TRUE; 1282 mumps->MatDestroy = B->ops->destroy; 1283 B->ops->destroy = MatDestroy_MUMPS; 1284 B->spptr = (void*)mumps; 1285 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1286 1287 *F = B; 1288 PetscFunctionReturn(0); 1289 } 1290 EXTERN_C_END 1291 1292 1293 EXTERN_C_BEGIN 1294 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1297 PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1298 { 1299 Mat B; 1300 PetscErrorCode ierr; 1301 Mat_MUMPS *mumps; 1302 PetscBool isSeqSBAIJ; 1303 1304 PetscFunctionBegin; 1305 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1306 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"); 1307 ierr = PetscTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1308 /* Create the factorization matrix */ 1309 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1310 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1311 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1312 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1313 if (isSeqSBAIJ) { 1314 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1315 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1316 } else { 1317 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1318 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1319 } 1320 1321 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1322 B->ops->view = MatView_MUMPS; 1323 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1324 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1325 B->factortype = MAT_FACTOR_CHOLESKY; 1326 if (A->spd_set && A->spd) mumps->sym = 1; 1327 else mumps->sym = 2; 1328 1329 mumps->isAIJ = PETSC_FALSE; 1330 mumps->MatDestroy = B->ops->destroy; 1331 B->ops->destroy = MatDestroy_MUMPS; 1332 B->spptr = (void*)mumps; 1333 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1334 1335 *F = B; 1336 PetscFunctionReturn(0); 1337 } 1338 EXTERN_C_END 1339 1340 EXTERN_C_BEGIN 1341 #undef __FUNCT__ 1342 #define __FUNCT__ "MatGetFactor_baij_mumps" 1343 PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 1344 { 1345 Mat B; 1346 PetscErrorCode ierr; 1347 Mat_MUMPS *mumps; 1348 PetscBool isSeqBAIJ; 1349 1350 PetscFunctionBegin; 1351 /* Create the factorization matrix */ 1352 ierr = PetscTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1353 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1354 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1355 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1356 if (isSeqBAIJ) { 1357 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,PETSC_NULL);CHKERRQ(ierr); 1358 } else { 1359 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1360 } 1361 1362 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1363 if (ftype == MAT_FACTOR_LU) { 1364 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1365 B->factortype = MAT_FACTOR_LU; 1366 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 1367 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 1368 mumps->sym = 0; 1369 } else { 1370 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 1371 } 1372 1373 B->ops->view = MatView_MUMPS; 1374 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1375 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl_MUMPS",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1376 1377 mumps->isAIJ = PETSC_TRUE; 1378 mumps->MatDestroy = B->ops->destroy; 1379 B->ops->destroy = MatDestroy_MUMPS; 1380 B->spptr = (void*)mumps; 1381 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1382 1383 *F = B; 1384 PetscFunctionReturn(0); 1385 } 1386 EXTERN_C_END 1387 1388