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,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 ((F)->factortype == MAT_FACTOR_LU){ */ 559 if(isMPIAIJ) { 560 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 561 } else { 562 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 563 } 564 F_diag->assembled = PETSC_TRUE; 565 if (lu->nSolve){ 566 ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr); 567 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 568 ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr); 569 } 570 } 571 (F)->assembled = PETSC_TRUE; 572 lu->matstruc = SAME_NONZERO_PATTERN; 573 lu->CleanUpMUMPS = PETSC_TRUE; 574 lu->nSolve = 0; 575 PetscFunctionReturn(0); 576 } 577 578 /* Note the Petsc r and c permutations are ignored */ 579 #undef __FUNCT__ 580 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 581 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 582 { 583 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 584 585 PetscFunctionBegin; 586 lu->sym = 0; 587 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 588 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 589 PetscFunctionReturn(0); 590 } 591 592 /* Note the Petsc r and c permutations are ignored */ 593 #undef __FUNCT__ 594 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 595 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 596 { 597 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 598 599 PetscFunctionBegin; 600 lu->sym = 0; 601 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 602 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 603 PetscFunctionReturn(0); 604 } 605 606 /* Note the Petsc r permutation is ignored */ 607 #undef __FUNCT__ 608 #define __FUNCT__ "MatCholeskyFactorSymbolic_SBAIJMUMPS" 609 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJMUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 610 { 611 Mat_MUMPS *lu = (Mat_MUMPS*)(F)->spptr; 612 613 PetscFunctionBegin; 614 lu->sym = 2; 615 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 616 (F)->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 617 #if !defined(PETSC_USE_COMPLEX) 618 (F)->ops->getinertia = MatGetInertia_SBAIJMUMPS; 619 #endif 620 PetscFunctionReturn(0); 621 } 622 623 #undef __FUNCT__ 624 #define __FUNCT__ "MatFactorInfo_MUMPS" 625 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 626 { 627 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 628 PetscErrorCode ierr; 629 630 PetscFunctionBegin; 631 /* check if matrix is mumps type */ 632 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 633 634 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 635 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 636 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 637 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 638 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 639 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 640 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 641 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 642 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 643 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 644 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 645 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 646 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 647 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 648 if (!lu->myid && lu->id.ICNTL(11)>0) { 649 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 650 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 651 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 652 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); 653 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 654 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 655 656 } 657 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 658 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 659 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 660 /* ICNTL(15-17) not used */ 661 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 662 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 663 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 664 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 665 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 666 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 667 668 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 669 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 670 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 671 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 672 673 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 674 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 675 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 676 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 677 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 678 679 /* infomation local to each processor */ 680 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr);} 681 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 682 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 683 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr);} 684 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 685 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 686 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr);} 687 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 688 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 689 690 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);} 691 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 692 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 693 694 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr);} 695 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 696 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 697 698 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr);} 699 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 700 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 701 702 if (!lu->myid){ /* information from the host */ 703 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 704 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 705 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 706 707 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 708 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 709 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 710 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 711 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 712 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 713 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 714 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 715 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 716 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 717 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 718 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 719 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 720 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); 721 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); 722 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); 723 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); 724 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 725 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); 726 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); 727 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 728 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 729 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 730 } 731 PetscFunctionReturn(0); 732 } 733 734 #undef __FUNCT__ 735 #define __FUNCT__ "MatView_MUMPS" 736 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 737 { 738 PetscErrorCode ierr; 739 PetscTruth iascii; 740 PetscViewerFormat format; 741 742 PetscFunctionBegin; 743 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 744 if (iascii) { 745 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 746 if (format == PETSC_VIEWER_ASCII_INFO){ 747 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 748 } 749 } 750 PetscFunctionReturn(0); 751 } 752 753 #undef __FUNCT__ 754 #define __FUNCT__ "MatGetInfo_MUMPS" 755 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 756 { 757 Mat_MUMPS *lu =(Mat_MUMPS*)A->spptr; 758 759 PetscFunctionBegin; 760 info->block_size = 1.0; 761 info->nz_allocated = lu->id.INFOG(20); 762 info->nz_used = lu->id.INFOG(20); 763 info->nz_unneeded = 0.0; 764 info->assemblies = 0.0; 765 info->mallocs = 0.0; 766 info->memory = 0.0; 767 info->fill_ratio_given = 0; 768 info->fill_ratio_needed = 0; 769 info->factor_mallocs = 0; 770 PetscFunctionReturn(0); 771 } 772 773 /*MC 774 MAT_SOLVER_MUMPS - A matrix type providing direct solvers (LU and Cholesky) for 775 distributed and sequential matrices via the external package MUMPS. 776 777 Works with MATAIJ and MATSBAIJ matrices 778 779 Options Database Keys: 780 + -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 781 . -mat_mumps_icntl_4 <0,...,4> - print level 782 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 783 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 784 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 785 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 786 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 787 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 788 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 789 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 790 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 791 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 792 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 793 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 794 795 Level: beginner 796 797 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 798 799 M*/ 800 801 EXTERN_C_BEGIN 802 #undef __FUNCT__ 803 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 804 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 805 { 806 PetscFunctionBegin; 807 *type = MAT_SOLVER_MUMPS; 808 PetscFunctionReturn(0); 809 } 810 EXTERN_C_END 811 812 EXTERN_C_BEGIN 813 /* 814 The seq and mpi versions of this function are the same 815 */ 816 #undef __FUNCT__ 817 #define __FUNCT__ "MatGetFactor_seqaij_mumps" 818 PetscErrorCode MatGetFactor_seqaij_mumps(Mat A,MatFactorType ftype,Mat *F) 819 { 820 Mat B; 821 PetscErrorCode ierr; 822 Mat_MUMPS *mumps; 823 824 PetscFunctionBegin; 825 /* Create the factorization matrix */ 826 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 827 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 828 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 829 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 830 831 B->ops->view = MatView_MUMPS; 832 B->ops->getinfo = MatGetInfo_MUMPS; 833 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 834 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 835 if (ftype == MAT_FACTOR_LU) { 836 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 837 B->factortype = MAT_FACTOR_LU; 838 } else { 839 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 840 B->factortype = MAT_FACTOR_CHOLESKY; 841 } 842 843 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 844 mumps->CleanUpMUMPS = PETSC_FALSE; 845 mumps->isAIJ = PETSC_TRUE; 846 mumps->scat_rhs = PETSC_NULL; 847 mumps->scat_sol = PETSC_NULL; 848 mumps->nSolve = 0; 849 mumps->MatDestroy = B->ops->destroy; 850 B->ops->destroy = MatDestroy_MUMPS; 851 B->spptr = (void*)mumps; 852 853 *F = B; 854 PetscFunctionReturn(0); 855 } 856 EXTERN_C_END 857 858 EXTERN_C_BEGIN 859 #undef __FUNCT__ 860 #define __FUNCT__ "MatGetFactor_seqsbaij_mumps" 861 PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 862 { 863 Mat B; 864 PetscErrorCode ierr; 865 Mat_MUMPS *mumps; 866 867 PetscFunctionBegin; 868 if (ftype != MAT_FACTOR_CHOLESKY) { 869 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 870 } 871 /* Create the factorization matrix */ 872 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 873 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 874 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 875 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 876 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 877 878 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 879 B->ops->view = MatView_MUMPS; 880 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 881 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 882 B->factortype = MAT_FACTOR_CHOLESKY; 883 884 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 885 mumps->CleanUpMUMPS = PETSC_FALSE; 886 mumps->isAIJ = PETSC_FALSE; 887 mumps->scat_rhs = PETSC_NULL; 888 mumps->scat_sol = PETSC_NULL; 889 mumps->nSolve = 0; 890 mumps->MatDestroy = B->ops->destroy; 891 B->ops->destroy = MatDestroy_MUMPS; 892 B->spptr = (void*)mumps; 893 894 *F = B; 895 PetscFunctionReturn(0); 896 } 897 EXTERN_C_END 898 899 EXTERN_C_BEGIN 900 #undef __FUNCT__ 901 #define __FUNCT__ "MatGetFactor_mpisbaij_mumps" 902 PetscErrorCode MatGetFactor_mpisbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 903 { 904 Mat B; 905 PetscErrorCode ierr; 906 Mat_MUMPS *mumps; 907 908 PetscFunctionBegin; 909 if (ftype != MAT_FACTOR_CHOLESKY) { 910 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 911 } 912 /* Create the factorization matrix */ 913 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 914 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 915 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 916 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 917 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 918 919 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 920 B->ops->view = MatView_MUMPS; 921 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 922 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 923 B->factortype = MAT_FACTOR_CHOLESKY; 924 925 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 926 mumps->CleanUpMUMPS = PETSC_FALSE; 927 mumps->isAIJ = PETSC_FALSE; 928 mumps->scat_rhs = PETSC_NULL; 929 mumps->scat_sol = PETSC_NULL; 930 mumps->nSolve = 0; 931 mumps->MatDestroy = B->ops->destroy; 932 B->ops->destroy = MatDestroy_MUMPS; 933 B->spptr = (void*)mumps; 934 935 *F = B; 936 PetscFunctionReturn(0); 937 } 938 EXTERN_C_END 939 940 EXTERN_C_BEGIN 941 #undef __FUNCT__ 942 #define __FUNCT__ "MatGetFactor_mpiaij_mumps" 943 PetscErrorCode MatGetFactor_mpiaij_mumps(Mat A,MatFactorType ftype,Mat *F) 944 { 945 Mat B; 946 PetscErrorCode ierr; 947 Mat_MUMPS *mumps; 948 949 PetscFunctionBegin; 950 /* Create the factorization matrix */ 951 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 952 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 953 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 954 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 955 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 956 957 if (ftype == MAT_FACTOR_LU) { 958 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 959 B->factortype = MAT_FACTOR_LU; 960 } else { 961 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 962 B->factortype = MAT_FACTOR_CHOLESKY; 963 } 964 965 B->ops->view = MatView_MUMPS; 966 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 967 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 968 969 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 970 mumps->CleanUpMUMPS = PETSC_FALSE; 971 mumps->isAIJ = PETSC_TRUE; 972 mumps->scat_rhs = PETSC_NULL; 973 mumps->scat_sol = PETSC_NULL; 974 mumps->nSolve = 0; 975 mumps->MatDestroy = B->ops->destroy; 976 B->ops->destroy = MatDestroy_MUMPS; 977 B->spptr = (void*)mumps; 978 979 *F = B; 980 PetscFunctionReturn(0); 981 } 982 EXTERN_C_END 983 984 EXTERN_C_BEGIN 985 #undef __FUNCT__ 986 #define __FUNCT__ "MatGetFactor_mpibaij_mumps" 987 PetscErrorCode MatGetFactor_mpibaij_mumps(Mat A,MatFactorType ftype,Mat *F) 988 { 989 Mat B; 990 PetscErrorCode ierr; 991 Mat_MUMPS *mumps; 992 993 PetscFunctionBegin; 994 /* Create the factorization matrix */ 995 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 996 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 997 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 998 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 999 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1000 1001 if (ftype == MAT_FACTOR_LU) { 1002 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1003 B->factortype = MAT_FACTOR_LU; 1004 } 1005 else { 1006 SETERRQ(PETSC_ERR_SUP,"Cholesky factorization for Petsc BAIJ matrices not supported yet\n"); 1007 } 1008 B->ops->view = MatView_MUMPS; 1009 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1010 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1011 1012 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1013 mumps->CleanUpMUMPS = PETSC_FALSE; 1014 mumps->isAIJ = PETSC_TRUE; 1015 mumps->scat_rhs = PETSC_NULL; 1016 mumps->scat_sol = PETSC_NULL; 1017 mumps->nSolve = 0; 1018 mumps->MatDestroy = B->ops->destroy; 1019 B->ops->destroy = MatDestroy_MUMPS; 1020 B->spptr = (void*)mumps; 1021 1022 *F = B; 1023 PetscFunctionReturn(0); 1024 } 1025 EXTERN_C_END 1026 1027 /* -------------------------------------------------------------------------------------------*/ 1028 /*@ 1029 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1030 1031 Collective on Mat 1032 1033 Input Parameters: 1034 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1035 . idx - index of MUMPS parameter array ICNTL() 1036 - icntl - value of MUMPS ICNTL(imumps) 1037 1038 Options Database: 1039 . -mat_mumps_icntl_<idx> <icntl> 1040 1041 Level: beginner 1042 1043 References: MUMPS Users' Guide 1044 1045 .seealso: MatGetFactor() 1046 @*/ 1047 #undef __FUNCT__ 1048 #define __FUNCT__ "MatMumpsSetIcntl" 1049 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt idx,PetscInt icntl) 1050 { 1051 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 1052 1053 PetscFunctionBegin; 1054 lu->id.ICNTL(idx) = icntl; 1055 PetscFunctionReturn(0); 1056 } 1057 1058