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