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 PetscTruth isAIJ,CleanUpMUMPS,mumpsview; 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 PetscMPIInt size=lu->size; 487 488 PetscFunctionBegin; 489 if (lu->CleanUpMUMPS) { 490 /* Terminate instance, deallocate memories */ 491 if (size > 1){ 492 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 493 ierr = VecScatterDestroy(lu->scat_rhs);CHKERRQ(ierr); 494 ierr = VecDestroy(lu->b_seq);CHKERRQ(ierr); 495 if (lu->nSolve && lu->scat_sol){ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);} 496 if (lu->nSolve && lu->x_seq){ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);} 497 } 498 ierr = PetscFree(lu->irn);CHKERRQ(ierr); 499 lu->id.job=JOB_END; 500 #if defined(PETSC_USE_COMPLEX) 501 zmumps_c(&lu->id); 502 #else 503 dmumps_c(&lu->id); 504 #endif 505 ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr); 506 } 507 /* clear composed functions */ 508 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatFactorGetSolverPackage_C","",PETSC_NULL);CHKERRQ(ierr); 509 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSetMumpsIcntl_C","",PETSC_NULL);CHKERRQ(ierr); 510 ierr = (lu->MatDestroy)(A);CHKERRQ(ierr); 511 PetscFunctionReturn(0); 512 } 513 514 #undef __FUNCT__ 515 #define __FUNCT__ "MatSolve_MUMPS" 516 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 517 { 518 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 519 PetscScalar *array; 520 Vec b_seq; 521 IS is_iden,is_petsc; 522 PetscErrorCode ierr; 523 PetscInt i; 524 525 PetscFunctionBegin; 526 lu->id.nrhs = 1; 527 b_seq = lu->b_seq; 528 if (lu->size > 1){ 529 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 530 ierr = VecScatterBegin(lu->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 531 ierr = VecScatterEnd(lu->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 532 if (!lu->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 533 } else { /* size == 1 */ 534 ierr = VecCopy(b,x);CHKERRQ(ierr); 535 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 536 } 537 if (!lu->myid) { /* define rhs on the host */ 538 lu->id.nrhs = 1; 539 #if defined(PETSC_USE_COMPLEX) 540 lu->id.rhs = (mumps_double_complex*)array; 541 #else 542 lu->id.rhs = array; 543 #endif 544 } 545 546 /* solve phase */ 547 /*-------------*/ 548 lu->id.job = JOB_SOLVE; 549 #if defined(PETSC_USE_COMPLEX) 550 zmumps_c(&lu->id); 551 #else 552 dmumps_c(&lu->id); 553 #endif 554 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)); 555 556 if (lu->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 557 if (!lu->nSolve){ /* create scatter scat_sol */ 558 ierr = ISCreateStride(PETSC_COMM_SELF,lu->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 559 for (i=0; i<lu->id.lsol_loc; i++){ 560 lu->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 561 } 562 ierr = ISCreateGeneral(PETSC_COMM_SELF,lu->id.lsol_loc,lu->id.isol_loc,&is_petsc);CHKERRQ(ierr); /* to */ 563 ierr = VecScatterCreate(lu->x_seq,is_iden,x,is_petsc,&lu->scat_sol);CHKERRQ(ierr); 564 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 565 ierr = ISDestroy(is_petsc);CHKERRQ(ierr); 566 } 567 ierr = VecScatterBegin(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 568 ierr = VecScatterEnd(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 569 } 570 lu->nSolve++; 571 PetscFunctionReturn(0); 572 } 573 574 #if !defined(PETSC_USE_COMPLEX) 575 /* 576 input: 577 F: numeric factor 578 output: 579 nneg: total number of negative pivots 580 nzero: 0 581 npos: (global dimension of F) - nneg 582 */ 583 584 #undef __FUNCT__ 585 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 586 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 587 { 588 Mat_MUMPS *lu =(Mat_MUMPS*)F->spptr; 589 PetscErrorCode ierr; 590 PetscMPIInt size; 591 592 PetscFunctionBegin; 593 ierr = MPI_Comm_size(((PetscObject)F)->comm,&size);CHKERRQ(ierr); 594 /* 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 */ 595 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)); 596 if (nneg){ 597 if (!lu->myid){ 598 *nneg = lu->id.INFOG(12); 599 } 600 ierr = MPI_Bcast(nneg,1,MPI_INT,0,lu->comm_mumps);CHKERRQ(ierr); 601 } 602 if (nzero) *nzero = 0; 603 if (npos) *npos = F->rmap->N - (*nneg); 604 PetscFunctionReturn(0); 605 } 606 #endif /* !defined(PETSC_USE_COMPLEX) */ 607 608 #undef __FUNCT__ 609 #define __FUNCT__ "MatFactorNumeric_MUMPS" 610 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 611 { 612 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 613 PetscErrorCode ierr; 614 MatReuse reuse; 615 Mat F_diag; 616 PetscTruth isMPIAIJ; 617 618 PetscFunctionBegin; 619 reuse = MAT_REUSE_MATRIX; 620 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 621 622 /* numerical factorization phase */ 623 /*-------------------------------*/ 624 lu->id.job = JOB_FACTNUMERIC; 625 if(!lu->id.ICNTL(18)) { 626 if (!lu->myid) { 627 #if defined(PETSC_USE_COMPLEX) 628 lu->id.a = (mumps_double_complex*)lu->val; 629 #else 630 lu->id.a = lu->val; 631 #endif 632 } 633 } else { 634 #if defined(PETSC_USE_COMPLEX) 635 lu->id.a_loc = (mumps_double_complex*)lu->val; 636 #else 637 lu->id.a_loc = lu->val; 638 #endif 639 } 640 #if defined(PETSC_USE_COMPLEX) 641 zmumps_c(&lu->id); 642 #else 643 dmumps_c(&lu->id); 644 #endif 645 if (lu->id.INFOG(1) < 0) { 646 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)); 647 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)); 648 } 649 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)); 650 651 if (lu->size > 1){ 652 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 653 if(isMPIAIJ) { 654 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 655 } else { 656 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 657 } 658 F_diag->assembled = PETSC_TRUE; 659 if (lu->nSolve){ 660 ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr); 661 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 662 ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr); 663 } 664 } 665 (F)->assembled = PETSC_TRUE; 666 lu->matstruc = SAME_NONZERO_PATTERN; 667 lu->CleanUpMUMPS = PETSC_TRUE; 668 lu->nSolve = 0; 669 670 if (lu->size > 1){ 671 /* distributed solution */ 672 lu->id.ICNTL(21) = 1; 673 if (!lu->nSolve){ 674 /* Create x_seq=sol_loc for repeated use */ 675 PetscInt lsol_loc; 676 PetscScalar *sol_loc; 677 lsol_loc = lu->id.INFO(23); /* length of sol_loc */ 678 ierr = PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&lu->id.isol_loc);CHKERRQ(ierr); 679 lu->id.lsol_loc = lsol_loc; 680 #if defined(PETSC_USE_COMPLEX) 681 lu->id.sol_loc = (mumps_double_complex*)sol_loc; 682 #else 683 lu->id.sol_loc = sol_loc; 684 #endif 685 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,lsol_loc,sol_loc,&lu->x_seq);CHKERRQ(ierr); 686 } 687 } 688 PetscFunctionReturn(0); 689 } 690 691 #undef __FUNCT__ 692 #define __FUNCT__ "PetscSetMUMPSOptions" 693 PetscErrorCode PetscSetMUMPSOptions(Mat F, Mat A) 694 { 695 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 696 PetscErrorCode ierr; 697 PetscInt icntl; 698 PetscTruth flg; 699 700 PetscFunctionBegin; 701 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 702 ierr = PetscOptionsTruth("-mat_mumps_view","View MUMPS parameters","None",lu->mumpsview,&lu->mumpsview,PETSC_NULL);CHKERRQ(ierr); 703 if (lu->size == 1){ 704 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 705 } else { 706 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 707 } 708 709 icntl=-1; 710 lu->id.ICNTL(4) = 0; /* level of printing; overwrite mumps default ICNTL(4)=2 */ 711 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 712 if ((flg && icntl > 0) || PetscLogPrintInfo) { 713 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 714 } else { /* no output */ 715 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 716 lu->id.ICNTL(2) = 0; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 717 lu->id.ICNTL(3) = 0; /* output stream for global information, default=6 */ 718 } 719 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); 720 icntl=-1; 721 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 722 if (flg) { 723 if (icntl== 1){ 724 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"); 725 } else { 726 lu->id.ICNTL(7) = icntl; 727 } 728 } 729 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); 730 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); 731 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); 732 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); 733 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); 734 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); 735 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); 736 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",lu->id.ICNTL(19),&lu->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr); 737 738 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); 739 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); 740 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); 741 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); 742 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); 743 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",lu->id.ICNTL(27),&lu->id.ICNTL(27),PETSC_NULL);CHKERRQ(ierr); 744 745 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 746 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); 747 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 748 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); 749 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); 750 PetscOptionsEnd(); 751 PetscFunctionReturn(0); 752 } 753 754 #undef __FUNCT__ 755 #define __FUNCT__ "PetscInitializeMUMPS" 756 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS* mumps) 757 { 758 PetscErrorCode ierr; 759 760 PetscFunctionBegin; 761 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &mumps->myid); 762 ierr = MPI_Comm_size(((PetscObject)A)->comm,&mumps->size);CHKERRQ(ierr); 763 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(mumps->comm_mumps));CHKERRQ(ierr); 764 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 765 766 mumps->id.job = JOB_INIT; 767 mumps->id.par = 1; /* host participates factorizaton and solve */ 768 mumps->id.sym = mumps->sym; 769 #if defined(PETSC_USE_COMPLEX) 770 zmumps_c(&mumps->id); 771 #else 772 dmumps_c(&mumps->id); 773 #endif 774 775 mumps->CleanUpMUMPS = PETSC_FALSE; 776 mumps->scat_rhs = PETSC_NULL; 777 mumps->scat_sol = PETSC_NULL; 778 mumps->nSolve = 0; 779 PetscFunctionReturn(0); 780 } 781 782 /* Note the Petsc r and c permutations are ignored */ 783 #undef __FUNCT__ 784 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 785 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 786 { 787 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 788 PetscErrorCode ierr; 789 MatReuse reuse; 790 Vec b; 791 IS is_iden; 792 const PetscInt M = A->rmap->N; 793 794 PetscFunctionBegin; 795 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 796 797 /* Set MUMPS options */ 798 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 799 800 reuse = MAT_INITIAL_MATRIX; 801 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 802 803 /* analysis phase */ 804 /*----------------*/ 805 lu->id.job = JOB_FACTSYMBOLIC; 806 lu->id.n = M; 807 switch (lu->id.ICNTL(18)){ 808 case 0: /* centralized assembled matrix input */ 809 if (!lu->myid) { 810 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 811 if (lu->id.ICNTL(6)>1){ 812 #if defined(PETSC_USE_COMPLEX) 813 lu->id.a = (mumps_double_complex*)lu->val; 814 #else 815 lu->id.a = lu->val; 816 #endif 817 } 818 } 819 break; 820 case 3: /* distributed assembled matrix input (size>1) */ 821 lu->id.nz_loc = lu->nz; 822 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 823 if (lu->id.ICNTL(6)>1) { 824 #if defined(PETSC_USE_COMPLEX) 825 lu->id.a_loc = (mumps_double_complex*)lu->val; 826 #else 827 lu->id.a_loc = lu->val; 828 #endif 829 } 830 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 831 if (!lu->myid){ 832 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 833 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 834 } else { 835 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 836 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 837 } 838 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 839 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 840 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 841 842 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 843 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 844 ierr = VecDestroy(b);CHKERRQ(ierr); 845 break; 846 } 847 #if defined(PETSC_USE_COMPLEX) 848 zmumps_c(&lu->id); 849 #else 850 dmumps_c(&lu->id); 851 #endif 852 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)); 853 854 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 855 F->ops->solve = MatSolve_MUMPS; 856 PetscFunctionReturn(0); 857 } 858 859 /* Note the Petsc r and c permutations are ignored */ 860 #undef __FUNCT__ 861 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 862 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 863 { 864 865 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 866 PetscErrorCode ierr; 867 MatReuse reuse; 868 Vec b; 869 IS is_iden; 870 const PetscInt M = A->rmap->N; 871 872 PetscFunctionBegin; 873 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 874 875 /* Set MUMPS options */ 876 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 877 878 reuse = MAT_INITIAL_MATRIX; 879 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 880 881 /* analysis phase */ 882 /*----------------*/ 883 lu->id.job = JOB_FACTSYMBOLIC; 884 lu->id.n = M; 885 switch (lu->id.ICNTL(18)){ 886 case 0: /* centralized assembled matrix input */ 887 if (!lu->myid) { 888 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 889 if (lu->id.ICNTL(6)>1){ 890 #if defined(PETSC_USE_COMPLEX) 891 lu->id.a = (mumps_double_complex*)lu->val; 892 #else 893 lu->id.a = lu->val; 894 #endif 895 } 896 } 897 break; 898 case 3: /* distributed assembled matrix input (size>1) */ 899 lu->id.nz_loc = lu->nz; 900 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 901 if (lu->id.ICNTL(6)>1) { 902 #if defined(PETSC_USE_COMPLEX) 903 lu->id.a_loc = (mumps_double_complex*)lu->val; 904 #else 905 lu->id.a_loc = lu->val; 906 #endif 907 } 908 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 909 if (!lu->myid){ 910 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 911 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 912 } else { 913 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 914 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 915 } 916 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 917 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 918 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 919 920 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 921 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 922 ierr = VecDestroy(b);CHKERRQ(ierr); 923 break; 924 } 925 #if defined(PETSC_USE_COMPLEX) 926 zmumps_c(&lu->id); 927 #else 928 dmumps_c(&lu->id); 929 #endif 930 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)); 931 932 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 933 F->ops->solve = MatSolve_MUMPS; 934 PetscFunctionReturn(0); 935 } 936 937 /* Note the Petsc r permutation is ignored */ 938 #undef __FUNCT__ 939 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 940 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 941 { 942 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 943 PetscErrorCode ierr; 944 MatReuse reuse; 945 Vec b; 946 IS is_iden; 947 const PetscInt M = A->rmap->N; 948 949 PetscFunctionBegin; 950 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 951 952 /* Set MUMPS options */ 953 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 954 955 reuse = MAT_INITIAL_MATRIX; 956 ierr = (*lu->ConvertToTriples)(A, 1 , reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 957 958 /* analysis phase */ 959 /*----------------*/ 960 lu->id.job = JOB_FACTSYMBOLIC; 961 lu->id.n = M; 962 switch (lu->id.ICNTL(18)){ 963 case 0: /* centralized assembled matrix input */ 964 if (!lu->myid) { 965 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 966 if (lu->id.ICNTL(6)>1){ 967 #if defined(PETSC_USE_COMPLEX) 968 lu->id.a = (mumps_double_complex*)lu->val; 969 #else 970 lu->id.a = lu->val; 971 #endif 972 } 973 } 974 break; 975 case 3: /* distributed assembled matrix input (size>1) */ 976 lu->id.nz_loc = lu->nz; 977 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 978 if (lu->id.ICNTL(6)>1) { 979 #if defined(PETSC_USE_COMPLEX) 980 lu->id.a_loc = (mumps_double_complex*)lu->val; 981 #else 982 lu->id.a_loc = lu->val; 983 #endif 984 } 985 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 986 if (!lu->myid){ 987 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 988 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 989 } else { 990 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 991 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 992 } 993 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 994 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 995 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 996 997 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 998 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 999 ierr = VecDestroy(b);CHKERRQ(ierr); 1000 break; 1001 } 1002 #if defined(PETSC_USE_COMPLEX) 1003 zmumps_c(&lu->id); 1004 #else 1005 dmumps_c(&lu->id); 1006 #endif 1007 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)); 1008 1009 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1010 F->ops->solve = MatSolve_MUMPS; 1011 #if !defined(PETSC_USE_COMPLEX) 1012 (F)->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1013 #endif 1014 PetscFunctionReturn(0); 1015 } 1016 1017 #undef __FUNCT__ 1018 #define __FUNCT__ "MatFactorInfo_MUMPS" 1019 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 1020 { 1021 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 1022 PetscErrorCode ierr; 1023 1024 PetscFunctionBegin; 1025 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 1026 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 1027 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 1028 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 1029 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 1030 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 1031 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 1032 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 1033 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 1034 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 1035 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 1036 if (lu->id.ICNTL(11)>0) { 1037 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 1038 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 1039 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 1040 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",lu->id.RINFOG(7),lu->id.RINFOG(8));CHKERRQ(ierr); 1041 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 1042 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 1043 1044 } 1045 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 1046 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 1047 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 1048 /* ICNTL(15-17) not used */ 1049 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 1050 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 1051 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 1052 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 1053 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 1054 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 1055 1056 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 1057 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 1058 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 1059 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 1060 1061 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 1062 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 1063 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 1064 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 1065 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 1066 1067 /* infomation local to each processor */ 1068 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1069 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 1070 ierr = PetscViewerFlush(viewer); 1071 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1072 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 1073 ierr = PetscViewerFlush(viewer); 1074 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1075 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 1076 ierr = PetscViewerFlush(viewer); 1077 1078 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1079 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 1080 ierr = PetscViewerFlush(viewer); 1081 1082 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1083 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 1084 ierr = PetscViewerFlush(viewer); 1085 1086 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1087 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 1088 ierr = PetscViewerFlush(viewer); 1089 1090 if (!lu->myid){ /* information from the host */ 1091 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 1092 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 1093 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 1094 1095 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 1096 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 1097 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 1098 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 1099 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 1100 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 1101 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 1102 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 1103 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 1104 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 1105 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 1106 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 1107 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 1108 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); 1109 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); 1110 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); 1111 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); 1112 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 1113 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); 1114 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); 1115 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 1116 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 1117 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 1118 } 1119 PetscFunctionReturn(0); 1120 } 1121 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "MatView_MUMPS" 1124 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1125 { 1126 PetscErrorCode ierr; 1127 PetscTruth iascii; 1128 PetscViewerFormat format; 1129 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1130 1131 PetscFunctionBegin; 1132 /* check if matrix is mumps type */ 1133 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1134 1135 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1136 if (iascii) { 1137 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1138 if (format == PETSC_VIEWER_ASCII_INFO || mumps->mumpsview){ 1139 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1140 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1141 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1142 if (mumps->mumpsview){ /* View all MUMPS parameters */ 1143 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 1144 } 1145 } 1146 } 1147 PetscFunctionReturn(0); 1148 } 1149 1150 #undef __FUNCT__ 1151 #define __FUNCT__ "MatGetInfo_MUMPS" 1152 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1153 { 1154 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1155 1156 PetscFunctionBegin; 1157 info->block_size = 1.0; 1158 info->nz_allocated = mumps->id.INFOG(20); 1159 info->nz_used = mumps->id.INFOG(20); 1160 info->nz_unneeded = 0.0; 1161 info->assemblies = 0.0; 1162 info->mallocs = 0.0; 1163 info->memory = 0.0; 1164 info->fill_ratio_given = 0; 1165 info->fill_ratio_needed = 0; 1166 info->factor_mallocs = 0; 1167 PetscFunctionReturn(0); 1168 } 1169 1170 /*MC 1171 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1172 distributed and sequential matrices via the external package MUMPS. 1173 1174 Works with MATAIJ and MATSBAIJ matrices 1175 1176 Options Database Keys: 1177 + -mat_mumps_icntl_4 <0,...,4> - print level 1178 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1179 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 1180 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1181 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1182 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1183 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1184 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1185 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1186 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1187 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1188 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1189 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1190 1191 Level: beginner 1192 1193 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1194 1195 M*/ 1196 1197 EXTERN_C_BEGIN 1198 #undef __FUNCT__ 1199 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1200 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1201 { 1202 PetscFunctionBegin; 1203 *type = MATSOLVERMUMPS; 1204 PetscFunctionReturn(0); 1205 } 1206 EXTERN_C_END 1207 1208 EXTERN_C_BEGIN 1209 /* MatGetFactor for Seq and MPI AIJ matrices */ 1210 #undef __FUNCT__ 1211 #define __FUNCT__ "MatGetFactor_aij_mumps" 1212 PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1213 { 1214 Mat B; 1215 PetscErrorCode ierr; 1216 Mat_MUMPS *mumps; 1217 PetscTruth isSeqAIJ; 1218 1219 PetscFunctionBegin; 1220 /* Create the factorization matrix */ 1221 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1222 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1223 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1224 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1225 if (isSeqAIJ) { 1226 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1227 } else { 1228 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1229 } 1230 1231 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1232 B->ops->view = MatView_MUMPS; 1233 B->ops->getinfo = MatGetInfo_MUMPS; 1234 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1235 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSetMumpsIcntl_C","MatSetMumpsIcntl",MatSetMumpsIcntl);CHKERRQ(ierr); 1236 if (ftype == MAT_FACTOR_LU) { 1237 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1238 B->factortype = MAT_FACTOR_LU; 1239 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1240 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1241 mumps->sym = 0; 1242 } else { 1243 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1244 B->factortype = MAT_FACTOR_CHOLESKY; 1245 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1246 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1247 if (A->spd_set && A->spd) mumps->sym = 1; 1248 else mumps->sym = 2; 1249 } 1250 1251 mumps->isAIJ = PETSC_TRUE; 1252 mumps->MatDestroy = B->ops->destroy; 1253 B->ops->destroy = MatDestroy_MUMPS; 1254 B->spptr = (void*)mumps; 1255 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1256 1257 *F = B; 1258 PetscFunctionReturn(0); 1259 } 1260 EXTERN_C_END 1261 1262 1263 EXTERN_C_BEGIN 1264 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1267 PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1268 { 1269 Mat B; 1270 PetscErrorCode ierr; 1271 Mat_MUMPS *mumps; 1272 PetscTruth isSeqSBAIJ; 1273 1274 PetscFunctionBegin; 1275 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1276 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"); 1277 ierr = PetscTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1278 /* Create the factorization matrix */ 1279 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1280 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1281 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1282 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1283 if (isSeqSBAIJ) { 1284 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1285 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1286 } else { 1287 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1288 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1289 } 1290 1291 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1292 B->ops->view = MatView_MUMPS; 1293 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1294 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSetMumpsIcntl_C","MatSetMumpsIcntl",MatSetMumpsIcntl);CHKERRQ(ierr); 1295 B->factortype = MAT_FACTOR_CHOLESKY; 1296 if (A->spd_set && A->spd) mumps->sym = 1; 1297 else mumps->sym = 2; 1298 1299 mumps->isAIJ = PETSC_FALSE; 1300 mumps->MatDestroy = B->ops->destroy; 1301 B->ops->destroy = MatDestroy_MUMPS; 1302 B->spptr = (void*)mumps; 1303 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 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 mumps->sym = 0; 1339 } else { 1340 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 1341 } 1342 1343 B->ops->view = MatView_MUMPS; 1344 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSetMumpsIcntl_C","MatSetMumpsIcntl",MatSetMumpsIcntl);CHKERRQ(ierr); 1346 1347 mumps->isAIJ = PETSC_TRUE; 1348 mumps->MatDestroy = B->ops->destroy; 1349 B->ops->destroy = MatDestroy_MUMPS; 1350 B->spptr = (void*)mumps; 1351 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1352 1353 *F = B; 1354 PetscFunctionReturn(0); 1355 } 1356 EXTERN_C_END 1357 1358 /* -------------------------------------------------------------------------------------------*/ 1359 #undef __FUNCT__ 1360 #define __FUNCT__ "MatSetMumpsIcntl" 1361 /*@ 1362 MatSetMumpsIcntl - Set MUMPS parameter ICNTL() 1363 1364 Logically Collective on Mat 1365 1366 Input Parameters: 1367 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1368 . icntl - index of MUMPS parameter array ICNTL() 1369 - ival - value of MUMPS ICNTL(icntl) 1370 1371 Options Database: 1372 . -mat_mumps_icntl_<icntl> <ival> 1373 1374 Level: beginner 1375 1376 References: MUMPS Users' Guide 1377 1378 .seealso: MatGetFactor() 1379 @*/ 1380 PetscErrorCode MatSetMumpsIcntl(Mat F,PetscInt icntl,PetscInt ival) 1381 { 1382 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 1383 1384 PetscFunctionBegin; 1385 PetscValidLogicalCollectiveInt(F,icntl,2); 1386 PetscValidLogicalCollectiveInt(F,ival,3); 1387 lu->id.ICNTL(icntl) = ival; 1388 PetscFunctionReturn(0); 1389 } 1390 1391