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