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