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 11 EXTERN_C_BEGIN 12 #if defined(PETSC_USE_COMPLEX) 13 #include "zmumps_c.h" 14 #else 15 #include "dmumps_c.h" 16 #endif 17 EXTERN_C_END 18 #define JOB_INIT -1 19 #define JOB_END -2 20 /* macros s.t. indices match MUMPS documentation */ 21 #define ICNTL(I) icntl[(I)-1] 22 #define CNTL(I) cntl[(I)-1] 23 #define INFOG(I) infog[(I)-1] 24 #define INFO(I) info[(I)-1] 25 #define RINFOG(I) rinfog[(I)-1] 26 #define RINFO(I) rinfo[(I)-1] 27 28 typedef struct { 29 #if defined(PETSC_USE_COMPLEX) 30 ZMUMPS_STRUC_C id; 31 #else 32 DMUMPS_STRUC_C id; 33 #endif 34 MatStructure matstruc; 35 PetscMPIInt myid,size; 36 PetscInt *irn,*jcn,sym,nSolve; 37 PetscScalar *val; 38 MPI_Comm comm_mumps; 39 VecScatter scat_rhs, scat_sol; 40 PetscTruth isAIJ,CleanUpMUMPS; 41 Vec b_seq,x_seq; 42 PetscErrorCode (*MatDestroy)(Mat); 43 } Mat_MUMPS; 44 45 EXTERN PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*); 46 47 /* convert Petsc mpiaij matrix to triples: row[nz], col[nz], val[nz] */ 48 /* 49 input: 50 A - matrix in mpiaij or mpisbaij (bs=1) format 51 shift - 0: C style output triple; 1: Fortran style output triple. 52 valOnly - FALSE: spaces are allocated and values are set for the triple 53 TRUE: only the values in v array are updated 54 output: 55 nnz - dim of r, c, and v (number of local nonzero entries of A) 56 r, c, v - row and col index, matrix values (matrix triples) 57 */ 58 PetscErrorCode MatConvertToTriples(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 59 { 60 PetscInt *ai, *aj, *bi, *bj, rstart,nz, *garray; 61 PetscErrorCode ierr; 62 PetscInt i,j,jj,jB,irow,m=A->rmap->n,*ajj,*bjj,countA,countB,colA_start,jcol; 63 PetscInt *row,*col; 64 PetscScalar *av, *bv,*val; 65 PetscTruth isAIJ; 66 67 PetscFunctionBegin; 68 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isAIJ);CHKERRQ(ierr); 69 if (isAIJ){ 70 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 71 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 72 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 73 nz = aa->nz + bb->nz; 74 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 75 garray = mat->garray; 76 av=aa->a; bv=bb->a; 77 78 } else { 79 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 80 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)(mat->A)->data; 81 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 82 if (A->rmap->bs > 1) SETERRQ1(PETSC_ERR_SUP," bs=%d is not supported yet\n", A->rmap->bs); 83 nz = aa->nz + bb->nz; 84 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 85 garray = mat->garray; 86 av=aa->a; bv=bb->a; 87 } 88 89 if (!valOnly){ 90 ierr = PetscMalloc(nz*sizeof(PetscInt) ,&row);CHKERRQ(ierr); 91 ierr = PetscMalloc(nz*sizeof(PetscInt),&col);CHKERRQ(ierr); 92 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 93 *r = row; *c = col; *v = val; 94 } else { 95 row = *r; col = *c; val = *v; 96 } 97 *nnz = nz; 98 99 jj = 0; irow = rstart; 100 for ( i=0; i<m; i++ ) { 101 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 102 countA = ai[i+1] - ai[i]; 103 countB = bi[i+1] - bi[i]; 104 bjj = bj + bi[i]; 105 106 /* get jB, the starting local col index for the 2nd B-part */ 107 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 108 j=-1; 109 do { 110 j++; 111 if (j == countB) break; 112 jcol = garray[bjj[j]]; 113 } while (jcol < colA_start); 114 jB = j; 115 116 /* B-part, smaller col index */ 117 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 118 for (j=0; j<jB; j++){ 119 jcol = garray[bjj[j]]; 120 if (!valOnly){ 121 row[jj] = irow + shift; col[jj] = jcol + shift; 122 123 } 124 val[jj++] = *bv++; 125 } 126 /* A-part */ 127 for (j=0; j<countA; j++){ 128 if (!valOnly){ 129 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 130 } 131 val[jj++] = *av++; 132 } 133 /* B-part, larger col index */ 134 for (j=jB; j<countB; j++){ 135 if (!valOnly){ 136 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 137 } 138 val[jj++] = *bv++; 139 } 140 irow++; 141 } 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "MatDestroy_MUMPS" 147 PetscErrorCode MatDestroy_MUMPS(Mat A) 148 { 149 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 150 PetscErrorCode ierr; 151 PetscMPIInt size=lu->size; 152 153 PetscFunctionBegin; 154 if (lu->CleanUpMUMPS) { 155 /* Terminate instance, deallocate memories */ 156 if (size > 1){ 157 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 158 ierr = VecScatterDestroy(lu->scat_rhs);CHKERRQ(ierr); 159 ierr = VecDestroy(lu->b_seq);CHKERRQ(ierr); 160 if (lu->nSolve && lu->scat_sol){ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);} 161 if (lu->nSolve && lu->x_seq){ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);} 162 ierr = PetscFree(lu->val);CHKERRQ(ierr); 163 } 164 lu->id.job=JOB_END; 165 #if defined(PETSC_USE_COMPLEX) 166 zmumps_c(&lu->id); 167 #else 168 dmumps_c(&lu->id); 169 #endif 170 ierr = PetscFree(lu->irn);CHKERRQ(ierr); 171 ierr = PetscFree(lu->jcn);CHKERRQ(ierr); 172 ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr); 173 } 174 /* clear composed functions */ 175 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatFactorGetSolverPackage_C","",PETSC_NULL);CHKERRQ(ierr); 176 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMumpsSetIcntl_C","",PETSC_NULL);CHKERRQ(ierr); 177 ierr = (lu->MatDestroy)(A);CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNCT__ 182 #define __FUNCT__ "MatSolve_MUMPS" 183 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 184 { 185 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 186 PetscScalar *array; 187 Vec x_seq; 188 IS is_iden,is_petsc; 189 PetscErrorCode ierr; 190 PetscInt i; 191 192 PetscFunctionBegin; 193 lu->id.nrhs = 1; 194 x_seq = lu->b_seq; 195 if (lu->size > 1){ 196 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 197 ierr = VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 198 ierr = VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 199 if (!lu->myid) {ierr = VecGetArray(x_seq,&array);CHKERRQ(ierr);} 200 } else { /* size == 1 */ 201 ierr = VecCopy(b,x);CHKERRQ(ierr); 202 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 203 } 204 if (!lu->myid) { /* define rhs on the host */ 205 lu->id.nrhs = 1; 206 #if defined(PETSC_USE_COMPLEX) 207 lu->id.rhs = (mumps_double_complex*)array; 208 #else 209 lu->id.rhs = array; 210 #endif 211 } 212 if (lu->size == 1){ 213 ierr = VecRestoreArray(x,&array);CHKERRQ(ierr); 214 } else if (!lu->myid){ 215 ierr = VecRestoreArray(x_seq,&array);CHKERRQ(ierr); 216 } 217 218 if (lu->size > 1){ 219 /* distributed solution */ 220 lu->id.ICNTL(21) = 1; 221 if (!lu->nSolve){ 222 /* Create x_seq=sol_loc for repeated use */ 223 PetscInt lsol_loc; 224 PetscScalar *sol_loc; 225 lsol_loc = lu->id.INFO(23); /* length of sol_loc */ 226 ierr = PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&lu->id.isol_loc);CHKERRQ(ierr); 227 lu->id.lsol_loc = lsol_loc; 228 #if defined(PETSC_USE_COMPLEX) 229 lu->id.sol_loc = (mumps_double_complex*)sol_loc; 230 #else 231 lu->id.sol_loc = sol_loc; 232 #endif 233 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,lsol_loc,sol_loc,&lu->x_seq);CHKERRQ(ierr); 234 } 235 } 236 237 /* solve phase */ 238 /*-------------*/ 239 lu->id.job = 3; 240 #if defined(PETSC_USE_COMPLEX) 241 zmumps_c(&lu->id); 242 #else 243 dmumps_c(&lu->id); 244 #endif 245 if (lu->id.INFOG(1) < 0) { 246 SETERRQ1(PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 247 } 248 249 if (lu->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 250 if (!lu->nSolve){ /* create scatter scat_sol */ 251 ierr = ISCreateStride(PETSC_COMM_SELF,lu->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 252 for (i=0; i<lu->id.lsol_loc; i++){ 253 lu->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 254 } 255 ierr = ISCreateGeneral(PETSC_COMM_SELF,lu->id.lsol_loc,lu->id.isol_loc,&is_petsc);CHKERRQ(ierr); /* to */ 256 ierr = VecScatterCreate(lu->x_seq,is_iden,x,is_petsc,&lu->scat_sol);CHKERRQ(ierr); 257 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 258 ierr = ISDestroy(is_petsc);CHKERRQ(ierr); 259 } 260 ierr = VecScatterBegin(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 261 ierr = VecScatterEnd(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 262 } 263 lu->nSolve++; 264 PetscFunctionReturn(0); 265 } 266 267 #if !defined(PETSC_USE_COMPLEX) 268 /* 269 input: 270 F: numeric factor 271 output: 272 nneg: total number of negative pivots 273 nzero: 0 274 npos: (global dimension of F) - nneg 275 */ 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 279 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 280 { 281 Mat_MUMPS *lu =(Mat_MUMPS*)F->spptr; 282 PetscErrorCode ierr; 283 PetscMPIInt size; 284 285 PetscFunctionBegin; 286 ierr = MPI_Comm_size(((PetscObject)F)->comm,&size);CHKERRQ(ierr); 287 /* 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 */ 288 if (size > 1 && lu->id.ICNTL(13) != 1){ 289 SETERRQ1(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)); 290 } 291 if (nneg){ 292 if (!lu->myid){ 293 *nneg = lu->id.INFOG(12); 294 } 295 ierr = MPI_Bcast(nneg,1,MPI_INT,0,lu->comm_mumps);CHKERRQ(ierr); 296 } 297 if (nzero) *nzero = 0; 298 if (npos) *npos = F->rmap->N - (*nneg); 299 PetscFunctionReturn(0); 300 } 301 #endif /* !defined(PETSC_USE_COMPLEX) */ 302 303 #undef __FUNCT__ 304 #define __FUNCT__ "MatFactorNumeric_MUMPS" 305 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 306 { 307 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 308 PetscErrorCode ierr; 309 PetscInt rnz,nnz,nz=0,i,M=A->rmap->N,*ai,*aj,icntl; 310 PetscTruth valOnly,flg; 311 Mat F_diag; 312 IS is_iden; 313 Vec b; 314 PetscTruth isSeqAIJ,isSeqSBAIJ; 315 316 PetscFunctionBegin; 317 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 318 ierr = PetscTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 319 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 320 (F)->ops->solve = MatSolve_MUMPS; 321 322 /* Initialize a MUMPS instance */ 323 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 324 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 325 lu->id.job = JOB_INIT; 326 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 327 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 328 329 /* Set mumps options */ 330 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 331 lu->id.par=1; /* host participates factorizaton and solve */ 332 lu->id.sym=lu->sym; 333 if (lu->sym == 2){ 334 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 335 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 336 } 337 #if defined(PETSC_USE_COMPLEX) 338 zmumps_c(&lu->id); 339 #else 340 dmumps_c(&lu->id); 341 #endif 342 343 if (isSeqAIJ || isSeqSBAIJ){ 344 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 345 } else { 346 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 347 } 348 349 icntl=-1; 350 lu->id.ICNTL(4) = 0; /* level of printing; overwrite mumps default ICNTL(4)=2 */ 351 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 352 if ((flg && icntl > 0) || PetscLogPrintInfo) { 353 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 354 } else { /* no output */ 355 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 356 lu->id.ICNTL(2) = 0; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 357 lu->id.ICNTL(3) = 0; /* output stream for global information, default=6 */ 358 } 359 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); 360 icntl=-1; 361 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 362 if (flg) { 363 if (icntl== 1){ 364 SETERRQ(PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n"); 365 } else { 366 lu->id.ICNTL(7) = icntl; 367 } 368 } 369 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); 370 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); 371 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); 372 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); 373 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); 374 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); 375 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); 376 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",lu->id.ICNTL(19),&lu->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr); 377 378 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); 379 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); 380 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); 381 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); 382 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); 383 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",lu->id.ICNTL(27),&lu->id.ICNTL(27),PETSC_NULL);CHKERRQ(ierr); 384 385 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 386 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); 387 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 388 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); 389 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); 390 PetscOptionsEnd(); 391 } 392 393 /* define matrix A */ 394 switch (lu->id.ICNTL(18)){ 395 case 0: /* centralized assembled matrix input (size=1) */ 396 if (!lu->myid) { 397 if (isSeqAIJ){ 398 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 399 nz = aa->nz; 400 ai = aa->i; aj = aa->j; lu->val = aa->a; 401 } else if (isSeqSBAIJ) { 402 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 403 nz = aa->nz; 404 ai = aa->i; aj = aa->j; lu->val = aa->a; 405 } else { 406 SETERRQ(PETSC_ERR_SUP,"No mumps factorization for this matrix type"); 407 } 408 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization, get irn and jcn */ 409 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr); 410 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr); 411 nz = 0; 412 for (i=0; i<M; i++){ 413 rnz = ai[i+1] - ai[i]; 414 while (rnz--) { /* Fortran row/col index! */ 415 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 416 } 417 } 418 } 419 } 420 break; 421 case 3: /* distributed assembled matrix input (size>1) */ 422 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 423 valOnly = PETSC_FALSE; 424 } else { 425 valOnly = PETSC_TRUE; /* only update mat values, not row and col index */ 426 } 427 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 428 break; 429 default: SETERRQ(PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 430 } 431 432 /* analysis phase */ 433 /*----------------*/ 434 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 435 lu->id.job = 1; 436 437 lu->id.n = M; 438 switch (lu->id.ICNTL(18)){ 439 case 0: /* centralized assembled matrix input */ 440 if (!lu->myid) { 441 lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 442 if (lu->id.ICNTL(6)>1){ 443 #if defined(PETSC_USE_COMPLEX) 444 lu->id.a = (mumps_double_complex*)lu->val; 445 #else 446 lu->id.a = lu->val; 447 #endif 448 } 449 } 450 break; 451 case 3: /* distributed assembled matrix input (size>1) */ 452 lu->id.nz_loc = nnz; 453 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 454 if (lu->id.ICNTL(6)>1) { 455 #if defined(PETSC_USE_COMPLEX) 456 lu->id.a_loc = (mumps_double_complex*)lu->val; 457 #else 458 lu->id.a_loc = lu->val; 459 #endif 460 } 461 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 462 if (!lu->myid){ 463 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 464 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 465 } else { 466 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 467 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 468 } 469 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 470 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 471 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 472 473 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 474 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 475 ierr = VecDestroy(b);CHKERRQ(ierr); 476 break; 477 } 478 #if defined(PETSC_USE_COMPLEX) 479 zmumps_c(&lu->id); 480 #else 481 dmumps_c(&lu->id); 482 #endif 483 if (lu->id.INFOG(1) < 0) { 484 SETERRQ1(PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 485 } 486 } 487 488 /* numerical factorization phase */ 489 /*-------------------------------*/ 490 lu->id.job = 2; 491 if(!lu->id.ICNTL(18)) { 492 if (!lu->myid) { 493 #if defined(PETSC_USE_COMPLEX) 494 lu->id.a = (mumps_double_complex*)lu->val; 495 #else 496 lu->id.a = lu->val; 497 #endif 498 } 499 } else { 500 #if defined(PETSC_USE_COMPLEX) 501 lu->id.a_loc = (mumps_double_complex*)lu->val; 502 #else 503 lu->id.a_loc = lu->val; 504 #endif 505 } 506 #if defined(PETSC_USE_COMPLEX) 507 zmumps_c(&lu->id); 508 #else 509 dmumps_c(&lu->id); 510 #endif 511 if (lu->id.INFOG(1) < 0) { 512 if (lu->id.INFO(1) == -13) { 513 SETERRQ1(PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d megabytes\n",lu->id.INFO(2)); 514 } else { 515 SETERRQ2(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)); 516 } 517 } 518 519 if (!lu->myid && lu->id.ICNTL(16) > 0){ 520 SETERRQ1(PETSC_ERR_LIB," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 521 } 522 523 if (lu->size > 1){ 524 if ((F)->factor == MAT_FACTOR_LU){ 525 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 526 } else { 527 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 528 } 529 F_diag->assembled = PETSC_TRUE; 530 if (lu->nSolve){ 531 ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr); 532 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 533 ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr); 534 } 535 } 536 (F)->assembled = PETSC_TRUE; 537 lu->matstruc = SAME_NONZERO_PATTERN; 538 lu->CleanUpMUMPS = PETSC_TRUE; 539 lu->nSolve = 0; 540 PetscFunctionReturn(0); 541 } 542 543 /* Note the Petsc r and c permutations are ignored */ 544 #undef __FUNCT__ 545 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 546 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 547 { 548 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 549 550 PetscFunctionBegin; 551 lu->sym = 0; 552 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 553 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 554 PetscFunctionReturn(0); 555 } 556 557 558 /* Note the Petsc r permutation is ignored */ 559 #undef __FUNCT__ 560 #define __FUNCT__ "MatCholeskyFactorSymbolic_SBAIJMUMPS" 561 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJMUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 562 { 563 Mat_MUMPS *lu = (Mat_MUMPS*)(F)->spptr; 564 565 PetscFunctionBegin; 566 lu->sym = 2; 567 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 568 (F)->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 569 #if !defined(PETSC_USE_COMPLEX) 570 (F)->ops->getinertia = MatGetInertia_SBAIJMUMPS; 571 #endif 572 PetscFunctionReturn(0); 573 } 574 575 #undef __FUNCT__ 576 #define __FUNCT__ "MatFactorInfo_MUMPS" 577 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 578 { 579 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 580 PetscErrorCode ierr; 581 582 PetscFunctionBegin; 583 /* check if matrix is mumps type */ 584 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 585 586 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 587 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 588 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 589 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 590 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 591 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 592 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 593 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 594 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 595 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 596 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 597 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 598 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 599 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 600 if (!lu->myid && lu->id.ICNTL(11)>0) { 601 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 602 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 603 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 604 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); 605 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 606 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 607 608 } 609 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 610 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 611 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 612 /* ICNTL(15-17) not used */ 613 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 614 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 615 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 616 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 617 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 618 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 619 620 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 621 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 622 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 623 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 624 625 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 626 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 627 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 628 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 629 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 630 631 /* infomation local to each processor */ 632 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr);} 633 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 634 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 635 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr);} 636 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 637 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 638 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr);} 639 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 640 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 641 642 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);} 643 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 644 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 645 646 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr);} 647 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 648 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 649 650 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr);} 651 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 652 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 653 654 if (!lu->myid){ /* information from the host */ 655 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 656 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 657 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 658 659 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 660 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 661 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 662 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 663 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 664 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 665 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 666 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 667 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 668 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 669 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 670 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 671 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 672 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); 673 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); 674 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); 675 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); 676 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 677 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); 678 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); 679 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 680 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 681 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 682 } 683 PetscFunctionReturn(0); 684 } 685 686 #undef __FUNCT__ 687 #define __FUNCT__ "MatView_MUMPS" 688 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 689 { 690 PetscErrorCode ierr; 691 PetscTruth iascii; 692 PetscViewerFormat format; 693 694 PetscFunctionBegin; 695 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 696 if (iascii) { 697 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 698 if (format == PETSC_VIEWER_ASCII_INFO){ 699 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 700 } 701 } 702 PetscFunctionReturn(0); 703 } 704 705 #undef __FUNCT__ 706 #define __FUNCT__ "MatGetInfo_MUMPS" 707 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 708 { 709 Mat_MUMPS *lu =(Mat_MUMPS*)A->spptr; 710 711 PetscFunctionBegin; 712 info->block_size = 1.0; 713 info->nz_allocated = lu->id.INFOG(20); 714 info->nz_used = lu->id.INFOG(20); 715 info->nz_unneeded = 0.0; 716 info->assemblies = 0.0; 717 info->mallocs = 0.0; 718 info->memory = 0.0; 719 info->fill_ratio_given = 0; 720 info->fill_ratio_needed = 0; 721 info->factor_mallocs = 0; 722 PetscFunctionReturn(0); 723 } 724 725 /*MC 726 MAT_SOLVER_MUMPS - A matrix type providing direct solvers (LU and Cholesky) for 727 distributed and sequential matrices via the external package MUMPS. 728 729 Works with MATAIJ and MATSBAIJ matrices 730 731 Options Database Keys: 732 + -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 733 . -mat_mumps_icntl_4 <0,...,4> - print level 734 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 735 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 736 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 737 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 738 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 739 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 740 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 741 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 742 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 743 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 744 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 745 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 746 747 Level: beginner 748 749 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 750 751 M*/ 752 753 EXTERN_C_BEGIN 754 #undef __FUNCT__ 755 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 756 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 757 { 758 PetscFunctionBegin; 759 *type = MAT_SOLVER_MUMPS; 760 PetscFunctionReturn(0); 761 } 762 EXTERN_C_END 763 764 EXTERN_C_BEGIN 765 /* 766 The seq and mpi versions of this function are the same 767 */ 768 #undef __FUNCT__ 769 #define __FUNCT__ "MatGetFactor_seqaij_mumps" 770 PetscErrorCode MatGetFactor_seqaij_mumps(Mat A,MatFactorType ftype,Mat *F) 771 { 772 Mat B; 773 PetscErrorCode ierr; 774 Mat_MUMPS *mumps; 775 776 PetscFunctionBegin; 777 if (ftype != MAT_FACTOR_LU) { 778 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with MUMPS Cholesky, use SBAIJ matrix"); 779 } 780 /* Create the factorization matrix */ 781 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 782 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 783 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 784 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 785 786 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 787 B->ops->view = MatView_MUMPS; 788 B->ops->getinfo = MatGetInfo_MUMPS; 789 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 790 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 791 B->factor = MAT_FACTOR_LU; 792 793 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 794 mumps->CleanUpMUMPS = PETSC_FALSE; 795 mumps->isAIJ = PETSC_TRUE; 796 mumps->scat_rhs = PETSC_NULL; 797 mumps->scat_sol = PETSC_NULL; 798 mumps->nSolve = 0; 799 mumps->MatDestroy = B->ops->destroy; 800 B->ops->destroy = MatDestroy_MUMPS; 801 B->spptr = (void*)mumps; 802 803 *F = B; 804 PetscFunctionReturn(0); 805 } 806 EXTERN_C_END 807 808 EXTERN_C_BEGIN 809 #undef __FUNCT__ 810 #define __FUNCT__ "MatGetFactor_mpiaij_mumps" 811 PetscErrorCode MatGetFactor_mpiaij_mumps(Mat A,MatFactorType ftype,Mat *F) 812 { 813 Mat B; 814 PetscErrorCode ierr; 815 Mat_MUMPS *mumps; 816 817 PetscFunctionBegin; 818 if (ftype != MAT_FACTOR_LU) { 819 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc AIJ matrices with MUMPS Cholesky, use SBAIJ matrix"); 820 } 821 /* Create the factorization matrix */ 822 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 823 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 824 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 825 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 826 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 827 828 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 829 B->ops->view = MatView_MUMPS; 830 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 831 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 832 B->factor = MAT_FACTOR_LU; 833 834 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 835 mumps->CleanUpMUMPS = PETSC_FALSE; 836 mumps->isAIJ = PETSC_TRUE; 837 mumps->scat_rhs = PETSC_NULL; 838 mumps->scat_sol = PETSC_NULL; 839 mumps->nSolve = 0; 840 mumps->MatDestroy = B->ops->destroy; 841 B->ops->destroy = MatDestroy_MUMPS; 842 B->spptr = (void*)mumps; 843 844 *F = B; 845 PetscFunctionReturn(0); 846 } 847 EXTERN_C_END 848 849 EXTERN_C_BEGIN 850 #undef __FUNCT__ 851 #define __FUNCT__ "MatGetFactor_seqsbaij_mumps" 852 PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 853 { 854 Mat B; 855 PetscErrorCode ierr; 856 Mat_MUMPS *mumps; 857 858 PetscFunctionBegin; 859 if (ftype != MAT_FACTOR_CHOLESKY) { 860 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 861 } 862 /* Create the factorization matrix */ 863 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 864 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 865 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 866 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 867 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 868 869 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 870 B->ops->view = MatView_MUMPS; 871 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 872 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 873 B->factor = MAT_FACTOR_CHOLESKY; 874 875 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 876 mumps->CleanUpMUMPS = PETSC_FALSE; 877 mumps->isAIJ = PETSC_TRUE; 878 mumps->scat_rhs = PETSC_NULL; 879 mumps->scat_sol = PETSC_NULL; 880 mumps->nSolve = 0; 881 mumps->MatDestroy = B->ops->destroy; 882 B->ops->destroy = MatDestroy_MUMPS; 883 B->spptr = (void*)mumps; 884 885 *F = B; 886 PetscFunctionReturn(0); 887 } 888 EXTERN_C_END 889 890 EXTERN_C_BEGIN 891 #undef __FUNCT__ 892 #define __FUNCT__ "MatGetFactor_mpisbaij_mumps" 893 PetscErrorCode MatGetFactor_mpisbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 894 { 895 Mat B; 896 PetscErrorCode ierr; 897 Mat_MUMPS *mumps; 898 899 PetscFunctionBegin; 900 if (ftype != MAT_FACTOR_CHOLESKY) { 901 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 902 } 903 /* Create the factorization matrix */ 904 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 905 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 906 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 907 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 908 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 909 910 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 911 B->ops->view = MatView_MUMPS; 912 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 913 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 914 B->factor = MAT_FACTOR_CHOLESKY; 915 916 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 917 mumps->CleanUpMUMPS = PETSC_FALSE; 918 mumps->isAIJ = PETSC_TRUE; 919 mumps->scat_rhs = PETSC_NULL; 920 mumps->scat_sol = PETSC_NULL; 921 mumps->nSolve = 0; 922 mumps->MatDestroy = B->ops->destroy; 923 B->ops->destroy = MatDestroy_MUMPS; 924 B->spptr = (void*)mumps; 925 926 *F = B; 927 PetscFunctionReturn(0); 928 } 929 EXTERN_C_END 930 931 /* -------------------------------------------------------------------------------------------*/ 932 /*@ 933 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 934 935 Collective on Mat 936 937 Input Parameters: 938 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 939 . idx - index of MUMPS parameter array ICNTL() 940 - icntl - value of MUMPS ICNTL(imumps) 941 942 Options Database: 943 . -mat_mumps_icntl_<idx> <icntl> 944 945 Level: beginner 946 947 References: MUMPS Users' Guide 948 949 .seealso: MatGetFactor() 950 @*/ 951 #undef __FUNCT__ 952 #define __FUNCT__ "MatMumpsSetIcntl" 953 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt idx,PetscInt icntl) 954 { 955 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 956 957 PetscFunctionBegin; 958 lu->id.ICNTL(idx) = icntl; 959 PetscFunctionReturn(0); 960 } 961 962