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