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