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