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