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; 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_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_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_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_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_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_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_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_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_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 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_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_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_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 /* check if matrix is mumps type */ 755 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 756 757 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 758 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 759 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 760 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 761 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 762 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 763 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 764 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 765 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 766 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 767 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 768 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 769 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 770 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 771 if (!lu->myid && lu->id.ICNTL(11)>0) { 772 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 773 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 774 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 775 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); 776 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 777 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 778 779 } 780 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 781 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 782 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 783 /* ICNTL(15-17) not used */ 784 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 785 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 786 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 787 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 788 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 789 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 790 791 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 792 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 793 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 794 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 795 796 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 797 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 798 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 799 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 800 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 801 802 /* infomation local to each processor */ 803 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr);} 804 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 805 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 806 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr);} 807 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 808 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 809 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr);} 810 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 811 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 812 813 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);} 814 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 815 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 816 817 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr);} 818 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 819 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 820 821 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr);} 822 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 823 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 824 825 if (!lu->myid){ /* information from the host */ 826 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 827 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 828 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 829 830 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 831 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 832 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 833 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 834 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 835 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 836 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 837 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 838 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 839 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 840 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 841 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 842 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 843 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); 844 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); 845 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); 846 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); 847 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 848 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); 849 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); 850 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 851 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 852 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 853 } 854 PetscFunctionReturn(0); 855 } 856 857 #undef __FUNCT__ 858 #define __FUNCT__ "MatView_MUMPS" 859 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 860 { 861 PetscErrorCode ierr; 862 PetscTruth iascii; 863 PetscViewerFormat format; 864 865 PetscFunctionBegin; 866 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 867 if (iascii) { 868 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 869 if (format == PETSC_VIEWER_ASCII_INFO){ 870 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 871 } 872 } 873 PetscFunctionReturn(0); 874 } 875 876 #undef __FUNCT__ 877 #define __FUNCT__ "MatGetInfo_MUMPS" 878 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 879 { 880 Mat_MUMPS *lu =(Mat_MUMPS*)A->spptr; 881 882 PetscFunctionBegin; 883 info->block_size = 1.0; 884 info->nz_allocated = lu->id.INFOG(20); 885 info->nz_used = lu->id.INFOG(20); 886 info->nz_unneeded = 0.0; 887 info->assemblies = 0.0; 888 info->mallocs = 0.0; 889 info->memory = 0.0; 890 info->fill_ratio_given = 0; 891 info->fill_ratio_needed = 0; 892 info->factor_mallocs = 0; 893 PetscFunctionReturn(0); 894 } 895 896 /*MC 897 MAT_SOLVER_MUMPS - A matrix type providing direct solvers (LU and Cholesky) for 898 distributed and sequential matrices via the external package MUMPS. 899 900 Works with MATAIJ and MATSBAIJ matrices 901 902 Options Database Keys: 903 + -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 904 . -mat_mumps_icntl_4 <0,...,4> - print level 905 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 906 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 907 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 908 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 909 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 910 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 911 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 912 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 913 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 914 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 915 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 916 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 917 918 Level: beginner 919 920 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 921 922 M*/ 923 924 EXTERN_C_BEGIN 925 #undef __FUNCT__ 926 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 927 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 928 { 929 PetscFunctionBegin; 930 *type = MAT_SOLVER_MUMPS; 931 PetscFunctionReturn(0); 932 } 933 EXTERN_C_END 934 935 EXTERN_C_BEGIN 936 /* 937 The seq and mpi versions of this function are the same 938 */ 939 #undef __FUNCT__ 940 #define __FUNCT__ "MatGetFactor_seqaij_mumps" 941 PetscErrorCode MatGetFactor_seqaij_mumps(Mat A,MatFactorType ftype,Mat *F) 942 { 943 Mat B; 944 PetscErrorCode ierr; 945 Mat_MUMPS *mumps; 946 947 PetscFunctionBegin; 948 /* Create the factorization matrix */ 949 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 950 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 951 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 952 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 953 954 B->ops->view = MatView_MUMPS; 955 B->ops->getinfo = MatGetInfo_MUMPS; 956 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 957 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 958 if (ftype == MAT_FACTOR_LU) { 959 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 960 B->factortype = MAT_FACTOR_LU; 961 } else { 962 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 963 B->factortype = MAT_FACTOR_CHOLESKY; 964 } 965 966 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 967 mumps->CleanUpMUMPS = PETSC_FALSE; 968 mumps->isAIJ = PETSC_TRUE; 969 mumps->scat_rhs = PETSC_NULL; 970 mumps->scat_sol = PETSC_NULL; 971 mumps->nSolve = 0; 972 mumps->MatDestroy = B->ops->destroy; 973 B->ops->destroy = MatDestroy_MUMPS; 974 B->spptr = (void*)mumps; 975 976 *F = B; 977 PetscFunctionReturn(0); 978 } 979 EXTERN_C_END 980 981 EXTERN_C_BEGIN 982 #undef __FUNCT__ 983 #define __FUNCT__ "MatGetFactor_seqsbaij_mumps" 984 PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 985 { 986 Mat B; 987 PetscErrorCode ierr; 988 Mat_MUMPS *mumps; 989 990 PetscFunctionBegin; 991 if (ftype != MAT_FACTOR_CHOLESKY) { 992 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 993 } 994 /* Create the factorization matrix */ 995 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 996 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 997 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 998 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 999 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1000 1001 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1002 B->ops->view = MatView_MUMPS; 1003 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1004 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1005 B->factortype = MAT_FACTOR_CHOLESKY; 1006 1007 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1008 mumps->CleanUpMUMPS = PETSC_FALSE; 1009 mumps->isAIJ = PETSC_FALSE; 1010 mumps->scat_rhs = PETSC_NULL; 1011 mumps->scat_sol = PETSC_NULL; 1012 mumps->nSolve = 0; 1013 mumps->MatDestroy = B->ops->destroy; 1014 B->ops->destroy = MatDestroy_MUMPS; 1015 B->spptr = (void*)mumps; 1016 1017 *F = B; 1018 PetscFunctionReturn(0); 1019 } 1020 EXTERN_C_END 1021 1022 EXTERN_C_BEGIN 1023 #undef __FUNCT__ 1024 #define __FUNCT__ "MatGetFactor_mpisbaij_mumps" 1025 PetscErrorCode MatGetFactor_mpisbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1026 { 1027 Mat B; 1028 PetscErrorCode ierr; 1029 Mat_MUMPS *mumps; 1030 1031 PetscFunctionBegin; 1032 if (ftype != MAT_FACTOR_CHOLESKY) { 1033 SETERRQ(PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1034 } 1035 /* Create the factorization matrix */ 1036 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1037 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1038 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1039 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1040 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1041 1042 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1043 B->ops->view = MatView_MUMPS; 1044 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1045 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1046 B->factortype = MAT_FACTOR_CHOLESKY; 1047 1048 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1049 mumps->CleanUpMUMPS = PETSC_FALSE; 1050 mumps->isAIJ = PETSC_FALSE; 1051 mumps->scat_rhs = PETSC_NULL; 1052 mumps->scat_sol = PETSC_NULL; 1053 mumps->nSolve = 0; 1054 mumps->MatDestroy = B->ops->destroy; 1055 B->ops->destroy = MatDestroy_MUMPS; 1056 B->spptr = (void*)mumps; 1057 1058 *F = B; 1059 PetscFunctionReturn(0); 1060 } 1061 EXTERN_C_END 1062 1063 EXTERN_C_BEGIN 1064 #undef __FUNCT__ 1065 #define __FUNCT__ "MatGetFactor_mpiaij_mumps" 1066 PetscErrorCode MatGetFactor_mpiaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1067 { 1068 Mat B; 1069 PetscErrorCode ierr; 1070 Mat_MUMPS *mumps; 1071 1072 PetscFunctionBegin; 1073 /* Create the factorization matrix */ 1074 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1075 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1076 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1077 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1078 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1079 1080 if (ftype == MAT_FACTOR_LU) { 1081 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1082 B->factortype = MAT_FACTOR_LU; 1083 } else { 1084 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SBAIJMUMPS; 1085 B->factortype = MAT_FACTOR_CHOLESKY; 1086 } 1087 1088 B->ops->view = MatView_MUMPS; 1089 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1090 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1091 1092 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1093 mumps->CleanUpMUMPS = PETSC_FALSE; 1094 mumps->isAIJ = PETSC_TRUE; 1095 mumps->scat_rhs = PETSC_NULL; 1096 mumps->scat_sol = PETSC_NULL; 1097 mumps->nSolve = 0; 1098 mumps->MatDestroy = B->ops->destroy; 1099 B->ops->destroy = MatDestroy_MUMPS; 1100 B->spptr = (void*)mumps; 1101 1102 *F = B; 1103 PetscFunctionReturn(0); 1104 } 1105 EXTERN_C_END 1106 1107 EXTERN_C_BEGIN 1108 #undef __FUNCT__ 1109 #define __FUNCT__ "MatGetFactor_mpibaij_mumps" 1110 PetscErrorCode MatGetFactor_mpibaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1111 { 1112 Mat B; 1113 PetscErrorCode ierr; 1114 Mat_MUMPS *mumps; 1115 1116 PetscFunctionBegin; 1117 /* Create the factorization matrix */ 1118 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1119 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1120 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1121 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1122 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1123 1124 if (ftype == MAT_FACTOR_LU) { 1125 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1126 B->factortype = MAT_FACTOR_LU; 1127 } 1128 else { 1129 SETERRQ(PETSC_ERR_SUP,"Cholesky factorization for Petsc BAIJ matrices not supported yet\n"); 1130 } 1131 B->ops->view = MatView_MUMPS; 1132 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1133 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1134 1135 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1136 mumps->CleanUpMUMPS = PETSC_FALSE; 1137 mumps->isAIJ = PETSC_TRUE; 1138 mumps->scat_rhs = PETSC_NULL; 1139 mumps->scat_sol = PETSC_NULL; 1140 mumps->nSolve = 0; 1141 mumps->MatDestroy = B->ops->destroy; 1142 B->ops->destroy = MatDestroy_MUMPS; 1143 B->spptr = (void*)mumps; 1144 1145 *F = B; 1146 PetscFunctionReturn(0); 1147 } 1148 EXTERN_C_END 1149 1150 /* -------------------------------------------------------------------------------------------*/ 1151 /*@ 1152 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1153 1154 Collective on Mat 1155 1156 Input Parameters: 1157 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1158 . idx - index of MUMPS parameter array ICNTL() 1159 - icntl - value of MUMPS ICNTL(imumps) 1160 1161 Options Database: 1162 . -mat_mumps_icntl_<idx> <icntl> 1163 1164 Level: beginner 1165 1166 References: MUMPS Users' Guide 1167 1168 .seealso: MatGetFactor() 1169 @*/ 1170 #undef __FUNCT__ 1171 #define __FUNCT__ "MatMumpsSetIcntl" 1172 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt idx,PetscInt icntl) 1173 { 1174 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 1175 1176 PetscFunctionBegin; 1177 lu->id.ICNTL(idx) = icntl; 1178 PetscFunctionReturn(0); 1179 } 1180 1181