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