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