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