1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the MUMPS sparse solver 5 */ 6 #include "../src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 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 #include "../src/mat/impls/baij/seq/baij.h" 11 #include "../src/mat/impls/baij/mpi/mpibaij.h" 12 13 EXTERN_C_BEGIN 14 #if defined(PETSC_USE_COMPLEX) 15 #include "zmumps_c.h" 16 #else 17 #include "dmumps_c.h" 18 #endif 19 EXTERN_C_END 20 #define JOB_INIT -1 21 #define JOB_END -2 22 /* macros s.t. indices match MUMPS documentation */ 23 #define ICNTL(I) icntl[(I)-1] 24 #define CNTL(I) cntl[(I)-1] 25 #define INFOG(I) infog[(I)-1] 26 #define INFO(I) info[(I)-1] 27 #define RINFOG(I) rinfog[(I)-1] 28 #define RINFO(I) rinfo[(I)-1] 29 30 typedef struct { 31 #if defined(PETSC_USE_COMPLEX) 32 ZMUMPS_STRUC_C id; 33 #else 34 DMUMPS_STRUC_C id; 35 #endif 36 MatStructure matstruc; 37 PetscMPIInt myid,size; 38 PetscInt *irn,*jcn,sym,nSolve; 39 PetscScalar *val; 40 MPI_Comm comm_mumps; 41 VecScatter scat_rhs, scat_sol; 42 PetscTruth isAIJ,CleanUpMUMPS,mumpsview; 43 Vec b_seq,x_seq; 44 PetscErrorCode (*MatDestroy)(Mat); 45 } Mat_MUMPS; 46 47 EXTERN PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*); 48 49 /* convert Petsc mpiaij matrix to triples: row[nz], col[nz], val[nz] */ 50 /* 51 input: 52 A - matrix in mpiaij or mpisbaij (bs=1) format 53 shift - 0: C style output triple; 1: Fortran style output triple. 54 valOnly - FALSE: spaces are allocated and values are set for the triple 55 TRUE: only the values in v array are updated 56 output: 57 nnz - dim of r, c, and v (number of local nonzero entries of A) 58 r, c, v - row and col index, matrix values (matrix triples) 59 */ 60 PetscErrorCode MatConvertToTriples(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 61 { 62 PetscInt *ai, *aj, *bi, *bj, rstart,nz, *garray; 63 PetscErrorCode ierr; 64 PetscInt i,j,jj,jB,irow,m=A->rmap->n,*ajj,*bjj,countA,countB,colA_start,jcol; 65 PetscInt *row,*col; 66 PetscScalar *av, *bv,*val; 67 PetscTruth isAIJ; 68 69 PetscFunctionBegin; 70 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isAIJ);CHKERRQ(ierr); 71 if (isAIJ){ 72 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 73 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 74 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 75 nz = aa->nz + bb->nz; 76 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 77 garray = mat->garray; 78 av=aa->a; bv=bb->a; 79 80 } else { 81 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 82 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)(mat->A)->data; 83 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 84 if (A->rmap->bs > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP," bs=%d is not supported yet\n", A->rmap->bs); 85 nz = aa->nz + bb->nz; 86 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 87 garray = mat->garray; 88 av=aa->a; bv=bb->a; 89 } 90 91 if (!valOnly){ 92 ierr = PetscMalloc(nz*sizeof(PetscInt) ,&row);CHKERRQ(ierr); 93 ierr = PetscMalloc(nz*sizeof(PetscInt),&col);CHKERRQ(ierr); 94 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 95 *r = row; *c = col; *v = val; 96 } else { 97 row = *r; col = *c; val = *v; 98 } 99 *nnz = nz; 100 101 jj = 0; irow = rstart; 102 for ( i=0; i<m; i++ ) { 103 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 104 countA = ai[i+1] - ai[i]; 105 countB = bi[i+1] - bi[i]; 106 bjj = bj + bi[i]; 107 108 /* get jB, the starting local col index for the 2nd B-part */ 109 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 110 j=-1; 111 do { 112 j++; 113 if (j == countB) break; 114 jcol = garray[bjj[j]]; 115 } while (jcol < colA_start); 116 jB = j; 117 118 /* B-part, smaller col index */ 119 colA_start = rstart + ajj[0]; /* the smallest col index for A */ 120 for (j=0; j<jB; j++){ 121 jcol = garray[bjj[j]]; 122 if (!valOnly){ 123 row[jj] = irow + shift; col[jj] = jcol + shift; 124 125 } 126 val[jj++] = *bv++; 127 } 128 /* A-part */ 129 for (j=0; j<countA; j++){ 130 if (!valOnly){ 131 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 132 } 133 val[jj++] = *av++; 134 } 135 /* B-part, larger col index */ 136 for (j=jB; j<countB; j++){ 137 if (!valOnly){ 138 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 139 } 140 val[jj++] = *bv++; 141 } 142 irow++; 143 } 144 PetscFunctionReturn(0); 145 } 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "MatDestroy_MUMPS" 149 PetscErrorCode MatDestroy_MUMPS(Mat A) 150 { 151 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 152 PetscErrorCode ierr; 153 PetscMPIInt size=lu->size; 154 155 PetscFunctionBegin; 156 if (lu->CleanUpMUMPS) { 157 /* Terminate instance, deallocate memories */ 158 if (size > 1){ 159 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 160 ierr = VecScatterDestroy(lu->scat_rhs);CHKERRQ(ierr); 161 ierr = VecDestroy(lu->b_seq);CHKERRQ(ierr); 162 if (lu->nSolve && lu->scat_sol){ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);} 163 if (lu->nSolve && lu->x_seq){ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);} 164 ierr = PetscFree(lu->val);CHKERRQ(ierr); 165 } 166 if( size == 1 && (A)->factortype == MAT_FACTOR_CHOLESKY && lu->isAIJ) { 167 ierr = PetscFree(lu->val);CHKERRQ(ierr); 168 } 169 lu->id.job=JOB_END; 170 #if defined(PETSC_USE_COMPLEX) 171 zmumps_c(&lu->id); 172 #else 173 dmumps_c(&lu->id); 174 #endif 175 ierr = PetscFree(lu->irn);CHKERRQ(ierr); 176 ierr = PetscFree(lu->jcn);CHKERRQ(ierr); 177 ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr); 178 } 179 /* clear composed functions */ 180 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatFactorGetSolverPackage_C","",PETSC_NULL);CHKERRQ(ierr); 181 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMumpsSetIcntl_C","",PETSC_NULL);CHKERRQ(ierr); 182 ierr = (lu->MatDestroy)(A);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "MatSolve_MUMPS" 188 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 189 { 190 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 191 PetscScalar *array; 192 Vec x_seq; 193 IS is_iden,is_petsc; 194 PetscErrorCode ierr; 195 PetscInt i; 196 197 PetscFunctionBegin; 198 lu->id.nrhs = 1; 199 x_seq = lu->b_seq; 200 if (lu->size > 1){ 201 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 202 ierr = VecScatterBegin(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 203 ierr = VecScatterEnd(lu->scat_rhs,b,x_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 204 if (!lu->myid) {ierr = VecGetArray(x_seq,&array);CHKERRQ(ierr);} 205 } else { /* size == 1 */ 206 ierr = VecCopy(b,x);CHKERRQ(ierr); 207 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 208 } 209 if (!lu->myid) { /* define rhs on the host */ 210 lu->id.nrhs = 1; 211 #if defined(PETSC_USE_COMPLEX) 212 lu->id.rhs = (mumps_double_complex*)array; 213 #else 214 lu->id.rhs = array; 215 #endif 216 } 217 if (lu->size == 1){ 218 ierr = VecRestoreArray(x,&array);CHKERRQ(ierr); 219 } else if (!lu->myid){ 220 ierr = VecRestoreArray(x_seq,&array);CHKERRQ(ierr); 221 } 222 223 if (lu->size > 1){ 224 /* distributed solution */ 225 lu->id.ICNTL(21) = 1; 226 if (!lu->nSolve){ 227 /* Create x_seq=sol_loc for repeated use */ 228 PetscInt lsol_loc; 229 PetscScalar *sol_loc; 230 lsol_loc = lu->id.INFO(23); /* length of sol_loc */ 231 ierr = PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&lu->id.isol_loc);CHKERRQ(ierr); 232 lu->id.lsol_loc = lsol_loc; 233 #if defined(PETSC_USE_COMPLEX) 234 lu->id.sol_loc = (mumps_double_complex*)sol_loc; 235 #else 236 lu->id.sol_loc = sol_loc; 237 #endif 238 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,lsol_loc,sol_loc,&lu->x_seq);CHKERRQ(ierr); 239 } 240 } 241 242 /* solve phase */ 243 /*-------------*/ 244 lu->id.job = 3; 245 #if defined(PETSC_USE_COMPLEX) 246 zmumps_c(&lu->id); 247 #else 248 dmumps_c(&lu->id); 249 #endif 250 if (lu->id.INFOG(1) < 0) { 251 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 252 } 253 254 if (lu->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 255 if (!lu->nSolve){ /* create scatter scat_sol */ 256 ierr = ISCreateStride(PETSC_COMM_SELF,lu->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 257 for (i=0; i<lu->id.lsol_loc; i++){ 258 lu->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 259 } 260 ierr = ISCreateGeneral(PETSC_COMM_SELF,lu->id.lsol_loc,lu->id.isol_loc,&is_petsc);CHKERRQ(ierr); /* to */ 261 ierr = VecScatterCreate(lu->x_seq,is_iden,x,is_petsc,&lu->scat_sol);CHKERRQ(ierr); 262 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 263 ierr = ISDestroy(is_petsc);CHKERRQ(ierr); 264 } 265 ierr = VecScatterBegin(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 266 ierr = VecScatterEnd(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 267 } 268 lu->nSolve++; 269 PetscFunctionReturn(0); 270 } 271 272 #if !defined(PETSC_USE_COMPLEX) 273 /* 274 input: 275 F: numeric factor 276 output: 277 nneg: total number of negative pivots 278 nzero: 0 279 npos: (global dimension of F) - nneg 280 */ 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 284 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 285 { 286 Mat_MUMPS *lu =(Mat_MUMPS*)F->spptr; 287 PetscErrorCode ierr; 288 PetscMPIInt size; 289 290 PetscFunctionBegin; 291 ierr = MPI_Comm_size(((PetscObject)F)->comm,&size);CHKERRQ(ierr); 292 /* 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 */ 293 if (size > 1 && lu->id.ICNTL(13) != 1){ 294 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"ICNTL(13)=%d. -mat_mumps_icntl_13 must be set as 1 for correct global matrix inertia\n",lu->id.INFOG(13)); 295 } 296 if (nneg){ 297 if (!lu->myid){ 298 *nneg = lu->id.INFOG(12); 299 } 300 ierr = MPI_Bcast(nneg,1,MPI_INT,0,lu->comm_mumps);CHKERRQ(ierr); 301 } 302 if (nzero) *nzero = 0; 303 if (npos) *npos = F->rmap->N - (*nneg); 304 PetscFunctionReturn(0); 305 } 306 #endif /* !defined(PETSC_USE_COMPLEX) */ 307 308 #undef __FUNCT__ 309 #define __FUNCT__ "MatFactorNumeric_MUMPS" 310 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 311 { 312 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 313 Mat newMat; 314 PetscErrorCode ierr; 315 PetscInt rnz,nnz,nz=0,i,M=A->rmap->N,*ai,*aj,*adiag,jidx; 316 PetscScalar *av; 317 PetscTruth valOnly; 318 Mat F_diag; 319 IS is_iden; 320 Vec b; 321 PetscTruth isSeqAIJ,isSeqSBAIJ,isMPIAIJ; 322 323 PetscFunctionBegin; 324 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 325 ierr = PetscTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 326 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 327 328 /* define matrix A */ 329 switch (lu->id.ICNTL(18)){ 330 case 0: /* centralized assembled matrix input (size=1) */ 331 if (!lu->myid) { 332 if (isSeqAIJ){ 333 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 334 nz = aa->nz; 335 ai = aa->i; aj = aa->j; adiag = aa->diag; av = aa->a; 336 } else if (isSeqSBAIJ) { 337 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 338 nz = aa->nz; 339 ai = aa->i; aj = aa->j; av = aa->a; 340 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No mumps factorization for this matrix type"); 341 342 if ((F)->factortype == MAT_FACTOR_CHOLESKY && isSeqAIJ) { 343 nz = 0; 344 for (i=0; i<M; i++){ 345 rnz = ai[i+1] - adiag[i]; 346 jidx = adiag[i]; 347 while (rnz--) { /* Fortran row/col index! */ 348 lu->val[nz] = av[jidx]; jidx++; nz++; 349 } 350 } 351 } else { 352 lu->val = av; 353 } 354 } 355 break; 356 case 3: /* distributed assembled matrix input (size>1) */ 357 valOnly = PETSC_TRUE; /* only update mat values, not row and col index */ 358 359 if(((F)->factortype == MAT_FACTOR_CHOLESKY) && isMPIAIJ) { 360 /* Create an SBAIJ matrix and use this matrix to set the lu values */ 361 ierr = MatConvert(A,MATMPISBAIJ,MAT_INITIAL_MATRIX,&newMat);CHKERRQ(ierr); 362 ierr = MatConvertToTriples(newMat,1,valOnly,&nnz,&lu->irn , &lu->jcn, &lu->val);CHKERRQ(ierr); 363 ierr = MatDestroy(newMat);CHKERRQ(ierr); 364 } 365 else { 366 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 367 } 368 break; 369 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 370 } 371 372 /* analysis phase */ 373 /*----------------*/ 374 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 375 lu->id.job = 1; 376 377 lu->id.n = M; 378 switch (lu->id.ICNTL(18)){ 379 case 0: /* centralized assembled matrix input */ 380 if (!lu->myid) { 381 lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 382 if (lu->id.ICNTL(6)>1){ 383 #if defined(PETSC_USE_COMPLEX) 384 lu->id.a = (mumps_double_complex*)lu->val; 385 #else 386 lu->id.a = lu->val; 387 #endif 388 } 389 } 390 break; 391 case 3: /* distributed assembled matrix input (size>1) */ 392 lu->id.nz_loc = nnz; 393 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 394 if (lu->id.ICNTL(6)>1) { 395 #if defined(PETSC_USE_COMPLEX) 396 lu->id.a_loc = (mumps_double_complex*)lu->val; 397 #else 398 lu->id.a_loc = lu->val; 399 #endif 400 } 401 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 402 if (!lu->myid){ 403 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 404 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 405 } else { 406 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 407 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 408 } 409 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 410 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 411 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 412 413 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 414 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 415 ierr = VecDestroy(b);CHKERRQ(ierr); 416 break; 417 } 418 #if defined(PETSC_USE_COMPLEX) 419 zmumps_c(&lu->id); 420 #else 421 dmumps_c(&lu->id); 422 #endif 423 if (lu->id.INFOG(1) < 0) { 424 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 425 } 426 } 427 428 /* numerical factorization phase */ 429 /*-------------------------------*/ 430 lu->id.job = 2; 431 if(!lu->id.ICNTL(18)) { 432 if (!lu->myid) { 433 #if defined(PETSC_USE_COMPLEX) 434 lu->id.a = (mumps_double_complex*)lu->val; 435 #else 436 lu->id.a = lu->val; 437 #endif 438 } 439 } else { 440 #if defined(PETSC_USE_COMPLEX) 441 lu->id.a_loc = (mumps_double_complex*)lu->val; 442 #else 443 lu->id.a_loc = lu->val; 444 #endif 445 } 446 #if defined(PETSC_USE_COMPLEX) 447 zmumps_c(&lu->id); 448 #else 449 dmumps_c(&lu->id); 450 #endif 451 if (lu->id.INFOG(1) < 0) { 452 if (lu->id.INFO(1) == -13) { 453 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d megabytes\n",lu->id.INFO(2)); 454 } else { 455 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: INFO(1)=%d, INFO(2)=%d\n",lu->id.INFO(1),lu->id.INFO(2)); 456 } 457 } 458 459 if (!lu->myid && lu->id.ICNTL(16) > 0){ 460 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 461 } 462 463 if (lu->size > 1){ 464 if(isMPIAIJ) { 465 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 466 } else { 467 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 468 } 469 F_diag->assembled = PETSC_TRUE; 470 if (lu->nSolve){ 471 ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr); 472 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 473 ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr); 474 } 475 } 476 (F)->assembled = PETSC_TRUE; 477 lu->matstruc = SAME_NONZERO_PATTERN; 478 lu->CleanUpMUMPS = PETSC_TRUE; 479 lu->nSolve = 0; 480 PetscFunctionReturn(0); 481 } 482 483 #undef __FUNCT__ 484 #define __FUNCT__ "PetscSetMUMPSOptions" 485 PetscErrorCode PetscSetMUMPSOptions(Mat F, Mat A) 486 { 487 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 488 PetscErrorCode ierr; 489 PetscInt icntl; 490 PetscTruth flg; 491 492 PetscFunctionBegin; 493 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 494 ierr = PetscOptionsTruth("-mat_mumps_view","View MUMPS parameters","None",lu->mumpsview,&lu->mumpsview,PETSC_NULL);CHKERRQ(ierr); 495 if (lu->size == 1){ 496 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 497 } else { 498 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 499 } 500 501 icntl=-1; 502 lu->id.ICNTL(4) = 0; /* level of printing; overwrite mumps default ICNTL(4)=2 */ 503 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 504 if ((flg && icntl > 0) || PetscLogPrintInfo) { 505 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 506 } else { /* no output */ 507 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 508 lu->id.ICNTL(2) = 0; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 509 lu->id.ICNTL(3) = 0; /* output stream for global information, default=6 */ 510 } 511 ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): column permutation and/or scaling to get a zero-free diagonal (0 to 7)","None",lu->id.ICNTL(6),&lu->id.ICNTL(6),PETSC_NULL);CHKERRQ(ierr); 512 icntl=-1; 513 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 514 if (flg) { 515 if (icntl== 1){ 516 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"pivot order be set by the user in PERM_IN -- not supported by the PETSc/MUMPS interface\n"); 517 } else { 518 lu->id.ICNTL(7) = icntl; 519 } 520 } 521 ierr = PetscOptionsInt("-mat_mumps_icntl_8","ICNTL(8): scaling strategy (-2 to 7 or 77)","None",lu->id.ICNTL(8),&lu->id.ICNTL(8),PETSC_NULL);CHKERRQ(ierr); 522 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); 523 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); 524 ierr = PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): statistics related to the linear system solved (via -ksp_view)","None",lu->id.ICNTL(11),&lu->id.ICNTL(11),PETSC_NULL);CHKERRQ(ierr); 525 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control: defines the ordering strategy with scaling constraints (0 to 3","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr); 526 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control: with or without ScaLAPACK","None",lu->id.ICNTL(13),&lu->id.ICNTL(13),PETSC_NULL);CHKERRQ(ierr); 527 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); 528 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",lu->id.ICNTL(19),&lu->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr); 529 530 ierr = PetscOptionsInt("-mat_mumps_icntl_22","ICNTL(22): in-core/out-of-core facility (0 or 1)","None",lu->id.ICNTL(22),&lu->id.ICNTL(22),PETSC_NULL);CHKERRQ(ierr); 531 ierr = PetscOptionsInt("-mat_mumps_icntl_23","ICNTL(23): max size of the working memory (MB) that can allocate per processor","None",lu->id.ICNTL(23),&lu->id.ICNTL(23),PETSC_NULL);CHKERRQ(ierr); 532 ierr = PetscOptionsInt("-mat_mumps_icntl_24","ICNTL(24): detection of null pivot rows (0 or 1)","None",lu->id.ICNTL(24),&lu->id.ICNTL(24),PETSC_NULL);CHKERRQ(ierr); 533 ierr = PetscOptionsInt("-mat_mumps_icntl_25","ICNTL(25): computation of a null space basis","None",lu->id.ICNTL(25),&lu->id.ICNTL(25),PETSC_NULL);CHKERRQ(ierr); 534 ierr = PetscOptionsInt("-mat_mumps_icntl_26","ICNTL(26): Schur options for right-hand side or solution vector","None",lu->id.ICNTL(26),&lu->id.ICNTL(26),PETSC_NULL);CHKERRQ(ierr); 535 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",lu->id.ICNTL(27),&lu->id.ICNTL(27),PETSC_NULL);CHKERRQ(ierr); 536 537 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 538 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); 539 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 540 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",lu->id.CNTL(4),&lu->id.CNTL(4),PETSC_NULL);CHKERRQ(ierr); 541 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",lu->id.CNTL(5),&lu->id.CNTL(5),PETSC_NULL);CHKERRQ(ierr); 542 PetscOptionsEnd(); 543 PetscFunctionReturn(0); 544 } 545 546 #undef __FUNCT__ 547 #define __FUNCT__ "PetscInitializeMUMPS" 548 PetscErrorCode PetscInitializeMUMPS(Mat F) 549 { 550 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 551 PetscErrorCode ierr; 552 PetscInt icntl; 553 PetscTruth flg; 554 555 PetscFunctionBegin; 556 lu->id.job = JOB_INIT; 557 lu->id.par=1; /* host participates factorizaton and solve */ 558 lu->id.sym=lu->sym; 559 if (lu->sym == 2){ 560 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 561 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 562 } 563 #if defined(PETSC_USE_COMPLEX) 564 zmumps_c(&lu->id); 565 #else 566 dmumps_c(&lu->id); 567 #endif 568 PetscFunctionReturn(0); 569 } 570 571 /* Note the Petsc r and c permutations are ignored */ 572 #undef __FUNCT__ 573 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 574 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 575 { 576 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 577 PetscErrorCode ierr; 578 PetscTruth isSeqAIJ,isMPIAIJ; 579 PetscInt nz=0,M=A->rmap->N,rnz,i,nnz; 580 const PetscInt *ai,*aj; 581 PetscTruth valOnly; 582 583 PetscFunctionBegin; 584 lu->sym = 0; 585 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 586 587 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 588 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 589 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 590 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 591 592 /* Initialize a MUMPS instance */ 593 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 594 /* Set MUMPS options */ 595 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 596 597 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 598 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 599 switch (lu->id.ICNTL(18)){ 600 case 0: /* centralized assembled matrix input (size=1) */ 601 if (!lu->myid) { 602 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 603 nz = aa->nz; 604 ai = aa->i; aj = aa->j; lu->val = aa->a; 605 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr); 606 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr); 607 nz = 0; 608 for (i=0; i<M; i++){ 609 rnz = ai[i+1] - ai[i]; 610 while (rnz--) { /* Fortran row/col index! */ 611 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 612 } 613 } 614 } 615 break; 616 case 3: /* distributed assembled matrix input (size>1) */ 617 valOnly = PETSC_FALSE; 618 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 619 break; 620 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 621 } 622 623 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 624 F->ops->solve = MatSolve_MUMPS; 625 PetscFunctionReturn(0); 626 } 627 628 /* Note the Petsc r and c permutations are ignored */ 629 #undef __FUNCT__ 630 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 631 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 632 { 633 634 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 635 PetscErrorCode ierr; 636 637 PetscFunctionBegin; 638 lu->sym = 0; 639 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 640 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 641 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 642 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 643 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 644 645 /* Initialize a MUMPS instance */ 646 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 647 /* Set MUMPS options */ 648 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 649 650 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 651 F->ops->solve = MatSolve_MUMPS; 652 PetscFunctionReturn(0); 653 } 654 655 /* Note the Petsc r permutation is ignored */ 656 #undef __FUNCT__ 657 #define __FUNCT__ "MatCholeskyFactorSymbolic_SBAIJMUMPS" 658 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJMUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 659 { 660 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 661 Mat newMat; 662 PetscErrorCode ierr; 663 PetscInt nz=0,M=A->rmap->N,rnz,i,nnz,jidx; 664 const PetscInt *ai,*aj,*adiag; 665 PetscScalar *av; 666 PetscTruth valOnly,isSeqAIJ,isMPIAIJ; 667 668 669 PetscFunctionBegin; 670 lu->sym = 2; 671 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 672 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 673 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 674 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 675 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 676 677 /* Initialize a MUMPS instance */ 678 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 679 /* Set MUMPS options */ 680 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 681 682 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 683 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 684 switch (lu->id.ICNTL(18)){ 685 case 0: /* centralized assembled matrix input (size=1) */ 686 if (!lu->myid) { 687 if(isSeqAIJ) { 688 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 689 nz = aa->nz; 690 ai = aa->i; aj = aa->j; adiag = aa->diag; av = aa->a; 691 } else { 692 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 693 nz = aa->nz; 694 ai = aa->i; aj = aa->j; adiag = aa->diag; av = aa->a; 695 } 696 if ((F)->factortype == MAT_FACTOR_CHOLESKY && isSeqAIJ) { 697 nz = M + (nz - M)/2; /* nz for upper triangle */ 698 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr); 699 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr); 700 ierr = PetscMalloc(nz*sizeof(PetscScalar),&lu->val);CHKERRQ(ierr); 701 nz = 0; 702 for (i=0; i<M; i++){ 703 rnz = ai[i+1] - adiag[i]; 704 jidx = adiag[i]; 705 while (rnz--) { /* Fortran row/col index! */ 706 lu->irn[nz] = i+1; lu->jcn[nz] = aj[jidx]+1; 707 lu->val[nz] = av[jidx]; jidx++; nz++; 708 } 709 } 710 } else { 711 lu->val = av; 712 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr); 713 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr); 714 nz = 0; 715 for (i=0; i<M; i++){ 716 rnz = ai[i+1] - ai[i]; 717 while (rnz--) { /* Fortran row/col index! */ 718 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 719 } 720 } 721 } 722 } 723 break; 724 case 3: /* distributed assembled matrix input (size>1) */ 725 valOnly = PETSC_FALSE; 726 if(((F)->factortype == MAT_FACTOR_CHOLESKY) && isMPIAIJ) { 727 /* Create an SBAIJ matrix and use this matrix to set the lu values */ 728 ierr = MatConvert(A,MATMPISBAIJ,MAT_INITIAL_MATRIX,&newMat);CHKERRQ(ierr); 729 ierr = MatConvertToTriples(newMat,1,valOnly,&nnz,&lu->irn , &lu->jcn, &lu->val);CHKERRQ(ierr); 730 ierr = MatDestroy(newMat);CHKERRQ(ierr); 731 } else { 732 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 733 } 734 break; 735 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 736 } 737 738 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 739 F->ops->solve = MatSolve_MUMPS; 740 #if !defined(PETSC_USE_COMPLEX) 741 (F)->ops->getinertia = MatGetInertia_SBAIJMUMPS; 742 #endif 743 PetscFunctionReturn(0); 744 } 745 746 #undef __FUNCT__ 747 #define __FUNCT__ "MatFactorInfo_MUMPS" 748 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 749 { 750 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 751 PetscErrorCode ierr; 752 753 PetscFunctionBegin; 754 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 755 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 756 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 757 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 758 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 759 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 760 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 761 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 762 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 763 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 764 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 765 if (!lu->myid && lu->id.ICNTL(11)>0) { 766 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 767 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 768 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 769 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); 770 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 771 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 772 773 } 774 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 775 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 776 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 777 /* ICNTL(15-17) not used */ 778 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 779 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 780 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 781 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 782 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 783 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 784 785 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 786 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 787 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 788 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 789 790 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 791 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 792 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 793 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 794 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 795 796 /* infomation local to each processor */ 797 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr);} 798 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 799 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 800 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr);} 801 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 802 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 803 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr);} 804 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 805 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 806 807 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr);} 808 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 809 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 810 811 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr);} 812 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 813 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 814 815 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr);} 816 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 817 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 818 819 if (!lu->myid){ /* information from the host */ 820 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 821 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 822 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 823 824 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 825 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 826 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 827 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 828 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 829 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 830 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 831 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 832 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 833 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 834 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 835 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 836 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 837 ierr = PetscViewerASCIIPrintf(viewer," INFOG(16) (estimated size (in MB) of all MUMPS internal data for factorization after analysis: value on the most memory consuming processor): %d \n",lu->id.INFOG(16));CHKERRQ(ierr); 838 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); 839 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); 840 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); 841 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 842 ierr = PetscViewerASCIIPrintf(viewer," INFOG(21) (size in MB of memory effectively used during factorization - value on the most memory consuming processor): %d \n",lu->id.INFOG(21));CHKERRQ(ierr); 843 ierr = PetscViewerASCIIPrintf(viewer," INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): %d \n",lu->id.INFOG(22));CHKERRQ(ierr); 844 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 845 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 846 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 847 } 848 PetscFunctionReturn(0); 849 } 850 851 #undef __FUNCT__ 852 #define __FUNCT__ "MatView_MUMPS" 853 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 854 { 855 PetscErrorCode ierr; 856 PetscTruth iascii; 857 PetscViewerFormat format; 858 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 859 860 PetscFunctionBegin; 861 /* check if matrix is mumps type */ 862 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 863 864 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 865 if (iascii) { 866 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 867 if (format == PETSC_VIEWER_ASCII_INFO || mumps->mumpsview){ 868 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 869 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 870 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 871 if (mumps->mumpsview){ /* View all MUMPS parameters */ 872 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 873 } 874 } 875 } 876 PetscFunctionReturn(0); 877 } 878 879 #undef __FUNCT__ 880 #define __FUNCT__ "MatGetInfo_MUMPS" 881 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 882 { 883 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 884 885 PetscFunctionBegin; 886 info->block_size = 1.0; 887 info->nz_allocated = mumps->id.INFOG(20); 888 info->nz_used = mumps->id.INFOG(20); 889 info->nz_unneeded = 0.0; 890 info->assemblies = 0.0; 891 info->mallocs = 0.0; 892 info->memory = 0.0; 893 info->fill_ratio_given = 0; 894 info->fill_ratio_needed = 0; 895 info->factor_mallocs = 0; 896 PetscFunctionReturn(0); 897 } 898 899 /*MC 900 MAT_SOLVER_MUMPS - A matrix type providing direct solvers (LU and Cholesky) for 901 distributed and sequential matrices via the external package MUMPS. 902 903 Works with MATAIJ and MATSBAIJ matrices 904 905 Options Database Keys: 906 + -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 907 . -mat_mumps_icntl_4 <0,...,4> - print level 908 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 909 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 910 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 911 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 912 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 913 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 914 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 915 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 916 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 917 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 918 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 919 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 920 921 Level: beginner 922 923 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 924 925 M*/ 926 927 EXTERN_C_BEGIN 928 #undef __FUNCT__ 929 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 930 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 931 { 932 PetscFunctionBegin; 933 *type = MAT_SOLVER_MUMPS; 934 PetscFunctionReturn(0); 935 } 936 EXTERN_C_END 937 938 EXTERN_C_BEGIN 939 /* 940 The seq and mpi versions of this function are the same 941 */ 942 #undef __FUNCT__ 943 #define __FUNCT__ "MatGetFactor_seqaij_mumps" 944 PetscErrorCode MatGetFactor_seqaij_mumps(Mat A,MatFactorType ftype,Mat *F) 945 { 946 Mat B; 947 PetscErrorCode ierr; 948 Mat_MUMPS *mumps; 949 950 PetscFunctionBegin; 951 /* Create the factorization matrix */ 952 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 953 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 954 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 955 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 956 957 B->ops->view = MatView_MUMPS; 958 B->ops->getinfo = MatGetInfo_MUMPS; 959 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 960 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 961 if (ftype == MAT_FACTOR_LU) { 962 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 963 B->factortype = MAT_FACTOR_LU; 964 } else { 965 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 966 B->factortype = MAT_FACTOR_CHOLESKY; 967 } 968 969 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 970 mumps->CleanUpMUMPS = PETSC_FALSE; 971 mumps->isAIJ = PETSC_TRUE; 972 mumps->scat_rhs = PETSC_NULL; 973 mumps->scat_sol = PETSC_NULL; 974 mumps->nSolve = 0; 975 mumps->MatDestroy = B->ops->destroy; 976 B->ops->destroy = MatDestroy_MUMPS; 977 B->spptr = (void*)mumps; 978 979 *F = B; 980 PetscFunctionReturn(0); 981 } 982 EXTERN_C_END 983 984 EXTERN_C_BEGIN 985 #undef __FUNCT__ 986 #define __FUNCT__ "MatGetFactor_seqsbaij_mumps" 987 PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 988 { 989 Mat B; 990 PetscErrorCode ierr; 991 Mat_MUMPS *mumps; 992 993 PetscFunctionBegin; 994 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 995 /* Create the factorization matrix */ 996 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 997 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 998 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 999 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1000 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1001 1002 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1003 B->ops->view = MatView_MUMPS; 1004 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1005 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1006 B->factortype = MAT_FACTOR_CHOLESKY; 1007 1008 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1009 mumps->CleanUpMUMPS = PETSC_FALSE; 1010 mumps->isAIJ = PETSC_FALSE; 1011 mumps->scat_rhs = PETSC_NULL; 1012 mumps->scat_sol = PETSC_NULL; 1013 mumps->nSolve = 0; 1014 mumps->MatDestroy = B->ops->destroy; 1015 B->ops->destroy = MatDestroy_MUMPS; 1016 B->spptr = (void*)mumps; 1017 1018 *F = B; 1019 PetscFunctionReturn(0); 1020 } 1021 EXTERN_C_END 1022 1023 EXTERN_C_BEGIN 1024 #undef __FUNCT__ 1025 #define __FUNCT__ "MatGetFactor_mpisbaij_mumps" 1026 PetscErrorCode MatGetFactor_mpisbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1027 { 1028 Mat B; 1029 PetscErrorCode ierr; 1030 Mat_MUMPS *mumps; 1031 1032 PetscFunctionBegin; 1033 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1034 /* Create the factorization matrix */ 1035 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1036 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1037 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1038 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1039 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1040 1041 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1042 B->ops->view = MatView_MUMPS; 1043 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1044 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1045 B->factortype = MAT_FACTOR_CHOLESKY; 1046 1047 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1048 mumps->CleanUpMUMPS = PETSC_FALSE; 1049 mumps->isAIJ = PETSC_FALSE; 1050 mumps->scat_rhs = PETSC_NULL; 1051 mumps->scat_sol = PETSC_NULL; 1052 mumps->nSolve = 0; 1053 mumps->MatDestroy = B->ops->destroy; 1054 B->ops->destroy = MatDestroy_MUMPS; 1055 B->spptr = (void*)mumps; 1056 1057 *F = B; 1058 PetscFunctionReturn(0); 1059 } 1060 EXTERN_C_END 1061 1062 EXTERN_C_BEGIN 1063 #undef __FUNCT__ 1064 #define __FUNCT__ "MatGetFactor_mpiaij_mumps" 1065 PetscErrorCode MatGetFactor_mpiaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1066 { 1067 Mat B; 1068 PetscErrorCode ierr; 1069 Mat_MUMPS *mumps; 1070 1071 PetscFunctionBegin; 1072 /* Create the factorization matrix */ 1073 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1074 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1075 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1076 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1077 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1078 1079 if (ftype == MAT_FACTOR_LU) { 1080 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1081 B->factortype = MAT_FACTOR_LU; 1082 } else { 1083 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1084 B->factortype = MAT_FACTOR_CHOLESKY; 1085 } 1086 1087 B->ops->view = MatView_MUMPS; 1088 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1089 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1090 1091 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1092 mumps->CleanUpMUMPS = PETSC_FALSE; 1093 mumps->isAIJ = PETSC_TRUE; 1094 mumps->scat_rhs = PETSC_NULL; 1095 mumps->scat_sol = PETSC_NULL; 1096 mumps->nSolve = 0; 1097 mumps->MatDestroy = B->ops->destroy; 1098 B->ops->destroy = MatDestroy_MUMPS; 1099 B->spptr = (void*)mumps; 1100 1101 *F = B; 1102 PetscFunctionReturn(0); 1103 } 1104 EXTERN_C_END 1105 1106 EXTERN_C_BEGIN 1107 #undef __FUNCT__ 1108 #define __FUNCT__ "MatGetFactor_mpibaij_mumps" 1109 PetscErrorCode MatGetFactor_mpibaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1110 { 1111 Mat B; 1112 PetscErrorCode ierr; 1113 Mat_MUMPS *mumps; 1114 1115 PetscFunctionBegin; 1116 /* Create the factorization matrix */ 1117 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1118 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1119 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1120 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1121 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1122 1123 if (ftype == MAT_FACTOR_LU) { 1124 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1125 B->factortype = MAT_FACTOR_LU; 1126 } 1127 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cholesky factorization for Petsc BAIJ matrices not supported yet\n"); 1128 B->ops->view = MatView_MUMPS; 1129 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1130 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1131 1132 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1133 mumps->CleanUpMUMPS = PETSC_FALSE; 1134 mumps->isAIJ = PETSC_TRUE; 1135 mumps->scat_rhs = PETSC_NULL; 1136 mumps->scat_sol = PETSC_NULL; 1137 mumps->nSolve = 0; 1138 mumps->MatDestroy = B->ops->destroy; 1139 B->ops->destroy = MatDestroy_MUMPS; 1140 B->spptr = (void*)mumps; 1141 1142 *F = B; 1143 PetscFunctionReturn(0); 1144 } 1145 EXTERN_C_END 1146 1147 /* -------------------------------------------------------------------------------------------*/ 1148 /*@ 1149 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1150 1151 Collective on Mat 1152 1153 Input Parameters: 1154 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1155 . idx - index of MUMPS parameter array ICNTL() 1156 - icntl - value of MUMPS ICNTL(imumps) 1157 1158 Options Database: 1159 . -mat_mumps_icntl_<idx> <icntl> 1160 1161 Level: beginner 1162 1163 References: MUMPS Users' Guide 1164 1165 .seealso: MatGetFactor() 1166 @*/ 1167 #undef __FUNCT__ 1168 #define __FUNCT__ "MatMumpsSetIcntl" 1169 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt idx,PetscInt icntl) 1170 { 1171 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 1172 1173 PetscFunctionBegin; 1174 lu->id.ICNTL(idx) = icntl; 1175 PetscFunctionReturn(0); 1176 } 1177 1178