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