1 /*$Id: mumps.c,v 1.10 2001/08/15 15:56:50 bsmith Exp $*/ 2 /* 3 Provides an interface to the MUMPS_4.2_beta sparse solver 4 */ 5 6 #include "src/mat/impls/aij/seq/aij.h" 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 RINFOG(I) rinfog[(I)-1] 25 26 typedef struct { 27 #if defined(PETSC_USE_COMPLEX) 28 ZMUMPS_STRUC_C id; 29 #else 30 DMUMPS_STRUC_C id; 31 #endif 32 MatStructure matstruc; 33 int myid,size,*irn,*jcn,sym; 34 PetscScalar *val; 35 MPI_Comm comm_mumps; 36 37 MatType basetype; 38 PetscTruth isAIJ,CleanUpMUMPS; 39 int (*MatView)(Mat,PetscViewer); 40 int (*MatAssemblyEnd)(Mat,MatAssemblyType); 41 int (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 42 int (*MatCholeskyFactorSymbolic)(Mat,IS,MatFactorInfo*,Mat*); 43 int (*MatDestroy)(Mat); 44 } Mat_AIJ_MUMPS; 45 46 /* convert Petsc mpiaij matrix to triples: row[nz], col[nz], val[nz] */ 47 /* 48 input: 49 A - matrix in mpiaij format 50 shift - 0: C style output triple; 1: Fortran style output triple. 51 valOnly - FALSE: spaces are allocated and values are set for the triple 52 TRUE: only the values in v array are updated 53 output: 54 nnz - dim of r, c, and v (number of local nonzero entries of A) 55 r, c, v - row and col index, matrix values (matrix triples) 56 */ 57 int MatConvertToTriples(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 58 { 59 int *ai, *aj, *bi, *bj, rstart,nz, *garray; 60 int ierr,i,j,jj,jB,irow,m=A->m,*ajj,*bjj,countA,countB,colA_start,jcol; 61 int *row,*col; 62 PetscScalar *av, *bv,*val; 63 Mat_AIJ_MUMPS *mumps = (Mat_AIJ_MUMPS *)A->spptr; 64 65 PetscFunctionBegin; 66 67 if (mumps->isAIJ){ 68 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 69 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 70 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 71 nz = aa->nz + bb->nz; 72 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= mat->rstart; 73 garray = mat->garray; 74 av=aa->a; bv=bb->a; 75 76 } else { 77 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 78 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)(mat->A)->data; 79 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 80 if (mat->bs > 1) SETERRQ1(PETSC_ERR_SUP," bs=%d is not supported yet\n", mat->bs); 81 nz = aa->s_nz + bb->nz; 82 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= mat->rstart; 83 garray = mat->garray; 84 av=aa->a; bv=bb->a; 85 } 86 87 if (!valOnly){ 88 ierr = PetscMalloc(nz*sizeof(int),&row);CHKERRQ(ierr); 89 ierr = PetscMalloc(nz*sizeof(int),&col);CHKERRQ(ierr); 90 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 91 *r = row; *c = col; *v = val; 92 } else { 93 row = *r; col = *c; val = *v; 94 } 95 *nnz = nz; 96 97 jj = 0; jB = 0; irow = rstart; 98 for ( i=0; i<m; i++ ) { 99 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 100 countA = ai[i+1] - ai[i]; 101 countB = bi[i+1] - bi[i]; 102 bjj = bj + bi[i]; 103 104 /* get jB, the starting local col index for the 2nd B-part */ 105 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 106 for (j=0; j<countB; j++){ 107 jcol = garray[bjj[j]]; 108 if (jcol > colA_start) { jB = j; break; } 109 if (j==countB-1) jB = countB; 110 } 111 112 /* B-part, smaller col index */ 113 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 114 for (j=0; j<jB; j++){ 115 jcol = garray[bjj[j]]; 116 if (!valOnly){ 117 row[jj] = irow + shift; col[jj] = jcol + shift; 118 } 119 val[jj++] = *bv++; 120 } 121 /* A-part */ 122 for (j=0; j<countA; j++){ 123 if (!valOnly){ 124 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 125 } 126 val[jj++] = *av++; 127 } 128 /* B-part, larger col index */ 129 for (j=jB; j<countB; j++){ 130 if (!valOnly){ 131 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 132 } 133 val[jj++] = *bv++; 134 } 135 irow++; 136 } 137 138 PetscFunctionReturn(0); 139 } 140 141 EXTERN_C_BEGIN 142 #undef __FUNCT__ 143 #define __FUNCT__ "MatConvert_MUMPS_Base" 144 int MatConvert_MUMPS_Base(Mat A,MatType type,Mat *newmat) { 145 /* This routine is only called to convert an unfactored PETSc-MUMPS matrix */ 146 /* to its base PETSc type, so we will ignore 'MatType type'. */ 147 int ierr; 148 Mat B=*newmat; 149 Mat_AIJ_MUMPS *lu=(Mat_AIJ_MUMPS*)A->spptr; 150 151 PetscFunctionBegin; 152 if (B != A) { 153 /* This routine was inherited from SeqAIJ. */ 154 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 155 } else { 156 157 B->ops->view = lu->MatView; 158 B->ops->assemblyend = lu->MatAssemblyEnd; 159 B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic; 160 B->ops->choleskyfactorsymbolic = lu->MatCholeskyFactorSymbolic; 161 B->ops->destroy = lu->MatDestroy; 162 ierr = PetscObjectChangeTypeName((PetscObject)B,lu->basetype);CHKERRQ(ierr); 163 ierr = PetscFree(lu);CHKERRQ(ierr); 164 } 165 *newmat = B; 166 PetscFunctionReturn(0); 167 } 168 EXTERN_C_END 169 170 #undef __FUNCT__ 171 #define __FUNCT__ "MatDestroy_AIJ_MUMPS" 172 int MatDestroy_AIJ_MUMPS(Mat A) 173 { 174 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)A->spptr; 175 int ierr,size=lu->size; 176 177 PetscFunctionBegin; 178 if (lu->CleanUpMUMPS) { 179 /* Terminate instance, deallocate memories */ 180 lu->id.job=JOB_END; 181 #if defined(PETSC_USE_COMPLEX) 182 zmumps_c(&lu->id); 183 #else 184 dmumps_c(&lu->id); 185 #endif 186 if (lu->irn) { 187 ierr = PetscFree(lu->irn);CHKERRQ(ierr); 188 } 189 if (lu->jcn) { 190 ierr = PetscFree(lu->jcn);CHKERRQ(ierr); 191 } 192 if (size>1 && lu->val) { 193 ierr = PetscFree(lu->val);CHKERRQ(ierr); 194 } 195 ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr); 196 } 197 ierr = MatConvert_MUMPS_Base(A,lu->basetype,&A);CHKERRQ(ierr); 198 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "MatFactorInfo_MUMPS" 204 int MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 205 { 206 Mat_AIJ_MUMPS *lu= (Mat_AIJ_MUMPS*)A->spptr; 207 int ierr; 208 209 PetscFunctionBegin; 210 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 211 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 212 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 214 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 215 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 216 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 217 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 218 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 219 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 220 if (lu->myid == 0 && lu->id.ICNTL(11)>0) { 221 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 222 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 223 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 224 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); 225 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 226 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 227 228 } 229 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 230 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 231 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (efficiency control): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 232 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(15) (efficiency control): %d \n",lu->id.ICNTL(15));CHKERRQ(ierr); 233 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 234 235 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 236 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 237 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 238 PetscFunctionReturn(0); 239 } 240 241 #undef __FUNCT__ 242 #define __FUNCT__ "MatView_AIJ_MUMPS" 243 int MatView_AIJ_MUMPS(Mat A,PetscViewer viewer) { 244 int ierr; 245 PetscTruth isascii; 246 PetscViewerFormat format; 247 Mat_AIJ_MUMPS *mumps=(Mat_AIJ_MUMPS*)(A->spptr); 248 249 PetscFunctionBegin; 250 ierr = (*mumps->MatView)(A,viewer);CHKERRQ(ierr); 251 252 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 253 if (isascii) { 254 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 255 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 256 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 257 } 258 } 259 PetscFunctionReturn(0); 260 } 261 262 #undef __FUNCT__ 263 #define __FUNCT__ "MatSolve_AIJ_MUMPS" 264 int MatSolve_AIJ_MUMPS(Mat A,Vec b,Vec x) 265 { 266 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)A->spptr; 267 PetscScalar *array; 268 Vec x_seq; 269 IS iden; 270 VecScatter scat; 271 int ierr; 272 273 PetscFunctionBegin; 274 if (lu->size > 1){ 275 if (!lu->myid){ 276 ierr = VecCreateSeq(PETSC_COMM_SELF,A->N,&x_seq);CHKERRQ(ierr); 277 ierr = ISCreateStride(PETSC_COMM_SELF,A->N,0,1,&iden);CHKERRQ(ierr); 278 } else { 279 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&x_seq);CHKERRQ(ierr); 280 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&iden);CHKERRQ(ierr); 281 } 282 ierr = VecScatterCreate(b,iden,x_seq,iden,&scat);CHKERRQ(ierr); 283 ierr = ISDestroy(iden);CHKERRQ(ierr); 284 285 ierr = VecScatterBegin(b,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 286 ierr = VecScatterEnd(b,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 287 if (!lu->myid) {ierr = VecGetArray(x_seq,&array);CHKERRQ(ierr);} 288 } else { /* size == 1 */ 289 ierr = VecCopy(b,x);CHKERRQ(ierr); 290 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 291 } 292 if (!lu->myid) { /* define rhs on the host */ 293 #if defined(PETSC_USE_COMPLEX) 294 lu->id.rhs = (mumps_double_complex*)array; 295 #else 296 lu->id.rhs = array; 297 #endif 298 } 299 300 /* solve phase */ 301 lu->id.job=3; 302 #if defined(PETSC_USE_COMPLEX) 303 zmumps_c(&lu->id); 304 #else 305 dmumps_c(&lu->id); 306 #endif 307 if (lu->id.INFOG(1) < 0) { 308 SETERRQ1(1,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 309 } 310 311 /* convert mumps solution x_seq to petsc mpi x */ 312 if (lu->size > 1) { 313 if (!lu->myid){ 314 ierr = VecRestoreArray(x_seq,&array);CHKERRQ(ierr); 315 } 316 ierr = VecScatterBegin(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 317 ierr = VecScatterEnd(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 318 ierr = VecScatterDestroy(scat);CHKERRQ(ierr); 319 ierr = VecDestroy(x_seq);CHKERRQ(ierr); 320 } else { 321 ierr = VecRestoreArray(x,&array);CHKERRQ(ierr); 322 } 323 324 PetscFunctionReturn(0); 325 } 326 327 #undef __FUNCT__ 328 #define __FUNCT__ "MatFactorNumeric_MPIAIJ_MUMPS" 329 int MatFactorNumeric_AIJ_MUMPS(Mat A,Mat *F) 330 { 331 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)(*F)->spptr; 332 Mat_AIJ_MUMPS *lua = (Mat_AIJ_MUMPS*)(A)->spptr; 333 int rnz,nnz,ierr,nz,i,M=A->M,*ai,*aj,icntl; 334 PetscTruth valOnly,flg; 335 336 PetscFunctionBegin; 337 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 338 (*F)->ops->solve = MatSolve_AIJ_MUMPS; 339 340 /* Initialize a MUMPS instance */ 341 ierr = MPI_Comm_rank(A->comm, &lu->myid); 342 ierr = MPI_Comm_size(A->comm,&lu->size);CHKERRQ(ierr); 343 lu->id.job = JOB_INIT; 344 ierr = MPI_Comm_dup(A->comm,&(lu->comm_mumps));CHKERRQ(ierr); 345 lu->id.comm_fortran = lu->comm_mumps; 346 347 /* Set mumps options */ 348 ierr = PetscOptionsBegin(A->comm,A->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 349 lu->id.par=1; /* host participates factorizaton and solve */ 350 lu->id.sym=lu->sym; 351 if (lu->sym == 2){ 352 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 353 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 354 } 355 #if defined(PETSC_USE_COMPLEX) 356 zmumps_c(&lu->id); 357 #else 358 dmumps_c(&lu->id); 359 #endif 360 361 if (lu->size == 1){ 362 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 363 } else { 364 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 365 } 366 367 icntl=-1; 368 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 369 if (flg && icntl > 0) { 370 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 371 } else { /* no output */ 372 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 373 lu->id.ICNTL(2) = -1; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 374 lu->id.ICNTL(3) = -1; /* output stream for global information, default=6 */ 375 lu->id.ICNTL(4) = 0; /* level of printing, 0,1,2,3,4, default=2 */ 376 } 377 ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): matrix prescaling (0 to 7)","None",lu->id.ICNTL(6),&lu->id.ICNTL(6),PETSC_NULL);CHKERRQ(ierr); 378 icntl=-1; 379 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 380 if (flg) { 381 if (icntl== 1){ 382 SETERRQ(PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n"); 383 } else { 384 lu->id.ICNTL(7) = icntl; 385 } 386 } 387 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); 388 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); 389 ierr = PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): error analysis, a positive value returns statistics (by -sles_view)","None",lu->id.ICNTL(11),&lu->id.ICNTL(11),PETSC_NULL);CHKERRQ(ierr); 390 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr); 391 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control","None",lu->id.ICNTL(13),&lu->id.ICNTL(13),PETSC_NULL);CHKERRQ(ierr); 392 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): efficiency control","None",lu->id.ICNTL(14),&lu->id.ICNTL(14),PETSC_NULL);CHKERRQ(ierr); 393 ierr = PetscOptionsInt("-mat_mumps_icntl_15","ICNTL(15): efficiency control","None",lu->id.ICNTL(15),&lu->id.ICNTL(15),PETSC_NULL);CHKERRQ(ierr); 394 395 /* 396 ierr = PetscOptionsInt("-mat_mumps_icntl_16","ICNTL(16): 1: rank detection; 2: rank detection and nullspace","None",lu->id.ICNTL(16),&icntl,&flg);CHKERRQ(ierr); 397 if (flg){ 398 if (icntl >-1 && icntl <3 ){ 399 if (lu->myid==0) lu->id.ICNTL(16) = icntl; 400 } else { 401 SETERRQ1(PETSC_ERR_SUP,"ICNTL(16)=%d -- not supported\n",icntl); 402 } 403 } 404 */ 405 406 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 407 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); 408 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 409 PetscOptionsEnd(); 410 } 411 412 /* define matrix A */ 413 switch (lu->id.ICNTL(18)){ 414 case 0: /* centralized assembled matrix input (size=1) */ 415 if (!lu->myid) { 416 if (lua->isAIJ){ 417 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 418 nz = aa->nz; 419 ai = aa->i; aj = aa->j; lu->val = aa->a; 420 } else { 421 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 422 nz = aa->s_nz; 423 ai = aa->i; aj = aa->j; lu->val = aa->a; 424 } 425 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization, get irn and jcn */ 426 ierr = PetscMalloc(nz*sizeof(int),&lu->irn);CHKERRQ(ierr); 427 ierr = PetscMalloc(nz*sizeof(int),&lu->jcn);CHKERRQ(ierr); 428 nz = 0; 429 for (i=0; i<M; i++){ 430 rnz = ai[i+1] - ai[i]; 431 while (rnz--) { /* Fortran row/col index! */ 432 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 433 } 434 } 435 } 436 } 437 break; 438 case 3: /* distributed assembled matrix input (size>1) */ 439 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 440 valOnly = PETSC_FALSE; 441 } else { 442 valOnly = PETSC_TRUE; /* only update mat values, not row and col index */ 443 } 444 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 445 break; 446 default: SETERRQ(PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 447 } 448 449 /* analysis phase */ 450 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 451 lu->id.n = M; 452 switch (lu->id.ICNTL(18)){ 453 case 0: /* centralized assembled matrix input */ 454 if (!lu->myid) { 455 lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 456 if (lu->id.ICNTL(6)>1){ 457 #if defined(PETSC_USE_COMPLEX) 458 lu->id.a = (mumps_double_complex*)lu->val; 459 #else 460 lu->id.a = lu->val; 461 #endif 462 } 463 } 464 break; 465 case 3: /* distributed assembled matrix input (size>1) */ 466 lu->id.nz_loc = nnz; 467 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 468 if (lu->id.ICNTL(6)>1) { 469 #if defined(PETSC_USE_COMPLEX) 470 lu->id.a_loc = (mumps_double_complex*)lu->val; 471 #else 472 lu->id.a_loc = lu->val; 473 #endif 474 } 475 break; 476 } 477 lu->id.job=1; 478 #if defined(PETSC_USE_COMPLEX) 479 zmumps_c(&lu->id); 480 #else 481 dmumps_c(&lu->id); 482 #endif 483 if (lu->id.INFOG(1) < 0) { 484 SETERRQ1(1,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 485 } 486 } 487 488 /* numerical factorization phase */ 489 if(lu->id.ICNTL(18) == 0) { 490 if (lu->myid == 0) { 491 #if defined(PETSC_USE_COMPLEX) 492 lu->id.a = (mumps_double_complex*)lu->val; 493 #else 494 lu->id.a = lu->val; 495 #endif 496 } 497 } else { 498 #if defined(PETSC_USE_COMPLEX) 499 lu->id.a_loc = (mumps_double_complex*)lu->val; 500 #else 501 lu->id.a_loc = lu->val; 502 #endif 503 } 504 lu->id.job=2; 505 #if defined(PETSC_USE_COMPLEX) 506 zmumps_c(&lu->id); 507 #else 508 dmumps_c(&lu->id); 509 #endif 510 if (lu->id.INFOG(1) < 0) { 511 SETERRQ1(1,"1, Error reported by MUMPS in numerical factorization phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 512 } 513 514 if (lu->myid==0 && lu->id.ICNTL(16) > 0){ 515 SETERRQ1(1," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 516 } 517 518 (*F)->assembled = PETSC_TRUE; 519 lu->matstruc = SAME_NONZERO_PATTERN; 520 lu->CleanUpMUMPS = PETSC_TRUE; 521 PetscFunctionReturn(0); 522 } 523 524 /* Note the Petsc r and c permutations are ignored */ 525 #undef __FUNCT__ 526 #define __FUNCT__ "MatLUFactorSymbolic_AIJ_MUMPS" 527 int MatLUFactorSymbolic_AIJ_MUMPS(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 528 { 529 Mat B; 530 Mat_AIJ_MUMPS *lu; 531 int ierr; 532 533 PetscFunctionBegin; 534 535 /* Create the factorization matrix */ 536 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 537 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 538 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 539 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 540 541 B->ops->lufactornumeric = MatFactorNumeric_AIJ_MUMPS; 542 B->factor = FACTOR_LU; 543 lu = (Mat_AIJ_MUMPS*)B->spptr; 544 lu->sym = 0; 545 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 546 547 *F = B; 548 PetscFunctionReturn(0); 549 } 550 551 /* Note the Petsc r permutation is ignored */ 552 #undef __FUNCT__ 553 #define __FUNCT__ "MatCholeskyFactorSymbolic_AIJ_MUMPS" 554 int MatCholeskyFactorSymbolic_AIJ_MUMPS(Mat A,IS r,MatFactorInfo *info,Mat *F) 555 { 556 Mat B; 557 Mat_AIJ_MUMPS *lu; 558 int ierr; 559 560 PetscFunctionBegin; 561 562 /* Create the factorization matrix */ 563 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 564 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 565 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 566 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 567 568 B->ops->choleskyfactornumeric = MatFactorNumeric_AIJ_MUMPS; 569 B->factor = FACTOR_CHOLESKY; 570 lu = (Mat_AIJ_MUMPS *)B->spptr; 571 lu->sym = 2; 572 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 573 574 *F = B; 575 PetscFunctionReturn(0); 576 } 577 578 #undef __FUNCT__ 579 #define __FUNCT__ "MatAssemblyEnd_AIJ_MUMPS" 580 int MatAssemblyEnd_AIJ_MUMPS(Mat A,MatAssemblyType mode) { 581 int ierr; 582 Mat_AIJ_MUMPS *mumps=(Mat_AIJ_MUMPS*)A->spptr; 583 584 PetscFunctionBegin; 585 ierr = (*mumps->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 586 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 587 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 588 A->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJ_MUMPS; 589 PetscFunctionReturn(0); 590 } 591 592 EXTERN_C_BEGIN 593 #undef __FUNCT__ 594 #define __FUNCT__ "MatConvert_AIJ_MUMPS" 595 int MatConvert_AIJ_MUMPS(Mat A,MatType newtype,Mat *newmat) { 596 int ierr,size; 597 MPI_Comm comm; 598 Mat B=*newmat; 599 Mat_AIJ_MUMPS *mumps; 600 601 PetscFunctionBegin; 602 if (B != A) { 603 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 604 } 605 606 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 607 ierr = PetscNew(Mat_AIJ_MUMPS,&mumps);CHKERRQ(ierr); 608 609 mumps->MatView = A->ops->view; 610 mumps->MatAssemblyEnd = A->ops->assemblyend; 611 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 612 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 613 mumps->MatDestroy = A->ops->destroy; 614 mumps->CleanUpMUMPS = PETSC_FALSE; 615 mumps->isAIJ = PETSC_TRUE; 616 617 B->spptr = (void *)mumps; 618 B->ops->view = MatView_AIJ_MUMPS; 619 B->ops->assemblyend = MatAssemblyEnd_AIJ_MUMPS; 620 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJ_MUMPS; 621 B->ops->destroy = MatDestroy_AIJ_MUMPS; 622 623 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 624 if (size == 1) { 625 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_aijmumps_C", 626 "MatConvert_AIJ_MUMPS",MatConvert_AIJ_MUMPS);CHKERRQ(ierr); 627 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_aijmumps_seqaij_C", 628 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 629 } else { 630 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpiaij_aijmumps_C", 631 "MatConvert_AIJ_MUMPS",MatConvert_AIJ_MUMPS);CHKERRQ(ierr); 632 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_aijmumps_mpiaij_C", 633 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 634 } 635 636 PetscLogInfo(0,"Using MUMPS for LU factorization and solves."); 637 ierr = PetscObjectChangeTypeName((PetscObject)B,newtype);CHKERRQ(ierr); 638 *newmat = B; 639 PetscFunctionReturn(0); 640 } 641 EXTERN_C_END 642 643 /*MC 644 MATAIJMUMPS - a matrix type providing direct solvers (LU) for distributed 645 and sequential matrices via the external package MUMPS. 646 647 If MUMPS is installed (see the manual for instructions 648 on how to declare the existence of external packages), 649 a matrix type can be constructed which invokes MUMPS solvers. 650 After calling MatCreate(...,A), simply call MatSetType(A,MATAIJMUMPS). 651 This matrix type is only supported for double precision real. 652 653 If created with a single process communicator, this matrix type inherits from MATSEQAIJ. 654 Otherwise, this matrix type inherits from MATMPIAIJ. Hence for single process communicators, 655 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 656 for communicators controlling multiple processes. It is recommended that you call both of 657 the above preallocation routines for simplicity. 658 659 Options Database Keys: 660 + -mat_type aijmumps 661 . -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 662 . -mat_mumps_icntl_4 <0,1,2,3,4> - print level 663 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 664 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 665 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 666 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 667 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -sles_view 668 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 669 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 670 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 671 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 672 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 673 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 674 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 675 676 Level: beginner 677 678 .seealso: MATSBAIJMUMPS 679 M*/ 680 681 EXTERN_C_BEGIN 682 #undef __FUNCT__ 683 #define __FUNCT__ "MatCreate_AIJ_MUMPS" 684 int MatCreate_AIJ_MUMPS(Mat A) { 685 int ierr,size; 686 MPI_Comm comm; 687 688 PetscFunctionBegin; 689 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 690 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 691 if (size == 1) { 692 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 693 } else { 694 ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr); 695 } 696 ierr = MatConvert_AIJ_MUMPS(A,MATAIJMUMPS,&A);CHKERRQ(ierr); 697 PetscFunctionReturn(0); 698 } 699 EXTERN_C_END 700 701 EXTERN_C_BEGIN 702 #undef __FUNCT__ 703 #define __FUNCT__ "MatLoad_AIJ_MUMPS" 704 int MatLoad_AIJ_MUMPS(PetscViewer viewer,MatType type,Mat *A) { 705 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 706 MPI_Comm comm; 707 708 PetscFunctionBegin; 709 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 710 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 711 if (size == 1) { 712 ierr = PetscFListFind(comm,MatLoadList,MATSEQAIJ,(void(**)(void))&r);CHKERRQ(ierr); 713 } else { 714 ierr = PetscFListFind(comm,MatLoadList,MATMPIAIJ,(void(**)(void))&r);CHKERRQ(ierr); 715 } 716 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 717 PetscFunctionReturn(0); 718 } 719 EXTERN_C_END 720 721 #undef __FUNCT__ 722 #define __FUNCT__ "MatAssemblyEnd_AIJ_MUMPS" 723 int MatAssemblyEnd_SBAIJ_MUMPS(Mat A,MatAssemblyType mode) { 724 int ierr; 725 Mat_AIJ_MUMPS *mumps=(Mat_AIJ_MUMPS*)A->spptr; 726 727 PetscFunctionBegin; 728 ierr = (*mumps->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 729 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 730 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 731 A->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJ_MUMPS; 732 PetscFunctionReturn(0); 733 } 734 735 EXTERN_C_BEGIN 736 #undef __FUNCT__ 737 #define __FUNCT__ "MatConvert_SBAIJ_MUMPS" 738 int MatConvert_SBAIJ_MUMPS(Mat A,MatType newtype,Mat *newmat) { 739 int ierr,size; 740 MPI_Comm comm; 741 Mat B=*newmat; 742 Mat_AIJ_MUMPS *mumps; 743 744 PetscFunctionBegin; 745 if (B != A) { 746 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 747 } 748 749 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 750 ierr = PetscNew(Mat_AIJ_MUMPS,&mumps);CHKERRQ(ierr); 751 752 mumps->MatView = A->ops->view; 753 mumps->MatAssemblyEnd = A->ops->assemblyend; 754 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 755 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 756 mumps->MatDestroy = A->ops->destroy; 757 mumps->CleanUpMUMPS = PETSC_FALSE; 758 mumps->isAIJ = PETSC_FALSE; 759 760 B->spptr = (void *)mumps; 761 B->ops->view = MatView_AIJ_MUMPS; 762 B->ops->assemblyend = MatAssemblyEnd_SBAIJ_MUMPS; 763 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJ_MUMPS; 764 B->ops->destroy = MatDestroy_AIJ_MUMPS; 765 766 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 767 if (size == 1) { 768 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_mumps_C", 769 "MatConvert_SBAIJ_MUMPS",MatConvert_SBAIJ_MUMPS);CHKERRQ(ierr); 770 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mumps_seqsbaij_C", 771 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 772 } else { 773 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpisbaij_mumps_C", 774 "MatConvert_SBAIJ_MUMPS",MatConvert_SBAIJ_MUMPS);CHKERRQ(ierr); 775 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mumps_mpisbaij_C", 776 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 777 } 778 779 PetscLogInfo(0,"Using MUMPS for Cholesky factorization and solves."); 780 ierr = PetscObjectChangeTypeName((PetscObject)B,newtype);CHKERRQ(ierr); 781 *newmat = B; 782 PetscFunctionReturn(0); 783 } 784 EXTERN_C_END 785 786 /*MC 787 MATSBAIJMUMPS - a symmetric matrix type providing direct solvers (Cholesky) for 788 distributed and sequential matrices via the external package MUMPS. 789 790 If MUMPS is installed (see the manual for instructions 791 on how to declare the existence of external packages), 792 a matrix type can be constructed which invokes MUMPS solvers. 793 After calling MatCreate(...,A), simply call MatSetType(A,MATSBAIJMUMPS). 794 This matrix type is only supported for double precision real. 795 796 If created with a single process communicator, this matrix type inherits from MATSEQSBAIJ. 797 Otherwise, this matrix type inherits from MATMPISBAIJ. Hence for single process communicators, 798 MatSeqSBAIJSetPreallocation is supported, and similarly MatMPISBAIJSetPreallocation is supported 799 for communicators controlling multiple processes. It is recommended that you call both of 800 the above preallocation routines for simplicity. 801 802 Options Database Keys: 803 + -mat_type aijmumps 804 . -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 805 . -mat_mumps_icntl_4 <0,...,4> - print level 806 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 807 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 808 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 809 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 810 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -sles_view 811 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 812 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 813 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 814 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 815 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 816 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 817 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 818 819 Level: beginner 820 821 .seealso: MATAIJMUMPS 822 M*/ 823 824 EXTERN_C_BEGIN 825 #undef __FUNCT__ 826 #define __FUNCT__ "MatCreate_SBAIJ_MUMPS" 827 int MatCreate_SBAIJ_MUMPS(Mat A) { 828 int ierr,size; 829 MPI_Comm comm; 830 831 PetscFunctionBegin; 832 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 833 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 834 if (size == 1) { 835 ierr = MatSetType(A,MATSEQSBAIJ);CHKERRQ(ierr); 836 } else { 837 ierr = MatSetType(A,MATMPISBAIJ);CHKERRQ(ierr); 838 } 839 ierr = MatConvert_SBAIJ_MUMPS(A,MATSBAIJMUMPS,&A);CHKERRQ(ierr); 840 PetscFunctionReturn(0); 841 } 842 EXTERN_C_END 843 844 EXTERN_C_BEGIN 845 #undef __FUNCT__ 846 #define __FUNCT__ "MatLoad_SBAIJ_MUMPS" 847 int MatLoad_SBAIJ_MUMPS(PetscViewer viewer,MatType type,Mat *A) { 848 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 849 MPI_Comm comm; 850 851 PetscFunctionBegin; 852 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 853 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 854 if (size == 1) { 855 ierr = PetscFListFind(comm,MatLoadList,MATSEQSBAIJ,(void(**)(void))&r);CHKERRQ(ierr); 856 } else { 857 ierr = PetscFListFind(comm,MatLoadList,MATMPISBAIJ,(void(**)(void))&r);CHKERRQ(ierr); 858 } 859 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 860 PetscFunctionReturn(0); 861 } 862 EXTERN_C_END 863