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 { 341 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No mumps factorization for this matrix type"); 342 } 343 if ((F)->factortype == MAT_FACTOR_CHOLESKY && isSeqAIJ) { 344 nz = 0; 345 for (i=0; i<M; i++){ 346 rnz = ai[i+1] - adiag[i]; 347 jidx = adiag[i]; 348 while (rnz--) { /* Fortran row/col index! */ 349 lu->val[nz] = av[jidx]; jidx++; nz++; 350 } 351 } 352 } else { 353 lu->val = av; 354 } 355 } 356 break; 357 case 3: /* distributed assembled matrix input (size>1) */ 358 valOnly = PETSC_TRUE; /* only update mat values, not row and col index */ 359 360 if(((F)->factortype == MAT_FACTOR_CHOLESKY) && isMPIAIJ) { 361 /* Create an SBAIJ matrix and use this matrix to set the lu values */ 362 ierr = MatConvert(A,MATMPISBAIJ,MAT_INITIAL_MATRIX,&newMat);CHKERRQ(ierr); 363 ierr = MatConvertToTriples(newMat,1,valOnly,&nnz,&lu->irn , &lu->jcn, &lu->val);CHKERRQ(ierr); 364 ierr = MatDestroy(newMat);CHKERRQ(ierr); 365 } 366 else { 367 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 368 } 369 break; 370 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 371 } 372 373 /* analysis phase */ 374 /*----------------*/ 375 if (lu->matstruc == DIFFERENT_NONZERO_PATTERN){ 376 lu->id.job = 1; 377 378 lu->id.n = M; 379 switch (lu->id.ICNTL(18)){ 380 case 0: /* centralized assembled matrix input */ 381 if (!lu->myid) { 382 lu->id.nz =nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 383 if (lu->id.ICNTL(6)>1){ 384 #if defined(PETSC_USE_COMPLEX) 385 lu->id.a = (mumps_double_complex*)lu->val; 386 #else 387 lu->id.a = lu->val; 388 #endif 389 } 390 } 391 break; 392 case 3: /* distributed assembled matrix input (size>1) */ 393 lu->id.nz_loc = nnz; 394 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 395 if (lu->id.ICNTL(6)>1) { 396 #if defined(PETSC_USE_COMPLEX) 397 lu->id.a_loc = (mumps_double_complex*)lu->val; 398 #else 399 lu->id.a_loc = lu->val; 400 #endif 401 } 402 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 403 if (!lu->myid){ 404 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 405 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 406 } else { 407 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 408 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 409 } 410 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 411 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 412 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 413 414 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 415 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 416 ierr = VecDestroy(b);CHKERRQ(ierr); 417 break; 418 } 419 #if defined(PETSC_USE_COMPLEX) 420 zmumps_c(&lu->id); 421 #else 422 dmumps_c(&lu->id); 423 #endif 424 if (lu->id.INFOG(1) < 0) { 425 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 426 } 427 } 428 429 /* numerical factorization phase */ 430 /*-------------------------------*/ 431 lu->id.job = 2; 432 if(!lu->id.ICNTL(18)) { 433 if (!lu->myid) { 434 #if defined(PETSC_USE_COMPLEX) 435 lu->id.a = (mumps_double_complex*)lu->val; 436 #else 437 lu->id.a = lu->val; 438 #endif 439 } 440 } else { 441 #if defined(PETSC_USE_COMPLEX) 442 lu->id.a_loc = (mumps_double_complex*)lu->val; 443 #else 444 lu->id.a_loc = lu->val; 445 #endif 446 } 447 #if defined(PETSC_USE_COMPLEX) 448 zmumps_c(&lu->id); 449 #else 450 dmumps_c(&lu->id); 451 #endif 452 if (lu->id.INFOG(1) < 0) { 453 if (lu->id.INFO(1) == -13) { 454 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)); 455 } else { 456 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)); 457 } 458 } 459 460 if (!lu->myid && lu->id.ICNTL(16) > 0){ 461 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 462 } 463 464 if (lu->size > 1){ 465 if(isMPIAIJ) { 466 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 467 } else { 468 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 469 } 470 F_diag->assembled = PETSC_TRUE; 471 if (lu->nSolve){ 472 ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr); 473 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 474 ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr); 475 } 476 } 477 (F)->assembled = PETSC_TRUE; 478 lu->matstruc = SAME_NONZERO_PATTERN; 479 lu->CleanUpMUMPS = PETSC_TRUE; 480 lu->nSolve = 0; 481 PetscFunctionReturn(0); 482 } 483 484 #undef __FUNCT__ 485 #define __FUNCT__ "PetscSetMUMPSOptions" 486 PetscErrorCode PetscSetMUMPSOptions(Mat F, Mat A) 487 { 488 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 489 PetscErrorCode ierr; 490 PetscInt icntl; 491 PetscTruth flg; 492 493 PetscFunctionBegin; 494 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 495 ierr = PetscOptionsTruth("-mat_mumps_view","View MUMPS parameters","None",lu->mumpsview,&lu->mumpsview,PETSC_NULL);CHKERRQ(ierr); 496 if (lu->size == 1){ 497 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 498 } else { 499 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 500 } 501 502 icntl=-1; 503 lu->id.ICNTL(4) = 0; /* level of printing; overwrite mumps default ICNTL(4)=2 */ 504 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 505 if ((flg && icntl > 0) || PetscLogPrintInfo) { 506 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 507 } else { /* no output */ 508 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 509 lu->id.ICNTL(2) = 0; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 510 lu->id.ICNTL(3) = 0; /* output stream for global information, default=6 */ 511 } 512 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); 513 icntl=-1; 514 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 515 if (flg) { 516 if (icntl== 1){ 517 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"); 518 } else { 519 lu->id.ICNTL(7) = icntl; 520 } 521 } 522 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); 523 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); 524 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); 525 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); 526 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); 527 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); 528 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); 529 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",lu->id.ICNTL(19),&lu->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr); 530 531 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); 532 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); 533 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); 534 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); 535 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); 536 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",lu->id.ICNTL(27),&lu->id.ICNTL(27),PETSC_NULL);CHKERRQ(ierr); 537 538 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 539 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); 540 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 541 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); 542 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); 543 PetscOptionsEnd(); 544 PetscFunctionReturn(0); 545 } 546 547 #undef __FUNCT__ 548 #define __FUNCT__ "PetscInitializeMUMPS" 549 PetscErrorCode PetscInitializeMUMPS(Mat F) 550 { 551 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 552 PetscErrorCode ierr; 553 PetscInt icntl; 554 PetscTruth flg; 555 556 PetscFunctionBegin; 557 lu->id.job = JOB_INIT; 558 lu->id.par=1; /* host participates factorizaton and solve */ 559 lu->id.sym=lu->sym; 560 if (lu->sym == 2){ 561 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 562 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 563 } 564 #if defined(PETSC_USE_COMPLEX) 565 zmumps_c(&lu->id); 566 #else 567 dmumps_c(&lu->id); 568 #endif 569 PetscFunctionReturn(0); 570 } 571 572 /* Note the Petsc r and c permutations are ignored */ 573 #undef __FUNCT__ 574 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 575 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 576 { 577 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 578 PetscErrorCode ierr; 579 PetscTruth isSeqAIJ,isMPIAIJ; 580 PetscInt nz=0,M=A->rmap->N,rnz,i,nnz; 581 const PetscInt *ai,*aj; 582 PetscTruth valOnly; 583 584 PetscFunctionBegin; 585 lu->sym = 0; 586 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 587 588 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 589 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 590 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 591 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 592 593 /* Initialize a MUMPS instance */ 594 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 595 /* Set MUMPS options */ 596 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 597 598 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 599 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 600 switch (lu->id.ICNTL(18)){ 601 case 0: /* centralized assembled matrix input (size=1) */ 602 if (!lu->myid) { 603 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 604 nz = aa->nz; 605 ai = aa->i; aj = aa->j; lu->val = aa->a; 606 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr); 607 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr); 608 nz = 0; 609 for (i=0; i<M; i++){ 610 rnz = ai[i+1] - ai[i]; 611 while (rnz--) { /* Fortran row/col index! */ 612 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 613 } 614 } 615 } 616 break; 617 case 3: /* distributed assembled matrix input (size>1) */ 618 valOnly = PETSC_FALSE; 619 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 620 break; 621 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 622 } 623 624 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 625 F->ops->solve = MatSolve_MUMPS; 626 PetscFunctionReturn(0); 627 } 628 629 /* Note the Petsc r and c permutations are ignored */ 630 #undef __FUNCT__ 631 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 632 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 633 { 634 635 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 636 PetscErrorCode ierr; 637 638 PetscFunctionBegin; 639 lu->sym = 0; 640 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 641 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 642 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 643 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 644 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 645 646 /* Initialize a MUMPS instance */ 647 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 648 /* Set MUMPS options */ 649 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 650 651 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 652 F->ops->solve = MatSolve_MUMPS; 653 PetscFunctionReturn(0); 654 } 655 656 /* Note the Petsc r permutation is ignored */ 657 #undef __FUNCT__ 658 #define __FUNCT__ "MatCholeskyFactorSymbolic_SBAIJMUMPS" 659 PetscErrorCode MatCholeskyFactorSymbolic_SBAIJMUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 660 { 661 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 662 Mat newMat; 663 PetscErrorCode ierr; 664 PetscInt nz=0,M=A->rmap->N,rnz,i,nnz,jidx; 665 const PetscInt *ai,*aj,*adiag; 666 PetscScalar *av; 667 PetscTruth valOnly,isSeqAIJ,isMPIAIJ; 668 669 670 PetscFunctionBegin; 671 lu->sym = 2; 672 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 673 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 674 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 675 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 676 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 677 678 /* Initialize a MUMPS instance */ 679 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 680 /* Set MUMPS options */ 681 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 682 683 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 684 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 685 switch (lu->id.ICNTL(18)){ 686 case 0: /* centralized assembled matrix input (size=1) */ 687 if (!lu->myid) { 688 if(isSeqAIJ) { 689 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 690 nz = aa->nz; 691 ai = aa->i; aj = aa->j; adiag = aa->diag; av = aa->a; 692 } else { 693 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 694 nz = aa->nz; 695 ai = aa->i; aj = aa->j; adiag = aa->diag; av = aa->a; 696 } 697 if ((F)->factortype == MAT_FACTOR_CHOLESKY && isSeqAIJ) { 698 nz = M + (nz - M)/2; /* nz for upper triangle */ 699 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr); 700 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr); 701 ierr = PetscMalloc(nz*sizeof(PetscScalar),&lu->val);CHKERRQ(ierr); 702 nz = 0; 703 for (i=0; i<M; i++){ 704 rnz = ai[i+1] - adiag[i]; 705 jidx = adiag[i]; 706 while (rnz--) { /* Fortran row/col index! */ 707 lu->irn[nz] = i+1; lu->jcn[nz] = aj[jidx]+1; 708 lu->val[nz] = av[jidx]; jidx++; nz++; 709 } 710 } 711 } else { 712 lu->val = av; 713 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->irn);CHKERRQ(ierr); 714 ierr = PetscMalloc(nz*sizeof(PetscInt),&lu->jcn);CHKERRQ(ierr); 715 nz = 0; 716 for (i=0; i<M; i++){ 717 rnz = ai[i+1] - ai[i]; 718 while (rnz--) { /* Fortran row/col index! */ 719 lu->irn[nz] = i+1; lu->jcn[nz] = (*aj)+1; aj++; nz++; 720 } 721 } 722 } 723 } 724 break; 725 case 3: /* distributed assembled matrix input (size>1) */ 726 valOnly = PETSC_FALSE; 727 if(((F)->factortype == MAT_FACTOR_CHOLESKY) && isMPIAIJ) { 728 /* Create an SBAIJ matrix and use this matrix to set the lu values */ 729 ierr = MatConvert(A,MATMPISBAIJ,MAT_INITIAL_MATRIX,&newMat);CHKERRQ(ierr); 730 ierr = MatConvertToTriples(newMat,1,valOnly,&nnz,&lu->irn , &lu->jcn, &lu->val);CHKERRQ(ierr); 731 ierr = MatDestroy(newMat);CHKERRQ(ierr); 732 } else { 733 ierr = MatConvertToTriples(A,1,valOnly, &nnz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 734 } 735 break; 736 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix input format is not supported by MUMPS."); 737 } 738 739 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 740 F->ops->solve = MatSolve_MUMPS; 741 #if !defined(PETSC_USE_COMPLEX) 742 (F)->ops->getinertia = MatGetInertia_SBAIJMUMPS; 743 #endif 744 PetscFunctionReturn(0); 745 } 746 747 #undef __FUNCT__ 748 #define __FUNCT__ "MatFactorInfo_MUMPS" 749 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 750 { 751 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 752 PetscErrorCode ierr; 753 754 PetscFunctionBegin; 755 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 756 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 757 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 758 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 759 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 760 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 761 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 762 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 763 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 764 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 765 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 766 if (!lu->myid && lu->id.ICNTL(11)>0) { 767 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 768 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 769 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 770 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); 771 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 772 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 773 774 } 775 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 776 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 777 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 778 /* ICNTL(15-17) not used */ 779 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 780 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 781 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 782 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 783 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 784 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 785 786 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 787 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 788 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 789 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 790 791 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 792 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 793 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 794 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 795 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 796 797 /* infomation local to each processor */ 798 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr);} 799 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 800 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 801 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr);} 802 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 803 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 804 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr);} 805 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 806 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 807 808 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);} 809 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 810 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 811 812 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr);} 813 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 814 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 815 816 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr);} 817 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 818 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 819 820 if (!lu->myid){ /* information from the host */ 821 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 822 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 823 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 824 825 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 826 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 827 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 828 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 829 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 830 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 831 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 832 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 833 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 834 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 835 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 836 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 837 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 838 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); 839 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); 840 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); 841 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); 842 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 843 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); 844 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); 845 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 846 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 847 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 848 } 849 PetscFunctionReturn(0); 850 } 851 852 #undef __FUNCT__ 853 #define __FUNCT__ "MatView_MUMPS" 854 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 855 { 856 PetscErrorCode ierr; 857 PetscTruth iascii; 858 PetscViewerFormat format; 859 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 860 861 PetscFunctionBegin; 862 /* check if matrix is mumps type */ 863 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 864 865 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 866 if (iascii) { 867 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 868 if (format == PETSC_VIEWER_ASCII_INFO || mumps->mumpsview){ 869 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 870 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 871 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 872 if (mumps->mumpsview){ /* View all MUMPS parameters */ 873 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 874 } 875 } 876 } 877 PetscFunctionReturn(0); 878 } 879 880 #undef __FUNCT__ 881 #define __FUNCT__ "MatGetInfo_MUMPS" 882 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 883 { 884 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 885 886 PetscFunctionBegin; 887 info->block_size = 1.0; 888 info->nz_allocated = mumps->id.INFOG(20); 889 info->nz_used = mumps->id.INFOG(20); 890 info->nz_unneeded = 0.0; 891 info->assemblies = 0.0; 892 info->mallocs = 0.0; 893 info->memory = 0.0; 894 info->fill_ratio_given = 0; 895 info->fill_ratio_needed = 0; 896 info->factor_mallocs = 0; 897 PetscFunctionReturn(0); 898 } 899 900 /*MC 901 MAT_SOLVER_MUMPS - A matrix type providing direct solvers (LU and Cholesky) for 902 distributed and sequential matrices via the external package MUMPS. 903 904 Works with MATAIJ and MATSBAIJ matrices 905 906 Options Database Keys: 907 + -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 908 . -mat_mumps_icntl_4 <0,...,4> - print level 909 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 910 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 911 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 912 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 913 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 914 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 915 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 916 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 917 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 918 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 919 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 920 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 921 922 Level: beginner 923 924 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 925 926 M*/ 927 928 EXTERN_C_BEGIN 929 #undef __FUNCT__ 930 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 931 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 932 { 933 PetscFunctionBegin; 934 *type = MAT_SOLVER_MUMPS; 935 PetscFunctionReturn(0); 936 } 937 EXTERN_C_END 938 939 EXTERN_C_BEGIN 940 /* 941 The seq and mpi versions of this function are the same 942 */ 943 #undef __FUNCT__ 944 #define __FUNCT__ "MatGetFactor_seqaij_mumps" 945 PetscErrorCode MatGetFactor_seqaij_mumps(Mat A,MatFactorType ftype,Mat *F) 946 { 947 Mat B; 948 PetscErrorCode ierr; 949 Mat_MUMPS *mumps; 950 951 PetscFunctionBegin; 952 /* Create the factorization matrix */ 953 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 954 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 955 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 956 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 957 958 B->ops->view = MatView_MUMPS; 959 B->ops->getinfo = MatGetInfo_MUMPS; 960 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 961 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 962 if (ftype == MAT_FACTOR_LU) { 963 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 964 B->factortype = MAT_FACTOR_LU; 965 } else { 966 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 967 B->factortype = MAT_FACTOR_CHOLESKY; 968 } 969 970 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 971 mumps->CleanUpMUMPS = PETSC_FALSE; 972 mumps->isAIJ = PETSC_TRUE; 973 mumps->scat_rhs = PETSC_NULL; 974 mumps->scat_sol = PETSC_NULL; 975 mumps->nSolve = 0; 976 mumps->MatDestroy = B->ops->destroy; 977 B->ops->destroy = MatDestroy_MUMPS; 978 B->spptr = (void*)mumps; 979 980 *F = B; 981 PetscFunctionReturn(0); 982 } 983 EXTERN_C_END 984 985 EXTERN_C_BEGIN 986 #undef __FUNCT__ 987 #define __FUNCT__ "MatGetFactor_seqsbaij_mumps" 988 PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 989 { 990 Mat B; 991 PetscErrorCode ierr; 992 Mat_MUMPS *mumps; 993 994 PetscFunctionBegin; 995 if (ftype != MAT_FACTOR_CHOLESKY) { 996 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 997 } 998 /* Create the factorization matrix */ 999 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1000 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1001 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1002 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1003 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1004 1005 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1006 B->ops->view = MatView_MUMPS; 1007 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1008 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1009 B->factortype = MAT_FACTOR_CHOLESKY; 1010 1011 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1012 mumps->CleanUpMUMPS = PETSC_FALSE; 1013 mumps->isAIJ = PETSC_FALSE; 1014 mumps->scat_rhs = PETSC_NULL; 1015 mumps->scat_sol = PETSC_NULL; 1016 mumps->nSolve = 0; 1017 mumps->MatDestroy = B->ops->destroy; 1018 B->ops->destroy = MatDestroy_MUMPS; 1019 B->spptr = (void*)mumps; 1020 1021 *F = B; 1022 PetscFunctionReturn(0); 1023 } 1024 EXTERN_C_END 1025 1026 EXTERN_C_BEGIN 1027 #undef __FUNCT__ 1028 #define __FUNCT__ "MatGetFactor_mpisbaij_mumps" 1029 PetscErrorCode MatGetFactor_mpisbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1030 { 1031 Mat B; 1032 PetscErrorCode ierr; 1033 Mat_MUMPS *mumps; 1034 1035 PetscFunctionBegin; 1036 if (ftype != MAT_FACTOR_CHOLESKY) { 1037 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1038 } 1039 /* Create the factorization matrix */ 1040 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1041 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1042 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1043 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1044 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1045 1046 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1047 B->ops->view = MatView_MUMPS; 1048 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1049 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1050 B->factortype = MAT_FACTOR_CHOLESKY; 1051 1052 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1053 mumps->CleanUpMUMPS = PETSC_FALSE; 1054 mumps->isAIJ = PETSC_FALSE; 1055 mumps->scat_rhs = PETSC_NULL; 1056 mumps->scat_sol = PETSC_NULL; 1057 mumps->nSolve = 0; 1058 mumps->MatDestroy = B->ops->destroy; 1059 B->ops->destroy = MatDestroy_MUMPS; 1060 B->spptr = (void*)mumps; 1061 1062 *F = B; 1063 PetscFunctionReturn(0); 1064 } 1065 EXTERN_C_END 1066 1067 EXTERN_C_BEGIN 1068 #undef __FUNCT__ 1069 #define __FUNCT__ "MatGetFactor_mpiaij_mumps" 1070 PetscErrorCode MatGetFactor_mpiaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1071 { 1072 Mat B; 1073 PetscErrorCode ierr; 1074 Mat_MUMPS *mumps; 1075 1076 PetscFunctionBegin; 1077 /* Create the factorization matrix */ 1078 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1079 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1080 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1081 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1082 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1083 1084 if (ftype == MAT_FACTOR_LU) { 1085 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1086 B->factortype = MAT_FACTOR_LU; 1087 } else { 1088 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1089 B->factortype = MAT_FACTOR_CHOLESKY; 1090 } 1091 1092 B->ops->view = MatView_MUMPS; 1093 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1094 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1095 1096 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1097 mumps->CleanUpMUMPS = PETSC_FALSE; 1098 mumps->isAIJ = PETSC_TRUE; 1099 mumps->scat_rhs = PETSC_NULL; 1100 mumps->scat_sol = PETSC_NULL; 1101 mumps->nSolve = 0; 1102 mumps->MatDestroy = B->ops->destroy; 1103 B->ops->destroy = MatDestroy_MUMPS; 1104 B->spptr = (void*)mumps; 1105 1106 *F = B; 1107 PetscFunctionReturn(0); 1108 } 1109 EXTERN_C_END 1110 1111 EXTERN_C_BEGIN 1112 #undef __FUNCT__ 1113 #define __FUNCT__ "MatGetFactor_mpibaij_mumps" 1114 PetscErrorCode MatGetFactor_mpibaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1115 { 1116 Mat B; 1117 PetscErrorCode ierr; 1118 Mat_MUMPS *mumps; 1119 1120 PetscFunctionBegin; 1121 /* Create the factorization matrix */ 1122 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1123 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1124 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1125 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1126 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1127 1128 if (ftype == MAT_FACTOR_LU) { 1129 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1130 B->factortype = MAT_FACTOR_LU; 1131 } 1132 else { 1133 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cholesky factorization for Petsc BAIJ matrices not supported yet\n"); 1134 } 1135 B->ops->view = MatView_MUMPS; 1136 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1137 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1138 1139 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1140 mumps->CleanUpMUMPS = PETSC_FALSE; 1141 mumps->isAIJ = PETSC_TRUE; 1142 mumps->scat_rhs = PETSC_NULL; 1143 mumps->scat_sol = PETSC_NULL; 1144 mumps->nSolve = 0; 1145 mumps->MatDestroy = B->ops->destroy; 1146 B->ops->destroy = MatDestroy_MUMPS; 1147 B->spptr = (void*)mumps; 1148 1149 *F = B; 1150 PetscFunctionReturn(0); 1151 } 1152 EXTERN_C_END 1153 1154 /* -------------------------------------------------------------------------------------------*/ 1155 /*@ 1156 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1157 1158 Collective on Mat 1159 1160 Input Parameters: 1161 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1162 . idx - index of MUMPS parameter array ICNTL() 1163 - icntl - value of MUMPS ICNTL(imumps) 1164 1165 Options Database: 1166 . -mat_mumps_icntl_<idx> <icntl> 1167 1168 Level: beginner 1169 1170 References: MUMPS Users' Guide 1171 1172 .seealso: MatGetFactor() 1173 @*/ 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "MatMumpsSetIcntl" 1176 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt idx,PetscInt icntl) 1177 { 1178 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 1179 1180 PetscFunctionBegin; 1181 lu->id.ICNTL(idx) = icntl; 1182 PetscFunctionReturn(0); 1183 } 1184 1185