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