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