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