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