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