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