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 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)(mat->A)->data; 79 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 80 if (mat->bs > 1) SETERRQ1(PETSC_ERR_SUP," bs=%d is not supported yet\n", mat->bs); 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 Mat_AIJ_MUMPS *lua = (Mat_AIJ_MUMPS*)(A)->spptr; 333 int rnz,nnz,ierr,nz,i,M=A->M,*ai,*aj,icntl; 334 PetscTruth valOnly,flg; 335 336 PetscFunctionBegin; 337 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 338 (*F)->ops->solve = MatSolve_AIJ_MUMPS; 339 340 /* Initialize a MUMPS instance */ 341 ierr = MPI_Comm_rank(A->comm, &lu->myid); 342 ierr = MPI_Comm_size(A->comm,&lu->size);CHKERRQ(ierr); 343 lu->id.job = JOB_INIT; 344 ierr = MPI_Comm_dup(A->comm,&(lu->comm_mumps));CHKERRQ(ierr); 345 lu->id.comm_fortran = lu->comm_mumps; 346 347 /* Set mumps options */ 348 ierr = PetscOptionsBegin(A->comm,A->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 349 lu->id.par=1; /* host participates factorizaton and solve */ 350 lu->id.sym=lu->sym; 351 if (lu->sym == 2){ 352 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 353 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 354 } 355 #if defined(PETSC_USE_COMPLEX) 356 zmumps_c(&lu->id); 357 #else 358 dmumps_c(&lu->id); 359 #endif 360 361 if (lu->size == 1){ 362 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 363 } else { 364 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 365 } 366 367 icntl=-1; 368 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 369 if (flg && icntl > 0) { 370 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 371 } else { /* no output */ 372 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 373 lu->id.ICNTL(2) = -1; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 374 lu->id.ICNTL(3) = -1; /* output stream for global information, default=6 */ 375 lu->id.ICNTL(4) = 0; /* level of printing, 0,1,2,3,4, default=2 */ 376 } 377 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); 378 icntl=-1; 379 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 380 if (flg) { 381 if (icntl== 1){ 382 SETERRQ(PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n"); 383 } else { 384 lu->id.ICNTL(7) = icntl; 385 } 386 } 387 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); 388 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); 389 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); 390 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr); 391 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control","None",lu->id.ICNTL(13),&lu->id.ICNTL(13),PETSC_NULL);CHKERRQ(ierr); 392 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): efficiency control","None",lu->id.ICNTL(14),&lu->id.ICNTL(14),PETSC_NULL);CHKERRQ(ierr); 393 ierr = PetscOptionsInt("-mat_mumps_icntl_15","ICNTL(15): efficiency control","None",lu->id.ICNTL(15),&lu->id.ICNTL(15),PETSC_NULL);CHKERRQ(ierr); 394 395 /* 396 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); 397 if (flg){ 398 if (icntl >-1 && icntl <3 ){ 399 if (lu->myid==0) lu->id.ICNTL(16) = icntl; 400 } else { 401 SETERRQ1(PETSC_ERR_SUP,"ICNTL(16)=%d -- not supported\n",icntl); 402 } 403 } 404 */ 405 406 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 407 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); 408 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 409 PetscOptionsEnd(); 410 } 411 412 /* define matrix A */ 413 switch (lu->id.ICNTL(18)){ 414 case 0: /* centralized assembled matrix input (size=1) */ 415 if (!lu->myid) { 416 if (lua->isAIJ){ 417 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 418 nz = aa->nz; 419 ai = aa->i; aj = aa->j; lu->val = aa->a; 420 } else { 421 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 422 nz = aa->s_nz; 423 ai = aa->i; aj = aa->j; lu->val = aa->a; 424 } 425 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization, get irn and jcn */ 426 ierr = PetscMalloc(nz*sizeof(int),&lu->irn);CHKERRQ(ierr); 427 ierr = PetscMalloc(nz*sizeof(int),&lu->jcn);CHKERRQ(ierr); 428 nz = 0; 429 for (i=0; i<M; i++){ 430 rnz = ai[i+1] - ai[i]; 431 while (rnz--) { /* Fortran row/col index! */ 432 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 433 } 434 } 435 } 436 } 437 break; 438 case 3: /* distributed assembled matrix input (size>1) */ 439 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 440 valOnly = PETSC_FALSE; 441 } else { 442 valOnly = PETSC_TRUE; /* only update mat values, not row and col index */ 443 } 444 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 445 break; 446 default: SETERRQ(PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 447 } 448 449 /* analysis phase */ 450 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 451 lu->id.n = M; 452 switch (lu->id.ICNTL(18)){ 453 case 0: /* centralized assembled matrix input */ 454 if (!lu->myid) { 455 lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 456 if (lu->id.ICNTL(6)>1){ 457 #if defined(PETSC_USE_COMPLEX) 458 lu->id.a = (mumps_double_complex*)lu->val; 459 #else 460 lu->id.a = lu->val; 461 #endif 462 } 463 } 464 break; 465 case 3: /* distributed assembled matrix input (size>1) */ 466 lu->id.nz_loc = nnz; 467 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 468 if (lu->id.ICNTL(6)>1) { 469 #if defined(PETSC_USE_COMPLEX) 470 lu->id.a_loc = (mumps_double_complex*)lu->val; 471 #else 472 lu->id.a_loc = lu->val; 473 #endif 474 } 475 break; 476 } 477 lu->id.job=1; 478 #if defined(PETSC_USE_COMPLEX) 479 zmumps_c(&lu->id); 480 #else 481 dmumps_c(&lu->id); 482 #endif 483 if (lu->id.INFOG(1) < 0) { 484 SETERRQ1(1,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 485 } 486 } 487 488 /* numerical factorization phase */ 489 if(lu->id.ICNTL(18) == 0) { 490 if (lu->myid == 0) { 491 #if defined(PETSC_USE_COMPLEX) 492 lu->id.a = (mumps_double_complex*)lu->val; 493 #else 494 lu->id.a = lu->val; 495 #endif 496 } 497 } else { 498 #if defined(PETSC_USE_COMPLEX) 499 lu->id.a_loc = (mumps_double_complex*)lu->val; 500 #else 501 lu->id.a_loc = lu->val; 502 #endif 503 } 504 lu->id.job=2; 505 #if defined(PETSC_USE_COMPLEX) 506 zmumps_c(&lu->id); 507 #else 508 dmumps_c(&lu->id); 509 #endif 510 if (lu->id.INFOG(1) < 0) { 511 SETERRQ1(1,"1, Error reported by MUMPS in numerical factorization phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 512 } 513 514 if (lu->myid==0 && lu->id.ICNTL(16) > 0){ 515 SETERRQ1(1," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 516 } 517 518 (*F)->assembled = PETSC_TRUE; 519 lu->matstruc = SAME_NONZERO_PATTERN; 520 PetscFunctionReturn(0); 521 } 522 523 /* Note the Petsc r and c permutations are ignored */ 524 #undef __FUNCT__ 525 #define __FUNCT__ "MatLUFactorSymbolic_AIJ_MUMPS" 526 int MatLUFactorSymbolic_AIJ_MUMPS(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 527 { 528 Mat B; 529 Mat_AIJ_MUMPS *lu; 530 int ierr; 531 532 PetscFunctionBegin; 533 534 /* Create the factorization matrix */ 535 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 536 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 537 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 538 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 539 540 B->ops->lufactornumeric = MatFactorNumeric_AIJ_MUMPS; 541 B->factor = FACTOR_LU; 542 lu = (Mat_AIJ_MUMPS*)B->spptr; 543 lu->sym = 0; 544 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 545 546 *F = B; 547 PetscFunctionReturn(0); 548 } 549 550 /* Note the Petsc r permutation is ignored */ 551 #undef __FUNCT__ 552 #define __FUNCT__ "MatCholeskyFactorSymbolic_AIJ_MUMPS" 553 int MatCholeskyFactorSymbolic_AIJ_MUMPS(Mat A,IS r,MatFactorInfo *info,Mat *F) 554 { 555 Mat B; 556 Mat_AIJ_MUMPS *lu; 557 int ierr; 558 559 PetscFunctionBegin; 560 561 /* Create the factorization matrix */ 562 ierr = MatCreate(A->comm,A->m,A->n,A->M,A->N,&B);CHKERRQ(ierr); 563 ierr = MatSetType(B,MATAIJMUMPS);CHKERRQ(ierr); 564 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 565 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 566 567 B->ops->choleskyfactornumeric = MatFactorNumeric_AIJ_MUMPS; 568 B->factor = FACTOR_CHOLESKY; 569 lu = (Mat_AIJ_MUMPS *)B->spptr; 570 lu->sym = 2; 571 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 572 573 *F = B; 574 PetscFunctionReturn(0); 575 } 576 577 #undef __FUNCT__ 578 #define __FUNCT__ "MatAssemblyEnd_AIJ_MUMPS" 579 int MatAssemblyEnd_AIJ_MUMPS(Mat A,MatAssemblyType mode) { 580 int ierr; 581 Mat_AIJ_MUMPS *mumps=(Mat_AIJ_MUMPS*)A->spptr; 582 583 PetscFunctionBegin; 584 ierr = (*mumps->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 585 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 586 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 587 A->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJ_MUMPS; 588 A->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJ_MUMPS; 589 PetscFunctionReturn(0); 590 } 591 592 EXTERN_C_BEGIN 593 #undef __FUNCT__ 594 #define __FUNCT__ "MatConvert_Base_MUMPS" 595 int MatConvert_Base_MUMPS(Mat A,MatType newtype,Mat *newmat) { 596 int ierr,size; 597 MPI_Comm comm; 598 Mat B=*newmat; 599 Mat_AIJ_MUMPS *mumps; 600 601 PetscFunctionBegin; 602 if (B != A) { 603 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 604 } 605 606 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 607 ierr = PetscNew(Mat_AIJ_MUMPS,&mumps);CHKERRQ(ierr); 608 609 mumps->MatView = A->ops->view; 610 mumps->MatAssemblyEnd = A->ops->assemblyend; 611 mumps->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 612 mumps->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic; 613 mumps->MatDestroy = A->ops->destroy; 614 mumps->CleanUpMUMPS = PETSC_FALSE; 615 616 B->spptr = (void *)mumps; 617 B->ops->view = MatView_AIJ_MUMPS; 618 B->ops->assemblyend = MatAssemblyEnd_AIJ_MUMPS; 619 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJ_MUMPS; 620 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJ_MUMPS; 621 B->ops->destroy = MatDestroy_AIJ_MUMPS; 622 623 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 624 /* This switch is brutal and should probably be changed, but I didn't want 4 routines. */ 625 if (newtype == MATAIJMUMPS) { 626 mumps->isAIJ = PETSC_TRUE; 627 if (size == 1) { 628 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_aijmumps_C", 629 "MatConvert_Base_MUMPS",MatConvert_Base_MUMPS);CHKERRQ(ierr); 630 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_aijmumps_seqaij_C", 631 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 632 } else { 633 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpiaij_aijmumps_C", 634 "MatConvert_Base_MUMPS",MatConvert_Base_MUMPS);CHKERRQ(ierr); 635 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_aijmumps_mpiaij_C", 636 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 637 } 638 } else { 639 mumps->isAIJ = PETSC_FALSE; 640 if (size == 1) { 641 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_mumps_C", 642 "MatConvert_Base_MUMPS",MatConvert_Base_MUMPS);CHKERRQ(ierr); 643 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mumps_seqsbaij_C", 644 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 645 } else { 646 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpisbaij_mumps_C", 647 "MatConvert_Base_MUMPS",MatConvert_Base_MUMPS);CHKERRQ(ierr); 648 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mumps_mpisbaij_C", 649 "MatConvert_MUMPS_Base",MatConvert_MUMPS_Base);CHKERRQ(ierr); 650 } 651 } 652 653 PetscLogInfo(0,"Using MUMPS for factorization and solves."); 654 ierr = PetscObjectChangeTypeName((PetscObject)B,newtype);CHKERRQ(ierr); 655 *newmat = B; 656 PetscFunctionReturn(0); 657 } 658 EXTERN_C_END 659 660 EXTERN_C_BEGIN 661 #undef __FUNCT__ 662 #define __FUNCT__ "MatCreate_AIJ_MUMPS" 663 int MatCreate_AIJ_MUMPS(Mat A) { 664 int ierr,size; 665 MPI_Comm comm; 666 667 PetscFunctionBegin; 668 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 669 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 670 if (size == 1) { 671 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 672 } else { 673 ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr); 674 } 675 ierr = MatConvert_Base_MUMPS(A,MATAIJMUMPS,&A);CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 EXTERN_C_END 679 680 EXTERN_C_BEGIN 681 #undef __FUNCT__ 682 #define __FUNCT__ "MatCreate_SBAIJ_MUMPS" 683 int MatCreate_SBAIJ_MUMPS(Mat A) { 684 int ierr,size; 685 MPI_Comm comm; 686 687 PetscFunctionBegin; 688 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 689 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);CHKERRQ(ierr); 690 if (size == 1) { 691 ierr = MatSetType(A,MATSEQSBAIJ);CHKERRQ(ierr); 692 } else { 693 ierr = MatSetType(A,MATMPISBAIJ);CHKERRQ(ierr); 694 } 695 ierr = MatConvert_Base_MUMPS(A,MATSBAIJMUMPS,&A);CHKERRQ(ierr); 696 PetscFunctionReturn(0); 697 } 698 EXTERN_C_END 699 700 EXTERN_C_BEGIN 701 #undef __FUNCT__ 702 #define __FUNCT__ "MatLoad_AIJ_MUMPS" 703 int MatLoad_AIJ_MUMPS(PetscViewer viewer,MatType type,Mat *A) { 704 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 705 MPI_Comm comm; 706 707 PetscFunctionBegin; 708 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 709 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 710 if (size == 1) { 711 ierr = PetscFListFind(comm,MatLoadList,MATSEQAIJ,(void(**)(void))&r);CHKERRQ(ierr); 712 } else { 713 ierr = PetscFListFind(comm,MatLoadList,MATMPIAIJ,(void(**)(void))&r);CHKERRQ(ierr); 714 } 715 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 716 PetscFunctionReturn(0); 717 } 718 EXTERN_C_END 719 720 EXTERN_C_BEGIN 721 #undef __FUNCT__ 722 #define __FUNCT__ "MatLoad_SBAIJ_MUMPS" 723 int MatLoad_SBAIJ_MUMPS(PetscViewer viewer,MatType type,Mat *A) { 724 int ierr,size,(*r)(PetscViewer,MatType,Mat*); 725 MPI_Comm comm; 726 727 PetscFunctionBegin; 728 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 729 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 730 if (size == 1) { 731 ierr = PetscFListFind(comm,MatLoadList,MATSEQSBAIJ,(void(**)(void))&r);CHKERRQ(ierr); 732 } else { 733 ierr = PetscFListFind(comm,MatLoadList,MATMPISBAIJ,(void(**)(void))&r);CHKERRQ(ierr); 734 } 735 ierr = (*r)(viewer,type,A);CHKERRQ(ierr); 736 PetscFunctionReturn(0); 737 } 738 EXTERN_C_END 739