1 /*$Id: mumps.c,v 1.10 2001/08/15 15:56:50 bsmith Exp $*/ 2 /* 3 Provides an interface to the MUMPS_4.2_beta sparse solver 4 */ 5 6 #include "src/mat/impls/aij/seq/aij.h" 7 #include "src/mat/impls/aij/mpi/mpiaij.h" 8 #include "src/mat/impls/sbaij/seq/sbaij.h" 9 #include "src/mat/impls/sbaij/mpi/mpisbaij.h" 10 11 EXTERN_C_BEGIN 12 #if defined(PETSC_USE_COMPLEX) 13 #include "zmumps_c.h" 14 #else 15 #include "dmumps_c.h" 16 #endif 17 EXTERN_C_END 18 #define JOB_INIT -1 19 #define JOB_END -2 20 /* macros s.t. indices match MUMPS documentation */ 21 #define ICNTL(I) icntl[(I)-1] 22 #define CNTL(I) cntl[(I)-1] 23 #define INFOG(I) infog[(I)-1] 24 #define RINFOG(I) rinfog[(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 MatType basetype; 38 PetscTruth isAIJ,CleanUpMUMPS; 39 int (*MatView)(Mat,PetscViewer); 40 int (*MatAssemblyEnd)(Mat,MatAssemblyType); 41 int (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 42 int (*MatCholeskyFactorSymbolic)(Mat,IS,MatFactorInfo*,Mat*); 43 int (*MatDestroy)(Mat); 44 } Mat_AIJ_MUMPS; 45 46 /* convert Petsc mpiaij matrix to triples: row[nz], col[nz], val[nz] */ 47 /* 48 input: 49 A - matrix in mpiaij format 50 shift - 0: C style output triple; 1: Fortran style output triple. 51 valOnly - FALSE: spaces are allocated and values are set for the triple 52 TRUE: only the values in v array are updated 53 output: 54 nnz - dim of r, c, and v (number of local nonzero entries of A) 55 r, c, v - row and col index, matrix values (matrix triples) 56 */ 57 int MatConvertToTriples(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 58 { 59 int *ai, *aj, *bi, *bj, rstart,nz, *garray; 60 int ierr,i,j,jj,jB,irow,m=A->m,*ajj,*bjj,countA,countB,colA_start,jcol; 61 int *row,*col; 62 PetscScalar *av, *bv,*val; 63 Mat_AIJ_MUMPS *mumps = (Mat_AIJ_MUMPS *)A->spptr; 64 65 PetscFunctionBegin; 66 67 if (mumps->isAIJ){ 68 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 69 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 70 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 71 nz = aa->nz + bb->nz; 72 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= mat->rstart; 73 garray = mat->garray; 74 av=aa->a; bv=bb->a; 75 76 } else { 77 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 78 if (mat->bs > 1) SETERRQ1(PETSC_ERR_SUP," bs=%d is not supported yet\n", mat->bs); 79 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)(mat->A)->data; 80 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 81 nz = aa->s_nz + bb->nz; 82 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= mat->rstart; 83 garray = mat->garray; 84 av=aa->a; bv=bb->a; 85 } 86 87 if (!valOnly){ 88 ierr = PetscMalloc(nz*sizeof(int),&row);CHKERRQ(ierr); 89 ierr = PetscMalloc(nz*sizeof(int),&col);CHKERRQ(ierr); 90 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 91 *r = row; *c = col; *v = val; 92 } else { 93 row = *r; col = *c; val = *v; 94 } 95 *nnz = nz; 96 97 jj = 0; jB = 0; irow = rstart; 98 for ( i=0; i<m; i++ ) { 99 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 100 countA = ai[i+1] - ai[i]; 101 countB = bi[i+1] - bi[i]; 102 bjj = bj + bi[i]; 103 104 /* get jB, the starting local col index for the 2nd B-part */ 105 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 106 for (j=0; j<countB; j++){ 107 jcol = garray[bjj[j]]; 108 if (jcol > colA_start) { jB = j; break; } 109 if (j==countB-1) jB = countB; 110 } 111 112 /* B-part, smaller col index */ 113 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 114 for (j=0; j<jB; j++){ 115 jcol = garray[bjj[j]]; 116 if (!valOnly){ 117 row[jj] = irow + shift; col[jj] = jcol + shift; 118 } 119 val[jj++] = *bv++; 120 } 121 /* A-part */ 122 for (j=0; j<countA; j++){ 123 if (!valOnly){ 124 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 125 } 126 val[jj++] = *av++; 127 } 128 /* B-part, larger col index */ 129 for (j=jB; j<countB; j++){ 130 if (!valOnly){ 131 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 132 } 133 val[jj++] = *bv++; 134 } 135 irow++; 136 } 137 138 PetscFunctionReturn(0); 139 } 140 141 EXTERN_C_BEGIN 142 #undef __FUNCT__ 143 #define __FUNCT__ "MatConvert_MUMPS_Base" 144 int MatConvert_MUMPS_Base(Mat A,MatType type,Mat *newmat) { 145 /* This routine is only called to convert an unfactored PETSc-MUMPS matrix */ 146 /* to its base PETSc type, so we will ignore 'MatType type'. */ 147 int ierr; 148 Mat B=*newmat; 149 Mat_AIJ_MUMPS *lu=(Mat_AIJ_MUMPS*)A->spptr; 150 151 PetscFunctionBegin; 152 if (B != A) { 153 /* This routine was inherited from SeqAIJ. */ 154 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 155 } else { 156 157 B->ops->view = lu->MatView; 158 B->ops->assemblyend = lu->MatAssemblyEnd; 159 B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic; 160 B->ops->choleskyfactorsymbolic = lu->MatCholeskyFactorSymbolic; 161 B->ops->destroy = lu->MatDestroy; 162 ierr = PetscObjectChangeTypeName((PetscObject)B,lu->basetype);CHKERRQ(ierr); 163 ierr = PetscFree(lu);CHKERRQ(ierr); 164 } 165 *newmat = B; 166 PetscFunctionReturn(0); 167 } 168 EXTERN_C_END 169 170 #undef __FUNCT__ 171 #define __FUNCT__ "MatDestroy_AIJ_MUMPS" 172 int MatDestroy_AIJ_MUMPS(Mat A) 173 { 174 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)A->spptr; 175 int ierr,size=lu->size; 176 177 PetscFunctionBegin; 178 if (lu->CleanUpMUMPS) { 179 /* Terminate instance, deallocate memories */ 180 lu->id.job=JOB_END; 181 #if defined(PETSC_USE_COMPLEX) 182 zmumps_c(&lu->id); 183 #else 184 dmumps_c(&lu->id); 185 #endif 186 if (lu->irn) { 187 ierr = PetscFree(lu->irn);CHKERRQ(ierr); 188 } 189 if (lu->jcn) { 190 ierr = PetscFree(lu->jcn);CHKERRQ(ierr); 191 } 192 if (size>1 && lu->val) { 193 ierr = PetscFree(lu->val);CHKERRQ(ierr); 194 } 195 ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr); 196 } 197 ierr = MatConvert_MUMPS_Base(A,lu->basetype,&A);CHKERRQ(ierr); 198 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "MatFactorInfo_MUMPS" 204 int MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 205 { 206 Mat_AIJ_MUMPS *lu= (Mat_AIJ_MUMPS*)A->spptr; 207 int ierr; 208 209 PetscFunctionBegin; 210 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 211 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 212 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 214 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 215 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 216 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 217 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 218 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 219 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 220 if (lu->myid == 0 && lu->id.ICNTL(11)>0) { 221 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 222 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 223 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 224 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); 225 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 226 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 227 228 } 229 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 230 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 231 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (efficiency control): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 232 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(15) (efficiency control): %d \n",lu->id.ICNTL(15));CHKERRQ(ierr); 233 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 234 235 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 236 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 237 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 238 PetscFunctionReturn(0); 239 } 240 241 #undef __FUNCT__ 242 #define __FUNCT__ "MatView_AIJ_MUMPS" 243 int MatView_AIJ_MUMPS(Mat A,PetscViewer viewer) { 244 int ierr; 245 PetscTruth isascii; 246 PetscViewerFormat format; 247 Mat_AIJ_MUMPS *mumps=(Mat_AIJ_MUMPS*)(A->spptr); 248 249 PetscFunctionBegin; 250 ierr = (*mumps->MatView)(A,viewer);CHKERRQ(ierr); 251 252 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 253 if (isascii) { 254 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 255 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 256 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 257 } 258 } 259 PetscFunctionReturn(0); 260 } 261 262 #undef __FUNCT__ 263 #define __FUNCT__ "MatSolve_AIJ_MUMPS" 264 int MatSolve_AIJ_MUMPS(Mat A,Vec b,Vec x) 265 { 266 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)A->spptr; 267 PetscScalar *array; 268 Vec x_seq; 269 IS iden; 270 VecScatter scat; 271 int ierr; 272 273 PetscFunctionBegin; 274 if (lu->size > 1){ 275 if (!lu->myid){ 276 ierr = VecCreateSeq(PETSC_COMM_SELF,A->N,&x_seq);CHKERRQ(ierr); 277 ierr = ISCreateStride(PETSC_COMM_SELF,A->N,0,1,&iden);CHKERRQ(ierr); 278 } else { 279 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&x_seq);CHKERRQ(ierr); 280 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&iden);CHKERRQ(ierr); 281 } 282 ierr = VecScatterCreate(b,iden,x_seq,iden,&scat);CHKERRQ(ierr); 283 ierr = ISDestroy(iden);CHKERRQ(ierr); 284 285 ierr = VecScatterBegin(b,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 286 ierr = VecScatterEnd(b,x_seq,INSERT_VALUES,SCATTER_FORWARD,scat);CHKERRQ(ierr); 287 if (!lu->myid) {ierr = VecGetArray(x_seq,&array);CHKERRQ(ierr);} 288 } else { /* size == 1 */ 289 ierr = VecCopy(b,x);CHKERRQ(ierr); 290 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 291 } 292 if (!lu->myid) { /* define rhs on the host */ 293 #if defined(PETSC_USE_COMPLEX) 294 lu->id.rhs = (mumps_double_complex*)array; 295 #else 296 lu->id.rhs = array; 297 #endif 298 } 299 300 /* solve phase */ 301 lu->id.job=3; 302 #if defined(PETSC_USE_COMPLEX) 303 zmumps_c(&lu->id); 304 #else 305 dmumps_c(&lu->id); 306 #endif 307 if (lu->id.INFOG(1) < 0) { 308 SETERRQ1(1,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 309 } 310 311 /* convert mumps solution x_seq to petsc mpi x */ 312 if (lu->size > 1) { 313 if (!lu->myid){ 314 ierr = VecRestoreArray(x_seq,&array);CHKERRQ(ierr); 315 } 316 ierr = VecScatterBegin(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 317 ierr = VecScatterEnd(x_seq,x,INSERT_VALUES,SCATTER_REVERSE,scat);CHKERRQ(ierr); 318 ierr = VecScatterDestroy(scat);CHKERRQ(ierr); 319 ierr = VecDestroy(x_seq);CHKERRQ(ierr); 320 } else { 321 ierr = VecRestoreArray(x,&array);CHKERRQ(ierr); 322 } 323 324 PetscFunctionReturn(0); 325 } 326 327 #undef __FUNCT__ 328 #define __FUNCT__ "MatFactorNumeric_MPIAIJ_MUMPS" 329 int MatFactorNumeric_AIJ_MUMPS(Mat A,Mat *F) 330 { 331 Mat_AIJ_MUMPS *lu = (Mat_AIJ_MUMPS*)(*F)->spptr; 332 int rnz,nnz,ierr,nz,i,M=A->M,*ai,*aj,icntl; 333 PetscTruth valOnly,flg; 334 335 PetscFunctionBegin; 336 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 337 (*F)->ops->solve = MatSolve_AIJ_MUMPS; 338 339 /* Initialize a MUMPS instance */ 340 ierr = MPI_Comm_rank(A->comm, &lu->myid); 341 ierr = MPI_Comm_size(A->comm,&lu->size);CHKERRQ(ierr); 342 lu->id.job = JOB_INIT; 343 ierr = MPI_Comm_dup(A->comm,&(lu->comm_mumps));CHKERRQ(ierr); 344 lu->id.comm_fortran = lu->comm_mumps; 345 346 /* Set mumps options */ 347 ierr = PetscOptionsBegin(A->comm,A->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 348 lu->id.par=1; /* host participates factorizaton and solve */ 349 lu->id.sym=lu->sym; 350 if (lu->sym == 2){ 351 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 352 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 353 } 354 #if defined(PETSC_USE_COMPLEX) 355 zmumps_c(&lu->id); 356 #else 357 dmumps_c(&lu->id); 358 #endif 359 360 if (lu->size == 1){ 361 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 362 } else { 363 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 364 } 365 366 icntl=-1; 367 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 368 if (flg && icntl > 0) { 369 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 370 } else { /* no output */ 371 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 372 lu->id.ICNTL(2) = -1; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 373 lu->id.ICNTL(3) = -1; /* output stream for global information, default=6 */ 374 lu->id.ICNTL(4) = 0; /* level of printing, 0,1,2,3,4, default=2 */ 375 } 376 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); 377 icntl=-1; 378 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 379 if (flg) { 380 if (icntl== 1){ 381 SETERRQ(PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n"); 382 } else { 383 lu->id.ICNTL(7) = icntl; 384 } 385 } 386 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); 387 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); 388 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); 389 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr); 390 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control","None",lu->id.ICNTL(13),&lu->id.ICNTL(13),PETSC_NULL);CHKERRQ(ierr); 391 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): efficiency control","None",lu->id.ICNTL(14),&lu->id.ICNTL(14),PETSC_NULL);CHKERRQ(ierr); 392 ierr = PetscOptionsInt("-mat_mumps_icntl_15","ICNTL(15): efficiency control","None",lu->id.ICNTL(15),&lu->id.ICNTL(15),PETSC_NULL);CHKERRQ(ierr); 393 394 /* 395 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); 396 if (flg){ 397 if (icntl >-1 && icntl <3 ){ 398 if (lu->myid==0) lu->id.ICNTL(16) = icntl; 399 } else { 400 SETERRQ1(PETSC_ERR_SUP,"ICNTL(16)=%d -- not supported\n",icntl); 401 } 402 } 403 */ 404 405 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 406 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); 407 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 408 PetscOptionsEnd(); 409 } 410 411 /* define matrix A */ 412 switch (lu->id.ICNTL(18)){ 413 case 0: /* centralized assembled matrix input (size=1) */ 414 if (!lu->myid) { 415 if (lu->isAIJ){ 416 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 417 nz = aa->nz; 418 ai = aa->i; aj = aa->j; lu->val = aa->a; 419 } else { 420 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 421 nz = aa->s_nz; 422 ai = aa->i; aj = aa->j; lu->val = aa->a; 423 } 424 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization, get irn and jcn */ 425 ierr = PetscMalloc(nz*sizeof(int),&lu->irn);CHKERRQ(ierr); 426 ierr = PetscMalloc(nz*sizeof(int),&lu->jcn);CHKERRQ(ierr); 427 nz = 0; 428 for (i=0; i<M; i++){ 429 rnz = ai[i+1] - ai[i]; 430 while (rnz--) { /* Fortran row/col index! */ 431 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 432 } 433 } 434 } 435 } 436 break; 437 case 3: /* distributed assembled matrix input (size>1) */ 438 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 439 valOnly = PETSC_FALSE; 440 } else { 441 valOnly = PETSC_TRUE; /* only update mat values, not row and col index */ 442 } 443 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 444 break; 445 default: SETERRQ(PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 446 } 447 448 /* analysis phase */ 449 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 450 lu->id.n = M; 451 switch (lu->id.ICNTL(18)){ 452 case 0: /* centralized assembled matrix input */ 453 if (!lu->myid) { 454 lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 455 if (lu->id.ICNTL(6)>1){ 456 #if defined(PETSC_USE_COMPLEX) 457 lu->id.a = (mumps_double_complex*)lu->val; 458 #else 459 lu->id.a = lu->val; 460 #endif 461 } 462 } 463 break; 464 case 3: /* distributed assembled matrix input (size>1) */ 465 lu->id.nz_loc = nnz; 466 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 467 if (lu->id.ICNTL(6)>1) { 468 #if defined(PETSC_USE_COMPLEX) 469 lu->id.a_loc = (mumps_double_complex*)lu->val; 470 #else 471 lu->id.a_loc = lu->val; 472 #endif 473 } 474 break; 475 } 476 lu->id.job=1; 477 #if defined(PETSC_USE_COMPLEX) 478 zmumps_c(&lu->id); 479 #else 480 dmumps_c(&lu->id); 481 #endif 482 if (lu->id.INFOG(1) < 0) { 483 SETERRQ1(1,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 484 } 485 } 486 487 /* numerical factorization phase */ 488 if(lu->id.ICNTL(18) == 0) { 489 if (lu->myid == 0) { 490 #if defined(PETSC_USE_COMPLEX) 491 lu->id.a = (mumps_double_complex*)lu->val; 492 #else 493 lu->id.a = lu->val; 494 #endif 495 } 496 } else { 497 #if defined(PETSC_USE_COMPLEX) 498 lu->id.a_loc = (mumps_double_complex*)lu->val; 499 #else 500 lu->id.a_loc = lu->val; 501 #endif 502 } 503 lu->id.job=2; 504 #if defined(PETSC_USE_COMPLEX) 505 zmumps_c(&lu->id); 506 #else 507 dmumps_c(&lu->id); 508 #endif 509 if (lu->id.INFOG(1) < 0) { 510 SETERRQ1(1,"1, Error reported by MUMPS in numerical factorization phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 511 } 512 513 if (lu->myid==0 && lu->id.ICNTL(16) > 0){ 514 SETERRQ1(1," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 515 } 516 517 (*F)->assembled = PETSC_TRUE; 518 lu->matstruc = SAME_NONZERO_PATTERN; 519 PetscFunctionReturn(0); 520 } 521 522 /* Note the Petsc r and c permutations are ignored */ 523 #undef __FUNCT__ 524 #define __FUNCT__ "MatLUFactorSymbolic_AIJ_MUMPS" 525 int MatLUFactorSymbolic_AIJ_MUMPS(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 526 { 527 Mat B; 528 Mat_AIJ_MUMPS *lu; 529 int ierr; 530 531 PetscFunctionBegin; 532 533 /* Create the factorization matrix */ 534 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 535 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 536 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 537 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 538 539 B->ops->lufactornumeric = MatFactorNumeric_AIJ_MUMPS; 540 B->factor = FACTOR_LU; 541 lu = (Mat_AIJ_MUMPS*)B->spptr; 542 lu->sym = 0; 543 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 544 545 *F = B; 546 PetscFunctionReturn(0); 547 } 548 549 /* Note the Petsc r permutation is ignored */ 550 #undef __FUNCT__ 551 #define __FUNCT__ "MatCholeskyFactorSymbolic_AIJ_MUMPS" 552 int MatCholeskyFactorSymbolic_AIJ_MUMPS(Mat A,IS r,MatFactorInfo *info,Mat *F) 553 { 554 Mat B; 555 Mat_AIJ_MUMPS *lu; 556 int ierr; 557 558 PetscFunctionBegin; 559 560 /* Create the factorization matrix */ 561 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 562 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 563 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 564 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 565 566 B->ops->choleskyfactornumeric = MatFactorNumeric_AIJ_MUMPS; 567 B->factor = FACTOR_CHOLESKY; 568 lu = (Mat_AIJ_MUMPS *)B->spptr; 569 lu->sym = 2; 570 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 571 572 *F = B; 573 PetscFunctionReturn(0); 574 } 575 576 #undef __FUNCT__ 577 #define __FUNCT__ "MatAssemblyEnd_AIJ_MUMPS" 578 int MatAssemblyEnd_AIJ_MUMPS(Mat A,MatAssemblyType mode) { 579 int ierr; 580 Mat_AIJ_MUMPS *mumps=(Mat_AIJ_MUMPS*)A->spptr; 581 582 PetscFunctionBegin; 583 ierr = (*mumps->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 584 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 585 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 586 A->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJ_MUMPS; 587 A->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJ_MUMPS; 588 PetscFunctionReturn(0); 589 } 590 591 EXTERN_C_BEGIN 592 #undef __FUNCT__ 593 #define __FUNCT__ "MatConvert_Base_MUMPS" 594 int MatConvert_Base_MUMPS(Mat A,MatType newtype,Mat *newmat) { 595 int ierr,size; 596 MPI_Comm comm; 597 Mat B=*newmat; 598 Mat_AIJ_MUMPS *mumps; 599 600 PetscFunctionBegin; 601 if (B != A) { 602 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 603 } 604 605 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 606 ierr = PetscNew(Mat_AIJ_MUMPS,&mumps);CHKERRQ(ierr); 607 608 mumps->MatView = A->ops->view; 609 mumps->MatAssemblyEnd = A->ops->assemblyend; 610 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 611 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 612 mumps->MatDestroy = A->ops->destroy; 613 mumps->CleanUpMUMPS = PETSC_FALSE; 614 615 A->spptr = (void *)mumps; 616 A->ops->view = MatView_AIJ_MUMPS; 617 A->ops->assemblyend = MatAssemblyEnd_AIJ_MUMPS; 618 A->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJ_MUMPS; 619 A->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJ_MUMPS; 620 A->ops->destroy = MatDestroy_AIJ_MUMPS; 621 622 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 623 if (newtype == MATAIJMUMPS) { /* This is brutal and should probably be changed, but I didn't want 4 routines. */ 624 if (size == 1) { 625 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_aijmumps_C", 626 "MatConvert_Base_MUMPS",MatConvert_Base_MUMPS);CHKERRQ(ierr); 627 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_aijmumps_seqaij_C", 628 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 629 } else { 630 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpiaij_aijmumps_C", 631 "MatConvert_Base_MUMPS",MatConvert_Base_MUMPS);CHKERRQ(ierr); 632 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_aijmumps_mpiaij_C", 633 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 634 } 635 } else { 636 if (size == 1) { 637 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_mumps_C", 638 "MatConvert_Base_MUMPS",MatConvert_Base_MUMPS);CHKERRQ(ierr); 639 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mumps_seqsbaij_C", 640 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 641 } else { 642 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpisbaij_mumps_C", 643 "MatConvert_Base_MUMPS",MatConvert_Base_MUMPS);CHKERRQ(ierr); 644 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mumps_mpisbaij_C", 645 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 646 } 647 } 648 649 PetscLogInfo(0,"Using MUMPS for factorization and solves."); 650 ierr = PetscObjectChangeTypeName((PetscObject)B,newtype);CHKERRQ(ierr); 651 *newmat = B; 652 PetscFunctionReturn(0); 653 } 654 EXTERN_C_END 655 656 EXTERN_C_BEGIN 657 #undef __FUNCT__ 658 #define __FUNCT__ "MatCreate_AIJ_MUMPS" 659 int MatCreate_AIJ_MUMPS(Mat A) { 660 int ierr,size; 661 MPI_Comm comm; 662 663 PetscFunctionBegin; 664 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 665 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 666 if (size == 1) { 667 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 668 } else { 669 ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr); 670 } 671 ierr = MatConvert_Base_MUMPS(A,MATAIJMUMPS,&A);CHKERRQ(ierr); 672 PetscFunctionReturn(0); 673 } 674 EXTERN_C_END 675 676 EXTERN_C_BEGIN 677 #undef __FUNCT__ 678 #define __FUNCT__ "MatCreate_SBAIJ_MUMPS" 679 int MatCreate_SBAIJ_MUMPS(Mat A) { 680 int ierr,size; 681 MPI_Comm comm; 682 683 PetscFunctionBegin; 684 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 685 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 686 if (size == 1) { 687 ierr = MatSetType(A,MATSEQSBAIJ);CHKERRQ(ierr); 688 } else { 689 ierr = MatSetType(A,MATMPISBAIJ);CHKERRQ(ierr); 690 } 691 ierr = MatConvert_Base_MUMPS(A,MATSBAIJMUMPS,&A);CHKERRQ(ierr); 692 PetscFunctionReturn(0); 693 } 694 EXTERN_C_END 695 696 EXTERN_C_BEGIN 697 #undef __FUNCT__ 698 #define __FUNCT__ "MatLoad_AIJ_MUMPS" 699 int MatLoad_AIJ_MUMPS(PetscViewer viewer,MatType type,Mat *A) { 700 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 701 MPI_Comm comm; 702 703 PetscFunctionBegin; 704 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 705 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 706 if (size == 1) { 707 ierr = PetscFListFind(comm,MatLoadList,MATSEQAIJ,(void(**)(void))&r);CHKERRQ(ierr); 708 } else { 709 ierr = PetscFListFind(comm,MatLoadList,MATMPIAIJ,(void(**)(void))&r);CHKERRQ(ierr); 710 } 711 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 712 PetscFunctionReturn(0); 713 } 714 EXTERN_C_END 715 716 EXTERN_C_BEGIN 717 #undef __FUNCT__ 718 #define __FUNCT__ "MatLoad_SBAIJ_MUMPS" 719 int MatLoad_SBAIJ_MUMPS(PetscViewer viewer,MatType type,Mat *A) { 720 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 721 MPI_Comm comm; 722 723 PetscFunctionBegin; 724 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 725 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 726 if (size == 1) { 727 ierr = PetscFListFind(comm,MatLoadList,MATSEQSBAIJ,(void(**)(void))&r);CHKERRQ(ierr); 728 } else { 729 ierr = PetscFListFind(comm,MatLoadList,MATMPISBAIJ,(void(**)(void))&r);CHKERRQ(ierr); 730 } 731 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 732 PetscFunctionReturn(0); 733 } 734 EXTERN_C_END 735