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 MatType basetype; 37 PetscTruth isAIJ,CleanUpMUMPS; 38 int (*MatDuplicate)(Mat,MatDuplicateOption,Mat*); 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_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 = PetscStrfree(mumps->basetype);CHKERRQ(ierr); 165 ierr = PetscFree(mumps);CHKERRQ(ierr); 166 *newmat = B; 167 PetscFunctionReturn(0); 168 } 169 EXTERN_C_END 170 171 #undef __FUNCT__ 172 #define __FUNCT__ "MatDestroy_MUMPS" 173 int MatDestroy_MUMPS(Mat A) { 174 Mat_MUMPS *lu=(Mat_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 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 206 int ierr; 207 208 PetscFunctionBegin; 209 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 210 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 211 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 212 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 214 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 215 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 216 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 217 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 218 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 219 if (lu->myid == 0 && lu->id.ICNTL(11)>0) { 220 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 221 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 222 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 223 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); 224 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 225 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 226 227 } 228 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 229 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 230 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (efficiency control): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 231 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(15) (efficiency control): %d \n",lu->id.ICNTL(15));CHKERRQ(ierr); 232 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 233 234 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 235 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 236 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "MatView_AIJMUMPS" 242 int MatView_AIJMUMPS(Mat A,PetscViewer viewer) { 243 int ierr; 244 PetscTruth isascii; 245 PetscViewerFormat format; 246 Mat_MUMPS *mumps=(Mat_MUMPS*)(A->spptr); 247 248 PetscFunctionBegin; 249 ierr = (*mumps->MatView)(A,viewer);CHKERRQ(ierr); 250 251 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 252 if (isascii) { 253 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 254 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 255 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 256 } 257 } 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "MatSolve_AIJMUMPS" 263 int MatSolve_AIJMUMPS(Mat A,Vec b,Vec x) { 264 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 265 PetscScalar *array; 266 Vec x_seq; 267 IS iden; 268 VecScatter scat; 269 int ierr; 270 271 PetscFunctionBegin; 272 if (lu->size > 1){ 273 if (!lu->myid){ 274 ierr = VecCreateSeq(PETSC_COMM_SELF,A->N,&x_seq);CHKERRQ(ierr); 275 ierr = ISCreateStride(PETSC_COMM_SELF,A->N,0,1,&iden);CHKERRQ(ierr); 276 } else { 277 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&x_seq);CHKERRQ(ierr); 278 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&iden);CHKERRQ(ierr); 279 } 280 ierr = VecScatterCreate(b,iden,x_seq,iden,&scat);CHKERRQ(ierr); 281 ierr = ISDestroy(iden);CHKERRQ(ierr); 282 283 ierr = VecScatterBegin(b,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 284 ierr = VecScatterEnd(b,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 285 if (!lu->myid) {ierr = VecGetArray(x_seq,&array);CHKERRQ(ierr);} 286 } else { /* size == 1 */ 287 ierr = VecCopy(b,x);CHKERRQ(ierr); 288 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 289 } 290 if (!lu->myid) { /* define rhs on the host */ 291 #if defined(PETSC_USE_COMPLEX) 292 lu->id.rhs = (mumps_double_complex*)array; 293 #else 294 lu->id.rhs = array; 295 #endif 296 } 297 298 /* solve phase */ 299 lu->id.job=3; 300 #if defined(PETSC_USE_COMPLEX) 301 zmumps_c(&lu->id); 302 #else 303 dmumps_c(&lu->id); 304 #endif 305 if (lu->id.INFOG(1) < 0) { 306 SETERRQ1(1,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 307 } 308 309 /* convert mumps solution x_seq to petsc mpi x */ 310 if (lu->size > 1) { 311 if (!lu->myid){ 312 ierr = VecRestoreArray(x_seq,&array);CHKERRQ(ierr); 313 } 314 ierr = VecScatterBegin(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 315 ierr = VecScatterEnd(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 316 ierr = VecScatterDestroy(scat);CHKERRQ(ierr); 317 ierr = VecDestroy(x_seq);CHKERRQ(ierr); 318 } else { 319 ierr = VecRestoreArray(x,&array);CHKERRQ(ierr); 320 } 321 322 PetscFunctionReturn(0); 323 } 324 325 #undef __FUNCT__ 326 #define __FUNCT__ "MatFactorNumeric_MPIAIJMUMPS" 327 int MatFactorNumeric_AIJMUMPS(Mat A,Mat *F) { 328 Mat_MUMPS *lu =(Mat_MUMPS*)(*F)->spptr; 329 Mat_MUMPS *lua=(Mat_MUMPS*)(A)->spptr; 330 int rnz,nnz,ierr,nz,i,M=A->M,*ai,*aj,icntl; 331 PetscTruth valOnly,flg; 332 333 PetscFunctionBegin; 334 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 335 (*F)->ops->solve = MatSolve_AIJMUMPS; 336 337 /* Initialize a MUMPS instance */ 338 ierr = MPI_Comm_rank(A->comm, &lu->myid); 339 ierr = MPI_Comm_size(A->comm,&lu->size);CHKERRQ(ierr); 340 lua->myid = lu->myid; lua->size = lu->size; 341 lu->id.job = JOB_INIT; 342 ierr = MPI_Comm_dup(A->comm,&(lu->comm_mumps));CHKERRQ(ierr); 343 lu->id.comm_fortran = lu->comm_mumps; 344 345 /* Set mumps options */ 346 ierr = PetscOptionsBegin(A->comm,A->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 347 lu->id.par=1; /* host participates factorizaton and solve */ 348 lu->id.sym=lu->sym; 349 if (lu->sym == 2){ 350 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 351 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 352 } 353 #if defined(PETSC_USE_COMPLEX) 354 zmumps_c(&lu->id); 355 #else 356 dmumps_c(&lu->id); 357 #endif 358 359 if (lu->size == 1){ 360 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 361 } else { 362 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 363 } 364 365 icntl=-1; 366 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 367 if (flg && icntl > 0) { 368 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 369 } else { /* no output */ 370 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 371 lu->id.ICNTL(2) = -1; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 372 lu->id.ICNTL(3) = -1; /* output stream for global information, default=6 */ 373 lu->id.ICNTL(4) = 0; /* level of printing, 0,1,2,3,4, default=2 */ 374 } 375 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); 376 icntl=-1; 377 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 378 if (flg) { 379 if (icntl== 1){ 380 SETERRQ(PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n"); 381 } else { 382 lu->id.ICNTL(7) = icntl; 383 } 384 } 385 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); 386 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); 387 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); 388 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr); 389 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control","None",lu->id.ICNTL(13),&lu->id.ICNTL(13),PETSC_NULL);CHKERRQ(ierr); 390 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): efficiency control","None",lu->id.ICNTL(14),&lu->id.ICNTL(14),PETSC_NULL);CHKERRQ(ierr); 391 ierr = PetscOptionsInt("-mat_mumps_icntl_15","ICNTL(15): efficiency control","None",lu->id.ICNTL(15),&lu->id.ICNTL(15),PETSC_NULL);CHKERRQ(ierr); 392 393 /* 394 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); 395 if (flg){ 396 if (icntl >-1 && icntl <3 ){ 397 if (lu->myid==0) lu->id.ICNTL(16) = icntl; 398 } else { 399 SETERRQ1(PETSC_ERR_SUP,"ICNTL(16)=%d -- not supported\n",icntl); 400 } 401 } 402 */ 403 404 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 405 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); 406 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 407 PetscOptionsEnd(); 408 } 409 410 /* define matrix A */ 411 switch (lu->id.ICNTL(18)){ 412 case 0: /* centralized assembled matrix input (size=1) */ 413 if (!lu->myid) { 414 if (lua->isAIJ){ 415 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 416 nz = aa->nz; 417 ai = aa->i; aj = aa->j; lu->val = aa->a; 418 } else { 419 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 420 nz = aa->s_nz; 421 ai = aa->i; aj = aa->j; lu->val = aa->a; 422 } 423 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization, get irn and jcn */ 424 ierr = PetscMalloc(nz*sizeof(int),&lu->irn);CHKERRQ(ierr); 425 ierr = PetscMalloc(nz*sizeof(int),&lu->jcn);CHKERRQ(ierr); 426 nz = 0; 427 for (i=0; i<M; i++){ 428 rnz = ai[i+1] - ai[i]; 429 while (rnz--) { /* Fortran row/col index! */ 430 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 431 } 432 } 433 } 434 } 435 break; 436 case 3: /* distributed assembled matrix input (size>1) */ 437 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 438 valOnly = PETSC_FALSE; 439 } else { 440 valOnly = PETSC_TRUE; /* only update mat values, not row and col index */ 441 } 442 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 443 break; 444 default: SETERRQ(PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 445 } 446 447 /* analysis phase */ 448 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 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 break; 474 } 475 lu->id.job=1; 476 #if defined(PETSC_USE_COMPLEX) 477 zmumps_c(&lu->id); 478 #else 479 dmumps_c(&lu->id); 480 #endif 481 if (lu->id.INFOG(1) < 0) { 482 SETERRQ1(1,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 483 } 484 } 485 486 /* numerical factorization phase */ 487 if(lu->id.ICNTL(18) == 0) { 488 if (lu->myid == 0) { 489 #if defined(PETSC_USE_COMPLEX) 490 lu->id.a = (mumps_double_complex*)lu->val; 491 #else 492 lu->id.a = lu->val; 493 #endif 494 } 495 } else { 496 #if defined(PETSC_USE_COMPLEX) 497 lu->id.a_loc = (mumps_double_complex*)lu->val; 498 #else 499 lu->id.a_loc = lu->val; 500 #endif 501 } 502 lu->id.job=2; 503 #if defined(PETSC_USE_COMPLEX) 504 zmumps_c(&lu->id); 505 #else 506 dmumps_c(&lu->id); 507 #endif 508 if (lu->id.INFOG(1) < 0) { 509 SETERRQ1(1,"1, Error reported by MUMPS in numerical factorization phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 510 } 511 512 if (lu->myid==0 && lu->id.ICNTL(16) > 0){ 513 SETERRQ1(1," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 514 } 515 516 (*F)->assembled = PETSC_TRUE; 517 lu->matstruc = SAME_NONZERO_PATTERN; 518 lu->CleanUpMUMPS = PETSC_TRUE; 519 PetscFunctionReturn(0); 520 } 521 522 /* Note the Petsc r and c permutations are ignored */ 523 #undef __FUNCT__ 524 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 525 int MatLUFactorSymbolic_AIJMUMPS(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) { 526 Mat B; 527 Mat_MUMPS *lu; 528 int ierr; 529 530 PetscFunctionBegin; 531 532 /* Create the factorization matrix */ 533 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 534 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 535 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 536 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 537 538 B->ops->lufactornumeric = MatFactorNumeric_AIJMUMPS; 539 B->factor = FACTOR_LU; 540 lu = (Mat_MUMPS*)B->spptr; 541 lu->sym = 0; 542 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 543 544 *F = B; 545 PetscFunctionReturn(0); 546 } 547 548 /* Note the Petsc r permutation is ignored */ 549 #undef __FUNCT__ 550 #define __FUNCT__ "MatCholeskyFactorSymbolic_SBAIJMUMPS" 551 int MatCholeskyFactorSymbolic_SBAIJMUMPS(Mat A,IS r,MatFactorInfo *info,Mat *F) { 552 Mat B; 553 Mat_MUMPS *lu; 554 int ierr; 555 556 PetscFunctionBegin; 557 558 /* Create the factorization matrix */ 559 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 560 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 561 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 562 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 563 564 B->ops->choleskyfactornumeric = MatFactorNumeric_AIJMUMPS; 565 B->factor = FACTOR_CHOLESKY; 566 lu = (Mat_MUMPS*)B->spptr; 567 lu->sym = 2; 568 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 569 570 *F = B; 571 PetscFunctionReturn(0); 572 } 573 574 #undef __FUNCT__ 575 #define __FUNCT__ "MatAssemblyEnd_AIJMUMPS" 576 int MatAssemblyEnd_AIJMUMPS(Mat A,MatAssemblyType mode) { 577 int ierr; 578 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 579 580 PetscFunctionBegin; 581 ierr = (*mumps->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 582 583 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 584 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 585 A->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 586 PetscFunctionReturn(0); 587 } 588 589 EXTERN_C_BEGIN 590 #undef __FUNCT__ 591 #define __FUNCT__ "MatConvert_AIJ_AIJMUMPS" 592 int MatConvert_AIJ_AIJMUMPS(Mat A,MatType newtype,Mat *newmat) { 593 int ierr,size; 594 MPI_Comm comm; 595 Mat B=*newmat; 596 Mat_MUMPS *mumps; 597 598 PetscFunctionBegin; 599 if (B != A) { 600 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 601 } 602 603 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 604 ierr = PetscNew(Mat_MUMPS,&mumps);CHKERRQ(ierr); 605 606 mumps->MatDuplicate = A->ops->duplicate; 607 mumps->MatView = A->ops->view; 608 mumps->MatAssemblyEnd = A->ops->assemblyend; 609 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 610 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 611 mumps->MatDestroy = A->ops->destroy; 612 mumps->CleanUpMUMPS = PETSC_FALSE; 613 mumps->isAIJ = PETSC_TRUE; 614 615 B->spptr = (void *)mumps; 616 B->ops->duplicate = MatDuplicate_AIJMUMPS; 617 B->ops->view = MatView_AIJMUMPS; 618 B->ops->assemblyend = MatAssemblyEnd_AIJMUMPS; 619 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 620 B->ops->destroy = MatDestroy_MUMPS; 621 622 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 623 if (size == 1) { 624 ierr = PetscStrallocpy(MATSEQAIJ,&(mumps->basetype));CHKERRQ(ierr); 625 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_aijmumps_C", 626 "MatConvert_AIJ_AIJMUMPS",MatConvert_AIJ_AIJMUMPS);CHKERRQ(ierr); 627 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_aijmumps_seqaij_C", 628 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 629 } else { 630 ierr = PetscStrallocpy(MATMPIAIJ,&(mumps->basetype));CHKERRQ(ierr); 631 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpiaij_aijmumps_C", 632 "MatConvert_AIJ_AIJMUMPS",MatConvert_AIJ_AIJMUMPS);CHKERRQ(ierr); 633 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_aijmumps_mpiaij_C", 634 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 635 } 636 637 PetscLogInfo(0,"Using MUMPS for LU factorization and solves."); 638 ierr = PetscObjectChangeTypeName((PetscObject)B,newtype);CHKERRQ(ierr); 639 *newmat = B; 640 PetscFunctionReturn(0); 641 } 642 EXTERN_C_END 643 644 #undef __FUNCT__ 645 #define __FUNCT__ "MatDuplicate_AIJMUMPS" 646 int MatDuplicate_AIJMUMPS(Mat A, MatDuplicateOption op, Mat *M) { 647 int ierr; 648 Mat_MUMPS *lu=(Mat_MUMPS *)A->spptr; 649 650 PetscFunctionBegin; 651 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 652 ierr = MatConvert_AIJ_AIJMUMPS(*M,MATAIJMUMPS,M);CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 656 /*MC 657 MATAIJMUMPS - MATAIJMUMPS = "aijmumps" - A matrix type providing direct solvers (LU) for distributed 658 and sequential matrices via the external package MUMPS. 659 660 If MUMPS is installed (see the manual for instructions 661 on how to declare the existence of external packages), 662 a matrix type can be constructed which invokes MUMPS solvers. 663 After calling MatCreate(...,A), simply call MatSetType(A,MATAIJMUMPS). 664 This matrix type is only supported for double precision real. 665 666 If created with a single process communicator, this matrix type inherits from MATSEQAIJ. 667 Otherwise, this matrix type inherits from MATMPIAIJ. Hence for single process communicators, 668 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 669 for communicators controlling multiple processes. It is recommended that you call both of 670 the above preallocation routines for simplicity. One can also call MatConvert for an inplace 671 conversion to or from the MATSEQAIJ or MATMPIAIJ type (depending on the communicator size) 672 without data copy. 673 674 Options Database Keys: 675 + -mat_type aijmumps - sets the matrix type to "aijmumps" during a call to MatSetFromOptions() 676 . -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 677 . -mat_mumps_icntl_4 <0,1,2,3,4> - print level 678 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 679 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 680 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 681 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 682 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -sles_view 683 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 684 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 685 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 686 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 687 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 688 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 689 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 690 691 Level: beginner 692 693 .seealso: MATSBAIJMUMPS 694 M*/ 695 696 EXTERN_C_BEGIN 697 #undef __FUNCT__ 698 #define __FUNCT__ "MatCreate_AIJMUMPS" 699 int MatCreate_AIJMUMPS(Mat A) { 700 int ierr,size; 701 MPI_Comm comm; 702 703 PetscFunctionBegin; 704 /* Change type name before calling MatSetType to force proper construction of SeqAIJ or MPIAIJ */ 705 /* and AIJMUMPS types */ 706 ierr = PetscObjectChangeTypeName((PetscObject)A,MATAIJMUMPS);CHKERRQ(ierr); 707 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 708 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 709 if (size == 1) { 710 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 711 } else { 712 ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr); 713 } 714 ierr = MatConvert_AIJ_AIJMUMPS(A,MATAIJMUMPS,&A);CHKERRQ(ierr); 715 PetscFunctionReturn(0); 716 } 717 EXTERN_C_END 718 719 #undef __FUNCT__ 720 #define __FUNCT__ "MatAssemblyEnd_SBAIJMUMPS" 721 int MatAssemblyEnd_SBAIJMUMPS(Mat A,MatAssemblyType mode) { 722 int ierr; 723 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 724 725 PetscFunctionBegin; 726 ierr = (*mumps->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 727 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 728 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 729 A->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 730 PetscFunctionReturn(0); 731 } 732 733 EXTERN_C_BEGIN 734 #undef __FUNCT__ 735 #define __FUNCT__ "MatConvert_SBAIJ_SBAIJMUMPS" 736 int MatConvert_SBAIJ_SBAIJMUMPS(Mat A,MatType newtype,Mat *newmat) { 737 int ierr,size; 738 MPI_Comm comm; 739 Mat B=*newmat; 740 Mat_MUMPS *mumps; 741 742 PetscFunctionBegin; 743 if (B != A) { 744 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 745 } 746 747 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 748 ierr = PetscNew(Mat_MUMPS,&mumps);CHKERRQ(ierr); 749 750 mumps->MatDuplicate = A->ops->duplicate; 751 mumps->MatView = A->ops->view; 752 mumps->MatAssemblyEnd = A->ops->assemblyend; 753 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 754 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 755 mumps->MatDestroy = A->ops->destroy; 756 mumps->CleanUpMUMPS = PETSC_FALSE; 757 mumps->isAIJ = PETSC_FALSE; 758 759 B->spptr = (void *)mumps; 760 B->ops->duplicate = MatDuplicate_SBAIJMUMPS; 761 B->ops->view = MatView_AIJMUMPS; 762 B->ops->assemblyend = MatAssemblyEnd_SBAIJMUMPS; 763 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 764 B->ops->destroy = MatDestroy_MUMPS; 765 766 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 767 if (size == 1) { 768 ierr = PetscStrallocpy(MATSEQSBAIJ,&(mumps->basetype));CHKERRQ(ierr); 769 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_sbaijmumps_C", 770 "MatConvert_SBAIJ_SBAIJMUMPS",MatConvert_SBAIJ_SBAIJMUMPS);CHKERRQ(ierr); 771 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_sbaijmumps_seqsbaij_C", 772 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 773 } else { 774 ierr = PetscStrallocpy(MATMPISBAIJ,&(mumps->basetype));CHKERRQ(ierr); 775 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpisbaij_sbaijmumps_C", 776 "MatConvert_SBAIJ_SBAIJMUMPS",MatConvert_SBAIJ_SBAIJMUMPS);CHKERRQ(ierr); 777 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_sbaijmumps_mpisbaij_C", 778 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 779 } 780 781 PetscLogInfo(0,"Using MUMPS for Cholesky factorization and solves."); 782 ierr = PetscObjectChangeTypeName((PetscObject)B,newtype);CHKERRQ(ierr); 783 *newmat = B; 784 PetscFunctionReturn(0); 785 } 786 EXTERN_C_END 787 788 #undef __FUNCT__ 789 #define __FUNCT__ "MatDuplicate_SBAIJMUMPS" 790 int MatDuplicate_SBAIJMUMPS(Mat A, MatDuplicateOption op, Mat *M) { 791 int ierr; 792 Mat_MUMPS *lu=(Mat_MUMPS *)A->spptr; 793 794 PetscFunctionBegin; 795 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 796 ierr = MatConvert_SBAIJ_SBAIJMUMPS(*M,MATSBAIJMUMPS,M);CHKERRQ(ierr); 797 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_MUMPS));CHKERRQ(ierr); 798 PetscFunctionReturn(0); 799 } 800 801 /*MC 802 MATSBAIJMUMPS - MATSBAIJMUMPS = "sbaijmumps" - A symmetric matrix type providing direct solvers (Cholesky) for 803 distributed and sequential matrices via the external package MUMPS. 804 805 If MUMPS is installed (see the manual for instructions 806 on how to declare the existence of external packages), 807 a matrix type can be constructed which invokes MUMPS solvers. 808 After calling MatCreate(...,A), simply call MatSetType(A,MATSBAIJMUMPS). 809 This matrix type is only supported for double precision real. 810 811 If created with a single process communicator, this matrix type inherits from MATSEQSBAIJ. 812 Otherwise, this matrix type inherits from MATMPISBAIJ. Hence for single process communicators, 813 MatSeqSBAIJSetPreallocation is supported, and similarly MatMPISBAIJSetPreallocation is supported 814 for communicators controlling multiple processes. It is recommended that you call both of 815 the above preallocation routines for simplicity. One can also call MatConvert for an inplace 816 conversion to or from the MATSEQSBAIJ or MATMPISBAIJ type (depending on the communicator size) 817 without data copy. 818 819 Options Database Keys: 820 + -mat_type sbaijmumps - sets the matrix type to "sbaijmumps" during a call to MatSetFromOptions() 821 . -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 822 . -mat_mumps_icntl_4 <0,...,4> - print level 823 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 824 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 825 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 826 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 827 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -sles_view 828 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 829 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 830 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 831 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 832 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 833 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 834 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 835 836 Level: beginner 837 838 .seealso: MATAIJMUMPS 839 M*/ 840 841 EXTERN_C_BEGIN 842 #undef __FUNCT__ 843 #define __FUNCT__ "MatCreate_SBAIJMUMPS" 844 int MatCreate_SBAIJMUMPS(Mat A) { 845 int ierr,size; 846 847 PetscFunctionBegin; 848 /* Change type name before calling MatSetType to force proper construction of SeqSBAIJ or MPISBAIJ */ 849 /* and SBAIJMUMPS types */ 850 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSBAIJMUMPS);CHKERRQ(ierr); 851 ierr = MPI_Comm_size(A->comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 852 if (size == 1) { 853 ierr = MatSetType(A,MATSEQSBAIJ);CHKERRQ(ierr); 854 } else { 855 ierr = MatSetType(A,MATMPISBAIJ);CHKERRQ(ierr); 856 } 857 ierr = MatConvert_SBAIJ_SBAIJMUMPS(A,MATSBAIJMUMPS,&A);CHKERRQ(ierr); 858 PetscFunctionReturn(0); 859 } 860 EXTERN_C_END 861