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