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 MatReuse reuse; 640 Mat F_diag; 641 PetscBool isMPIAIJ; 642 643 PetscFunctionBegin; 644 reuse = MAT_REUSE_MATRIX; 645 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 646 647 /* numerical factorization phase */ 648 /*-------------------------------*/ 649 lu->id.job = JOB_FACTNUMERIC; 650 if(!lu->id.ICNTL(18)) { 651 if (!lu->myid) { 652 #if defined(PETSC_USE_COMPLEX) 653 lu->id.a = (mumps_double_complex*)lu->val; 654 #else 655 lu->id.a = lu->val; 656 #endif 657 } 658 } else { 659 #if defined(PETSC_USE_COMPLEX) 660 lu->id.a_loc = (mumps_double_complex*)lu->val; 661 #else 662 lu->id.a_loc = lu->val; 663 #endif 664 } 665 #if defined(PETSC_USE_COMPLEX) 666 zmumps_c(&lu->id); 667 #else 668 dmumps_c(&lu->id); 669 #endif 670 if (lu->id.INFOG(1) < 0) { 671 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)); 672 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)); 673 } 674 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)); 675 676 if (lu->size > 1){ 677 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 678 if(isMPIAIJ) { 679 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 680 } else { 681 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 682 } 683 F_diag->assembled = PETSC_TRUE; 684 if (lu->nSolve){ 685 ierr = VecScatterDestroy(&lu->scat_sol);CHKERRQ(ierr); 686 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 687 ierr = VecDestroy(&lu->x_seq);CHKERRQ(ierr); 688 } 689 } 690 (F)->assembled = PETSC_TRUE; 691 lu->matstruc = SAME_NONZERO_PATTERN; 692 lu->CleanUpMUMPS = PETSC_TRUE; 693 lu->nSolve = 0; 694 695 if (lu->size > 1){ 696 /* distributed solution */ 697 if (!lu->nSolve){ 698 /* Create x_seq=sol_loc for repeated use */ 699 PetscInt lsol_loc; 700 PetscScalar *sol_loc; 701 lsol_loc = lu->id.INFO(23); /* length of sol_loc */ 702 ierr = PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&lu->id.isol_loc);CHKERRQ(ierr); 703 lu->id.lsol_loc = lsol_loc; 704 #if defined(PETSC_USE_COMPLEX) 705 lu->id.sol_loc = (mumps_double_complex*)sol_loc; 706 #else 707 lu->id.sol_loc = sol_loc; 708 #endif 709 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,lsol_loc,sol_loc,&lu->x_seq);CHKERRQ(ierr); 710 } 711 } 712 PetscFunctionReturn(0); 713 } 714 715 /* Sets MUMPS options from the options database */ 716 #undef __FUNCT__ 717 #define __FUNCT__ "PetscSetMUMPSFromOptions" 718 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 719 { 720 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 721 PetscErrorCode ierr; 722 PetscInt icntl; 723 PetscBool flg; 724 725 PetscFunctionBegin; 726 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 727 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 728 if (flg) mumps->id.ICNTL(1) = icntl; 729 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); 730 if (flg) mumps->id.ICNTL(2) = icntl; 731 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); 732 if (flg) mumps->id.ICNTL(3) = icntl; 733 734 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 735 if (flg) mumps->id.ICNTL(4) = icntl; 736 if (mumps->id.ICNTL(4) || PetscLogPrintInfo ) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 737 738 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); 739 if (flg) mumps->id.ICNTL(6) = icntl; 740 741 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); 742 if (flg) { 743 if (icntl== 1 && mumps->size > 1){ 744 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"); 745 } else { 746 mumps->id.ICNTL(7) = icntl; 747 } 748 } 749 750 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); 751 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); 752 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); 753 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); 754 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); 755 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); 756 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr); 757 758 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); 759 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); 760 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); 761 if (mumps->id.ICNTL(24)){ 762 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 763 } 764 765 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); 766 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); 767 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),PETSC_NULL);CHKERRQ(ierr); 768 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); 769 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); 770 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); 771 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); 772 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),PETSC_NULL);CHKERRQ(ierr); 773 774 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 775 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); 776 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 777 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); 778 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); 779 PetscOptionsEnd(); 780 PetscFunctionReturn(0); 781 } 782 783 #undef __FUNCT__ 784 #define __FUNCT__ "PetscInitializeMUMPS" 785 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS* mumps) 786 { 787 PetscErrorCode ierr; 788 789 PetscFunctionBegin; 790 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &mumps->myid); 791 ierr = MPI_Comm_size(((PetscObject)A)->comm,&mumps->size);CHKERRQ(ierr); 792 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(mumps->comm_mumps));CHKERRQ(ierr); 793 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 794 795 mumps->id.job = JOB_INIT; 796 mumps->id.par = 1; /* host participates factorizaton and solve */ 797 mumps->id.sym = mumps->sym; 798 #if defined(PETSC_USE_COMPLEX) 799 zmumps_c(&mumps->id); 800 #else 801 dmumps_c(&mumps->id); 802 #endif 803 804 mumps->CleanUpMUMPS = PETSC_FALSE; 805 mumps->scat_rhs = PETSC_NULL; 806 mumps->scat_sol = PETSC_NULL; 807 mumps->nSolve = 0; 808 809 /* set PETSc-MUMPS default options - override MUMPS default */ 810 mumps->id.ICNTL(3) = 0; 811 mumps->id.ICNTL(4) = 0; 812 if (mumps->size == 1){ 813 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 814 } else { 815 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 816 mumps->id.ICNTL(21) = 1; /* distributed solution */ 817 } 818 PetscFunctionReturn(0); 819 } 820 821 /* Note the Petsc r and c permutations are ignored */ 822 #undef __FUNCT__ 823 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 824 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 825 { 826 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 827 PetscErrorCode ierr; 828 MatReuse reuse; 829 Vec b; 830 IS is_iden; 831 const PetscInt M = A->rmap->N; 832 833 PetscFunctionBegin; 834 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 835 836 /* Set MUMPS options from the options database */ 837 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 838 839 reuse = MAT_INITIAL_MATRIX; 840 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 841 842 /* analysis phase */ 843 /*----------------*/ 844 lu->id.job = JOB_FACTSYMBOLIC; 845 lu->id.n = M; 846 switch (lu->id.ICNTL(18)){ 847 case 0: /* centralized assembled matrix input */ 848 if (!lu->myid) { 849 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 850 if (lu->id.ICNTL(6)>1){ 851 #if defined(PETSC_USE_COMPLEX) 852 lu->id.a = (mumps_double_complex*)lu->val; 853 #else 854 lu->id.a = lu->val; 855 #endif 856 } 857 if (lu->id.ICNTL(7) == 1){ /* use user-provide matrix ordering */ 858 if (!lu->myid) { 859 const PetscInt *idx; 860 PetscInt i,*perm_in; 861 ierr = PetscMalloc(M*sizeof(PetscInt),&perm_in);CHKERRQ(ierr); 862 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 863 lu->id.perm_in = perm_in; 864 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 865 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 866 } 867 } 868 } 869 break; 870 case 3: /* distributed assembled matrix input (size>1) */ 871 lu->id.nz_loc = lu->nz; 872 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 873 if (lu->id.ICNTL(6)>1) { 874 #if defined(PETSC_USE_COMPLEX) 875 lu->id.a_loc = (mumps_double_complex*)lu->val; 876 #else 877 lu->id.a_loc = lu->val; 878 #endif 879 } 880 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 881 if (!lu->myid){ 882 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 883 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 884 } else { 885 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 886 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 887 } 888 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 889 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 890 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 891 892 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 893 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 894 ierr = VecDestroy(&b);CHKERRQ(ierr); 895 break; 896 } 897 #if defined(PETSC_USE_COMPLEX) 898 zmumps_c(&lu->id); 899 #else 900 dmumps_c(&lu->id); 901 #endif 902 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)); 903 904 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 905 F->ops->solve = MatSolve_MUMPS; 906 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 907 F->ops->matsolve = MatMatSolve_MUMPS; 908 PetscFunctionReturn(0); 909 } 910 911 /* Note the Petsc r and c permutations are ignored */ 912 #undef __FUNCT__ 913 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 914 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 915 { 916 917 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 918 PetscErrorCode ierr; 919 MatReuse reuse; 920 Vec b; 921 IS is_iden; 922 const PetscInt M = A->rmap->N; 923 924 PetscFunctionBegin; 925 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 926 927 /* Set MUMPS options from the options database */ 928 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 929 930 reuse = MAT_INITIAL_MATRIX; 931 ierr = (*lu->ConvertToTriples)(A, 1, reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 932 933 /* analysis phase */ 934 /*----------------*/ 935 lu->id.job = JOB_FACTSYMBOLIC; 936 lu->id.n = M; 937 switch (lu->id.ICNTL(18)){ 938 case 0: /* centralized assembled matrix input */ 939 if (!lu->myid) { 940 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 941 if (lu->id.ICNTL(6)>1){ 942 #if defined(PETSC_USE_COMPLEX) 943 lu->id.a = (mumps_double_complex*)lu->val; 944 #else 945 lu->id.a = lu->val; 946 #endif 947 } 948 } 949 break; 950 case 3: /* distributed assembled matrix input (size>1) */ 951 lu->id.nz_loc = lu->nz; 952 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 953 if (lu->id.ICNTL(6)>1) { 954 #if defined(PETSC_USE_COMPLEX) 955 lu->id.a_loc = (mumps_double_complex*)lu->val; 956 #else 957 lu->id.a_loc = lu->val; 958 #endif 959 } 960 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 961 if (!lu->myid){ 962 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 963 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 964 } else { 965 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 966 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 967 } 968 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 969 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 970 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 971 972 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 973 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 974 ierr = VecDestroy(&b);CHKERRQ(ierr); 975 break; 976 } 977 #if defined(PETSC_USE_COMPLEX) 978 zmumps_c(&lu->id); 979 #else 980 dmumps_c(&lu->id); 981 #endif 982 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)); 983 984 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 985 F->ops->solve = MatSolve_MUMPS; 986 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 987 PetscFunctionReturn(0); 988 } 989 990 /* Note the Petsc r permutation and factor info are ignored */ 991 #undef __FUNCT__ 992 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 993 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 994 { 995 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 996 PetscErrorCode ierr; 997 MatReuse reuse; 998 Vec b; 999 IS is_iden; 1000 const PetscInt M = A->rmap->N; 1001 1002 PetscFunctionBegin; 1003 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 1004 1005 /* Set MUMPS options from the options database */ 1006 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1007 1008 reuse = MAT_INITIAL_MATRIX; 1009 ierr = (*lu->ConvertToTriples)(A, 1 , reuse, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 1010 1011 /* analysis phase */ 1012 /*----------------*/ 1013 lu->id.job = JOB_FACTSYMBOLIC; 1014 lu->id.n = M; 1015 switch (lu->id.ICNTL(18)){ 1016 case 0: /* centralized assembled matrix input */ 1017 if (!lu->myid) { 1018 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 1019 if (lu->id.ICNTL(6)>1){ 1020 #if defined(PETSC_USE_COMPLEX) 1021 lu->id.a = (mumps_double_complex*)lu->val; 1022 #else 1023 lu->id.a = lu->val; 1024 #endif 1025 } 1026 } 1027 break; 1028 case 3: /* distributed assembled matrix input (size>1) */ 1029 lu->id.nz_loc = lu->nz; 1030 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 1031 if (lu->id.ICNTL(6)>1) { 1032 #if defined(PETSC_USE_COMPLEX) 1033 lu->id.a_loc = (mumps_double_complex*)lu->val; 1034 #else 1035 lu->id.a_loc = lu->val; 1036 #endif 1037 } 1038 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1039 if (!lu->myid){ 1040 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 1041 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1042 } else { 1043 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 1044 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1045 } 1046 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 1047 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 1048 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 1049 1050 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 1051 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1052 ierr = VecDestroy(&b);CHKERRQ(ierr); 1053 break; 1054 } 1055 #if defined(PETSC_USE_COMPLEX) 1056 zmumps_c(&lu->id); 1057 #else 1058 dmumps_c(&lu->id); 1059 #endif 1060 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)); 1061 1062 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1063 F->ops->solve = MatSolve_MUMPS; 1064 F->ops->solvetranspose = MatSolve_MUMPS; 1065 #if !defined(PETSC_USE_COMPLEX) 1066 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1067 #else 1068 F->ops->getinertia = PETSC_NULL; 1069 #endif 1070 PetscFunctionReturn(0); 1071 } 1072 1073 #undef __FUNCT__ 1074 #define __FUNCT__ "MatView_MUMPS" 1075 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1076 { 1077 PetscErrorCode ierr; 1078 PetscBool iascii; 1079 PetscViewerFormat format; 1080 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 1081 1082 PetscFunctionBegin; 1083 /* check if matrix is mumps type */ 1084 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1085 1086 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1087 if (iascii) { 1088 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1089 if (format == PETSC_VIEWER_ASCII_INFO){ 1090 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1091 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",lu->id.sym);CHKERRQ(ierr); 1092 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",lu->id.par);CHKERRQ(ierr); 1093 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 1094 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 1095 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 1096 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 1097 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 1098 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 1099 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 1100 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 1101 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 1102 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 1103 if (lu->id.ICNTL(11)>0) { 1104 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 1105 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 1106 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 1107 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",lu->id.RINFOG(7),lu->id.RINFOG(8));CHKERRQ(ierr); 1108 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 1109 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 1110 } 1111 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 1112 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 1113 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 1114 /* ICNTL(15-17) not used */ 1115 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 1116 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 1117 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 1118 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 1119 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 1120 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 1121 1122 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 1123 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 1124 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 1125 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 1126 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",lu->id.ICNTL(28));CHKERRQ(ierr); 1127 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",lu->id.ICNTL(29));CHKERRQ(ierr); 1128 1129 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",lu->id.ICNTL(30));CHKERRQ(ierr); 1130 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",lu->id.ICNTL(31));CHKERRQ(ierr); 1131 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",lu->id.ICNTL(33));CHKERRQ(ierr); 1132 1133 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 1134 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 1135 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 1136 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 1137 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 1138 1139 /* infomation local to each processor */ 1140 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1141 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1142 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 1143 ierr = PetscViewerFlush(viewer); 1144 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1145 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 1146 ierr = PetscViewerFlush(viewer); 1147 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1148 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 1149 ierr = PetscViewerFlush(viewer); 1150 1151 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1152 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 1153 ierr = PetscViewerFlush(viewer); 1154 1155 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1156 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 1157 ierr = PetscViewerFlush(viewer); 1158 1159 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1160 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 1161 ierr = PetscViewerFlush(viewer); 1162 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1163 1164 if (!lu->myid){ /* information from the host */ 1165 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 1166 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 1167 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 1168 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); 1169 1170 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 1171 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 1172 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 1173 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 1174 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 1175 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 1176 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 1177 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 1178 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 1179 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 1180 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 1181 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 1182 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 1183 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); 1184 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); 1185 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); 1186 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); 1187 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 1188 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); 1189 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); 1190 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 1191 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 1192 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 1193 } 1194 } 1195 } 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "MatGetInfo_MUMPS" 1201 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1202 { 1203 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1204 1205 PetscFunctionBegin; 1206 info->block_size = 1.0; 1207 info->nz_allocated = mumps->id.INFOG(20); 1208 info->nz_used = mumps->id.INFOG(20); 1209 info->nz_unneeded = 0.0; 1210 info->assemblies = 0.0; 1211 info->mallocs = 0.0; 1212 info->memory = 0.0; 1213 info->fill_ratio_given = 0; 1214 info->fill_ratio_needed = 0; 1215 info->factor_mallocs = 0; 1216 PetscFunctionReturn(0); 1217 } 1218 1219 /* -------------------------------------------------------------------------------------------*/ 1220 #undef __FUNCT__ 1221 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1222 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1223 { 1224 Mat_MUMPS *lu =(Mat_MUMPS*)F->spptr; 1225 1226 PetscFunctionBegin; 1227 lu->id.ICNTL(icntl) = ival; 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "MatMumpsSetIcntl" 1233 /*@ 1234 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1235 1236 Logically Collective on Mat 1237 1238 Input Parameters: 1239 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1240 . icntl - index of MUMPS parameter array ICNTL() 1241 - ival - value of MUMPS ICNTL(icntl) 1242 1243 Options Database: 1244 . -mat_mumps_icntl_<icntl> <ival> 1245 1246 Level: beginner 1247 1248 References: MUMPS Users' Guide 1249 1250 .seealso: MatGetFactor() 1251 @*/ 1252 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1253 { 1254 PetscErrorCode ierr; 1255 1256 PetscFunctionBegin; 1257 PetscValidLogicalCollectiveInt(F,icntl,2); 1258 PetscValidLogicalCollectiveInt(F,ival,3); 1259 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1260 PetscFunctionReturn(0); 1261 } 1262 1263 /*MC 1264 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1265 distributed and sequential matrices via the external package MUMPS. 1266 1267 Works with MATAIJ and MATSBAIJ matrices 1268 1269 Options Database Keys: 1270 + -mat_mumps_icntl_4 <0,...,4> - print level 1271 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1272 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec) 1273 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1274 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1275 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1276 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1277 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1278 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1279 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1280 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1281 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1282 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1283 1284 Level: beginner 1285 1286 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1287 1288 M*/ 1289 1290 EXTERN_C_BEGIN 1291 #undef __FUNCT__ 1292 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1293 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1294 { 1295 PetscFunctionBegin; 1296 *type = MATSOLVERMUMPS; 1297 PetscFunctionReturn(0); 1298 } 1299 EXTERN_C_END 1300 1301 EXTERN_C_BEGIN 1302 /* MatGetFactor for Seq and MPI AIJ matrices */ 1303 #undef __FUNCT__ 1304 #define __FUNCT__ "MatGetFactor_aij_mumps" 1305 PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1306 { 1307 Mat B; 1308 PetscErrorCode ierr; 1309 Mat_MUMPS *mumps; 1310 PetscBool isSeqAIJ; 1311 1312 PetscFunctionBegin; 1313 /* Create the factorization matrix */ 1314 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1315 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1316 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1317 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1318 if (isSeqAIJ) { 1319 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1320 } else { 1321 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1322 } 1323 1324 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1325 B->ops->view = MatView_MUMPS; 1326 B->ops->getinfo = MatGetInfo_MUMPS; 1327 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1328 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl_MUMPS",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1329 if (ftype == MAT_FACTOR_LU) { 1330 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1331 B->factortype = MAT_FACTOR_LU; 1332 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1333 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1334 mumps->sym = 0; 1335 } else { 1336 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1337 B->factortype = MAT_FACTOR_CHOLESKY; 1338 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1339 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1340 if (A->spd_set && A->spd) mumps->sym = 1; 1341 else mumps->sym = 2; 1342 } 1343 1344 mumps->isAIJ = PETSC_TRUE; 1345 mumps->Destroy = B->ops->destroy; 1346 B->ops->destroy = MatDestroy_MUMPS; 1347 B->spptr = (void*)mumps; 1348 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1349 1350 *F = B; 1351 PetscFunctionReturn(0); 1352 } 1353 EXTERN_C_END 1354 1355 1356 EXTERN_C_BEGIN 1357 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1358 #undef __FUNCT__ 1359 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1360 PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1361 { 1362 Mat B; 1363 PetscErrorCode ierr; 1364 Mat_MUMPS *mumps; 1365 PetscBool isSeqSBAIJ; 1366 1367 PetscFunctionBegin; 1368 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1369 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"); 1370 ierr = PetscTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1371 /* Create the factorization matrix */ 1372 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1373 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1374 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1375 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1376 if (isSeqSBAIJ) { 1377 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1378 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1379 } else { 1380 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1381 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1382 } 1383 1384 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1385 B->ops->view = MatView_MUMPS; 1386 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1387 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1388 B->factortype = MAT_FACTOR_CHOLESKY; 1389 if (A->spd_set && A->spd) mumps->sym = 1; 1390 else mumps->sym = 2; 1391 1392 mumps->isAIJ = PETSC_FALSE; 1393 mumps->Destroy = B->ops->destroy; 1394 B->ops->destroy = MatDestroy_MUMPS; 1395 B->spptr = (void*)mumps; 1396 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1397 1398 *F = B; 1399 PetscFunctionReturn(0); 1400 } 1401 EXTERN_C_END 1402 1403 EXTERN_C_BEGIN 1404 #undef __FUNCT__ 1405 #define __FUNCT__ "MatGetFactor_baij_mumps" 1406 PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 1407 { 1408 Mat B; 1409 PetscErrorCode ierr; 1410 Mat_MUMPS *mumps; 1411 PetscBool isSeqBAIJ; 1412 1413 PetscFunctionBegin; 1414 /* Create the factorization matrix */ 1415 ierr = PetscTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1416 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1417 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1418 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1419 if (isSeqBAIJ) { 1420 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,PETSC_NULL);CHKERRQ(ierr); 1421 } else { 1422 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1423 } 1424 1425 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1426 if (ftype == MAT_FACTOR_LU) { 1427 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1428 B->factortype = MAT_FACTOR_LU; 1429 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 1430 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 1431 mumps->sym = 0; 1432 } else { 1433 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 1434 } 1435 1436 B->ops->view = MatView_MUMPS; 1437 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1438 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl_MUMPS",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1439 1440 mumps->isAIJ = PETSC_TRUE; 1441 mumps->Destroy = B->ops->destroy; 1442 B->ops->destroy = MatDestroy_MUMPS; 1443 B->spptr = (void*)mumps; 1444 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1445 1446 *F = B; 1447 PetscFunctionReturn(0); 1448 } 1449 EXTERN_C_END 1450 1451