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