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