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 #if defined(PETSC_USE_REAL_SINGLE) 12 #include <cmumps_c.h> 13 #else 14 #include <zmumps_c.h> 15 #endif 16 #else 17 #if defined(PETSC_USE_REAL_SINGLE) 18 #include <smumps_c.h> 19 #else 20 #include <dmumps_c.h> 21 #endif 22 #endif 23 EXTERN_C_END 24 #define JOB_INIT -1 25 #define JOB_FACTSYMBOLIC 1 26 #define JOB_FACTNUMERIC 2 27 #define JOB_SOLVE 3 28 #define JOB_END -2 29 30 /* calls to MUMPS */ 31 #if defined(PETSC_USE_COMPLEX) 32 #if defined(PETSC_USE_REAL_SINGLE) 33 #define PetscMUMPS_c cmumps_c 34 #else 35 #define PetscMUMPS_c zmumps_c 36 #endif 37 #else 38 #if defined(PETSC_USE_REAL_SINGLE) 39 #define PetscMUMPS_c smumps_c 40 #else 41 #define PetscMUMPS_c dmumps_c 42 #endif 43 #endif 44 45 46 /* macros s.t. indices match MUMPS documentation */ 47 #define ICNTL(I) icntl[(I)-1] 48 #define CNTL(I) cntl[(I)-1] 49 #define INFOG(I) infog[(I)-1] 50 #define INFO(I) info[(I)-1] 51 #define RINFOG(I) rinfog[(I)-1] 52 #define RINFO(I) rinfo[(I)-1] 53 54 typedef struct { 55 #if defined(PETSC_USE_COMPLEX) 56 #if defined(PETSC_USE_REAL_SINGLE) 57 CMUMPS_STRUC_C id; 58 #else 59 ZMUMPS_STRUC_C id; 60 #endif 61 #else 62 #if defined(PETSC_USE_REAL_SINGLE) 63 SMUMPS_STRUC_C id; 64 #else 65 DMUMPS_STRUC_C id; 66 #endif 67 #endif 68 69 MatStructure matstruc; 70 PetscMPIInt myid,size; 71 PetscInt *irn,*jcn,nz,sym; 72 PetscScalar *val; 73 MPI_Comm comm_mumps; 74 PetscBool isAIJ,CleanUpMUMPS; 75 PetscInt ICNTL9_pre; /* check if ICNTL(9) is changed from previous MatSolve */ 76 VecScatter scat_rhs, scat_sol; /* used by MatSolve() */ 77 Vec b_seq,x_seq; 78 VecScatter scat_rhss, scat_sols; /* used by MatMatSolve() */ 79 Vec bb_mpi,bb_seq,xx_mpi; 80 81 PetscErrorCode (*Destroy)(Mat); 82 PetscErrorCode (*ConvertToTriples)(Mat, int, MatReuse, int*, int**, int**, PetscScalar**); 83 } Mat_MUMPS; 84 85 extern PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*); 86 87 /* 88 MatConvertToTriples_A_B - convert Petsc matrix to triples: row[nz], col[nz], val[nz] 89 90 input: 91 A - matrix in aij,baij or sbaij (bs=1) format 92 shift - 0: C style output triple; 1: Fortran style output triple. 93 reuse - MAT_INITIAL_MATRIX: spaces are allocated and values are set for the triple 94 MAT_REUSE_MATRIX: only the values in v array are updated 95 output: 96 nnz - dim of r, c, and v (number of local nonzero entries of A) 97 r, c, v - row and col index, matrix values (matrix triples) 98 99 The returned values r, c, and sometimes v are obtained in a single PetscMalloc(). Then in MatDestroy_MUMPS() it is 100 freed with PetscFree((mumps->irn); This is not ideal code, the fact that v is ONLY sometimes part of mumps->irn means 101 that the PetscMalloc() cannot easily be replaced with a PetscMalloc3(). 102 103 */ 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "MatConvertToTriples_seqaij_seqaij" 107 PetscErrorCode MatConvertToTriples_seqaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 108 { 109 const PetscInt *ai,*aj,*ajj,M=A->rmap->n; 110 PetscInt nz,rnz,i,j; 111 PetscErrorCode ierr; 112 PetscInt *row,*col; 113 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 114 115 PetscFunctionBegin; 116 *v=aa->a; 117 if (reuse == MAT_INITIAL_MATRIX) { 118 nz = aa->nz; 119 ai = aa->i; 120 aj = aa->j; 121 *nnz = nz; 122 ierr = PetscMalloc1(2*nz, &row);CHKERRQ(ierr); 123 col = row + nz; 124 125 nz = 0; 126 for (i=0; i<M; i++) { 127 rnz = ai[i+1] - ai[i]; 128 ajj = aj + ai[i]; 129 for (j=0; j<rnz; j++) { 130 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 131 } 132 } 133 *r = row; *c = col; 134 } 135 PetscFunctionReturn(0); 136 } 137 138 #undef __FUNCT__ 139 #define __FUNCT__ "MatConvertToTriples_seqbaij_seqaij" 140 PetscErrorCode MatConvertToTriples_seqbaij_seqaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 141 { 142 Mat_SeqBAIJ *aa=(Mat_SeqBAIJ*)A->data; 143 const PetscInt *ai,*aj,*ajj,bs2 = aa->bs2; 144 PetscInt bs,M,nz,idx=0,rnz,i,j,k,m; 145 PetscErrorCode ierr; 146 PetscInt *row,*col; 147 148 PetscFunctionBegin; 149 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 150 M = A->rmap->N/bs; 151 *v = aa->a; 152 if (reuse == MAT_INITIAL_MATRIX) { 153 ai = aa->i; aj = aa->j; 154 nz = bs2*aa->nz; 155 *nnz = nz; 156 ierr = PetscMalloc1(2*nz, &row);CHKERRQ(ierr); 157 col = row + nz; 158 159 for (i=0; i<M; i++) { 160 ajj = aj + ai[i]; 161 rnz = ai[i+1] - ai[i]; 162 for (k=0; k<rnz; k++) { 163 for (j=0; j<bs; j++) { 164 for (m=0; m<bs; m++) { 165 row[idx] = i*bs + m + shift; 166 col[idx++] = bs*(ajj[k]) + j + shift; 167 } 168 } 169 } 170 } 171 *r = row; *c = col; 172 } 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "MatConvertToTriples_seqsbaij_seqsbaij" 178 PetscErrorCode MatConvertToTriples_seqsbaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 179 { 180 const PetscInt *ai, *aj,*ajj,M=A->rmap->n; 181 PetscInt nz,rnz,i,j; 182 PetscErrorCode ierr; 183 PetscInt *row,*col; 184 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data; 185 186 PetscFunctionBegin; 187 *v = aa->a; 188 if (reuse == MAT_INITIAL_MATRIX) { 189 nz = aa->nz; 190 ai = aa->i; 191 aj = aa->j; 192 *v = aa->a; 193 *nnz = nz; 194 ierr = PetscMalloc1(2*nz, &row);CHKERRQ(ierr); 195 col = row + nz; 196 197 nz = 0; 198 for (i=0; i<M; i++) { 199 rnz = ai[i+1] - ai[i]; 200 ajj = aj + ai[i]; 201 for (j=0; j<rnz; j++) { 202 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 203 } 204 } 205 *r = row; *c = col; 206 } 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "MatConvertToTriples_seqaij_seqsbaij" 212 PetscErrorCode MatConvertToTriples_seqaij_seqsbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 213 { 214 const PetscInt *ai,*aj,*ajj,*adiag,M=A->rmap->n; 215 PetscInt nz,rnz,i,j; 216 const PetscScalar *av,*v1; 217 PetscScalar *val; 218 PetscErrorCode ierr; 219 PetscInt *row,*col; 220 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 221 222 PetscFunctionBegin; 223 ai =aa->i; aj=aa->j;av=aa->a; 224 adiag=aa->diag; 225 if (reuse == MAT_INITIAL_MATRIX) { 226 /* count nz in the uppper triangular part of A */ 227 nz = 0; 228 for (i=0; i<M; i++) nz += ai[i+1] - adiag[i]; 229 *nnz = nz; 230 231 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 232 col = row + nz; 233 val = (PetscScalar*)(col + nz); 234 235 nz = 0; 236 for (i=0; i<M; i++) { 237 rnz = ai[i+1] - adiag[i]; 238 ajj = aj + adiag[i]; 239 v1 = av + adiag[i]; 240 for (j=0; j<rnz; j++) { 241 row[nz] = i+shift; col[nz] = ajj[j] + shift; val[nz++] = v1[j]; 242 } 243 } 244 *r = row; *c = col; *v = val; 245 } else { 246 nz = 0; val = *v; 247 for (i=0; i <M; i++) { 248 rnz = ai[i+1] - adiag[i]; 249 ajj = aj + adiag[i]; 250 v1 = av + adiag[i]; 251 for (j=0; j<rnz; j++) { 252 val[nz++] = v1[j]; 253 } 254 } 255 } 256 PetscFunctionReturn(0); 257 } 258 259 #undef __FUNCT__ 260 #define __FUNCT__ "MatConvertToTriples_mpisbaij_mpisbaij" 261 PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 262 { 263 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 264 PetscErrorCode ierr; 265 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 266 PetscInt *row,*col; 267 const PetscScalar *av, *bv,*v1,*v2; 268 PetscScalar *val; 269 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 270 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)(mat->A)->data; 271 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 272 273 PetscFunctionBegin; 274 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 275 av=aa->a; bv=bb->a; 276 277 garray = mat->garray; 278 279 if (reuse == MAT_INITIAL_MATRIX) { 280 nz = aa->nz + bb->nz; 281 *nnz = nz; 282 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 283 col = row + nz; 284 val = (PetscScalar*)(col + nz); 285 286 *r = row; *c = col; *v = val; 287 } else { 288 row = *r; col = *c; val = *v; 289 } 290 291 jj = 0; irow = rstart; 292 for (i=0; i<m; i++) { 293 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 294 countA = ai[i+1] - ai[i]; 295 countB = bi[i+1] - bi[i]; 296 bjj = bj + bi[i]; 297 v1 = av + ai[i]; 298 v2 = bv + bi[i]; 299 300 /* A-part */ 301 for (j=0; j<countA; j++) { 302 if (reuse == MAT_INITIAL_MATRIX) { 303 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 304 } 305 val[jj++] = v1[j]; 306 } 307 308 /* B-part */ 309 for (j=0; j < countB; j++) { 310 if (reuse == MAT_INITIAL_MATRIX) { 311 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 312 } 313 val[jj++] = v2[j]; 314 } 315 irow++; 316 } 317 PetscFunctionReturn(0); 318 } 319 320 #undef __FUNCT__ 321 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpiaij" 322 PetscErrorCode MatConvertToTriples_mpiaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 323 { 324 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 325 PetscErrorCode ierr; 326 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 327 PetscInt *row,*col; 328 const PetscScalar *av, *bv,*v1,*v2; 329 PetscScalar *val; 330 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 331 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(mat->A)->data; 332 Mat_SeqAIJ *bb = (Mat_SeqAIJ*)(mat->B)->data; 333 334 PetscFunctionBegin; 335 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 336 av=aa->a; bv=bb->a; 337 338 garray = mat->garray; 339 340 if (reuse == MAT_INITIAL_MATRIX) { 341 nz = aa->nz + bb->nz; 342 *nnz = nz; 343 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 344 col = row + nz; 345 val = (PetscScalar*)(col + nz); 346 347 *r = row; *c = col; *v = val; 348 } else { 349 row = *r; col = *c; val = *v; 350 } 351 352 jj = 0; irow = rstart; 353 for (i=0; i<m; i++) { 354 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 355 countA = ai[i+1] - ai[i]; 356 countB = bi[i+1] - bi[i]; 357 bjj = bj + bi[i]; 358 v1 = av + ai[i]; 359 v2 = bv + bi[i]; 360 361 /* A-part */ 362 for (j=0; j<countA; j++) { 363 if (reuse == MAT_INITIAL_MATRIX) { 364 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 365 } 366 val[jj++] = v1[j]; 367 } 368 369 /* B-part */ 370 for (j=0; j < countB; j++) { 371 if (reuse == MAT_INITIAL_MATRIX) { 372 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 373 } 374 val[jj++] = v2[j]; 375 } 376 irow++; 377 } 378 PetscFunctionReturn(0); 379 } 380 381 #undef __FUNCT__ 382 #define __FUNCT__ "MatConvertToTriples_mpibaij_mpiaij" 383 PetscErrorCode MatConvertToTriples_mpibaij_mpiaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 384 { 385 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data; 386 Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)(mat->A)->data; 387 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 388 const PetscInt *ai = aa->i, *bi = bb->i, *aj = aa->j, *bj = bb->j,*ajj, *bjj; 389 const PetscInt *garray = mat->garray,mbs=mat->mbs,rstart=A->rmap->rstart; 390 const PetscInt bs2=mat->bs2; 391 PetscErrorCode ierr; 392 PetscInt bs,nz,i,j,k,n,jj,irow,countA,countB,idx; 393 PetscInt *row,*col; 394 const PetscScalar *av=aa->a, *bv=bb->a,*v1,*v2; 395 PetscScalar *val; 396 397 PetscFunctionBegin; 398 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 399 if (reuse == MAT_INITIAL_MATRIX) { 400 nz = bs2*(aa->nz + bb->nz); 401 *nnz = nz; 402 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 403 col = row + nz; 404 val = (PetscScalar*)(col + nz); 405 406 *r = row; *c = col; *v = val; 407 } else { 408 row = *r; col = *c; val = *v; 409 } 410 411 jj = 0; irow = rstart; 412 for (i=0; i<mbs; i++) { 413 countA = ai[i+1] - ai[i]; 414 countB = bi[i+1] - bi[i]; 415 ajj = aj + ai[i]; 416 bjj = bj + bi[i]; 417 v1 = av + bs2*ai[i]; 418 v2 = bv + bs2*bi[i]; 419 420 idx = 0; 421 /* A-part */ 422 for (k=0; k<countA; k++) { 423 for (j=0; j<bs; j++) { 424 for (n=0; n<bs; n++) { 425 if (reuse == MAT_INITIAL_MATRIX) { 426 row[jj] = irow + n + shift; 427 col[jj] = rstart + bs*ajj[k] + j + shift; 428 } 429 val[jj++] = v1[idx++]; 430 } 431 } 432 } 433 434 idx = 0; 435 /* B-part */ 436 for (k=0; k<countB; k++) { 437 for (j=0; j<bs; j++) { 438 for (n=0; n<bs; n++) { 439 if (reuse == MAT_INITIAL_MATRIX) { 440 row[jj] = irow + n + shift; 441 col[jj] = bs*garray[bjj[k]] + j + shift; 442 } 443 val[jj++] = v2[idx++]; 444 } 445 } 446 } 447 irow += bs; 448 } 449 PetscFunctionReturn(0); 450 } 451 452 #undef __FUNCT__ 453 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpisbaij" 454 PetscErrorCode MatConvertToTriples_mpiaij_mpisbaij(Mat A,int shift,MatReuse reuse,int *nnz,int **r, int **c, PetscScalar **v) 455 { 456 const PetscInt *ai, *aj,*adiag, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 457 PetscErrorCode ierr; 458 PetscInt rstart,nz,nza,nzb,i,j,jj,irow,countA,countB; 459 PetscInt *row,*col; 460 const PetscScalar *av, *bv,*v1,*v2; 461 PetscScalar *val; 462 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 463 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(mat->A)->data; 464 Mat_SeqAIJ *bb =(Mat_SeqAIJ*)(mat->B)->data; 465 466 PetscFunctionBegin; 467 ai=aa->i; aj=aa->j; adiag=aa->diag; 468 bi=bb->i; bj=bb->j; garray = mat->garray; 469 av=aa->a; bv=bb->a; 470 471 rstart = A->rmap->rstart; 472 473 if (reuse == MAT_INITIAL_MATRIX) { 474 nza = 0; /* num of upper triangular entries in mat->A, including diagonals */ 475 nzb = 0; /* num of upper triangular entries in mat->B */ 476 for (i=0; i<m; i++) { 477 nza += (ai[i+1] - adiag[i]); 478 countB = bi[i+1] - bi[i]; 479 bjj = bj + bi[i]; 480 for (j=0; j<countB; j++) { 481 if (garray[bjj[j]] > rstart) nzb++; 482 } 483 } 484 485 nz = nza + nzb; /* total nz of upper triangular part of mat */ 486 *nnz = nz; 487 ierr = PetscMalloc((2*nz*sizeof(PetscInt)+nz*sizeof(PetscScalar)), &row);CHKERRQ(ierr); 488 col = row + nz; 489 val = (PetscScalar*)(col + nz); 490 491 *r = row; *c = col; *v = val; 492 } else { 493 row = *r; col = *c; val = *v; 494 } 495 496 jj = 0; irow = rstart; 497 for (i=0; i<m; i++) { 498 ajj = aj + adiag[i]; /* ptr to the beginning of the diagonal of this row */ 499 v1 = av + adiag[i]; 500 countA = ai[i+1] - adiag[i]; 501 countB = bi[i+1] - bi[i]; 502 bjj = bj + bi[i]; 503 v2 = bv + bi[i]; 504 505 /* A-part */ 506 for (j=0; j<countA; j++) { 507 if (reuse == MAT_INITIAL_MATRIX) { 508 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 509 } 510 val[jj++] = v1[j]; 511 } 512 513 /* B-part */ 514 for (j=0; j < countB; j++) { 515 if (garray[bjj[j]] > rstart) { 516 if (reuse == MAT_INITIAL_MATRIX) { 517 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 518 } 519 val[jj++] = v2[j]; 520 } 521 } 522 irow++; 523 } 524 PetscFunctionReturn(0); 525 } 526 527 #undef __FUNCT__ 528 #define __FUNCT__ "MatGetDiagonal_MUMPS" 529 PetscErrorCode MatGetDiagonal_MUMPS(Mat A,Vec v) 530 { 531 PetscFunctionBegin; 532 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: MUMPS factor"); 533 PetscFunctionReturn(0); 534 } 535 536 #undef __FUNCT__ 537 #define __FUNCT__ "MatDestroy_MUMPS" 538 PetscErrorCode MatDestroy_MUMPS(Mat A) 539 { 540 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 541 PetscErrorCode ierr; 542 543 PetscFunctionBegin; 544 if (mumps->CleanUpMUMPS) { 545 /* Terminate instance, deallocate memories */ 546 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 547 ierr = VecScatterDestroy(&mumps->scat_rhs);CHKERRQ(ierr); 548 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 549 ierr = VecDestroy(&mumps->b_seq);CHKERRQ(ierr); 550 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 551 ierr = PetscFree(mumps->id.perm_in);CHKERRQ(ierr); 552 ierr = PetscFree(mumps->irn);CHKERRQ(ierr); 553 554 /* created in MatMatSolve_MUMPS() */ 555 if (mumps->scat_rhss) { 556 ierr = VecScatterDestroy(&mumps->scat_rhss);CHKERRQ(ierr); 557 ierr = VecScatterDestroy(&mumps->scat_sols);CHKERRQ(ierr); 558 ierr = VecDestroy(&mumps->bb_seq);CHKERRQ(ierr); 559 ierr = VecDestroy(&mumps->bb_mpi);CHKERRQ(ierr); 560 ierr = VecDestroy(&mumps->xx_mpi);CHKERRQ(ierr); 561 } 562 563 mumps->id.job = JOB_END; 564 PetscMUMPS_c(&mumps->id); 565 ierr = MPI_Comm_free(&(mumps->comm_mumps));CHKERRQ(ierr); 566 } 567 if (mumps->Destroy) { 568 ierr = (mumps->Destroy)(A);CHKERRQ(ierr); 569 } 570 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 571 572 /* clear composed functions */ 573 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);CHKERRQ(ierr); 574 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetIcntl_C",NULL);CHKERRQ(ierr); 575 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetIcntl_C",NULL);CHKERRQ(ierr); 576 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsSetCntl_C",NULL);CHKERRQ(ierr); 577 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetCntl_C",NULL);CHKERRQ(ierr); 578 579 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfo_C",NULL);CHKERRQ(ierr); 580 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetInfog_C",NULL);CHKERRQ(ierr); 581 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfo_C",NULL);CHKERRQ(ierr); 582 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMumpsGetRinfog_C",NULL);CHKERRQ(ierr); 583 PetscFunctionReturn(0); 584 } 585 586 #undef __FUNCT__ 587 #define __FUNCT__ "MatSolve_MUMPS" 588 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 589 { 590 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 591 PetscScalar *array; 592 Vec b_seq; 593 IS is_iden,is_petsc; 594 PetscErrorCode ierr; 595 PetscInt i; 596 static PetscBool cite1 = PETSC_FALSE,cite2 = PETSC_FALSE; 597 598 PetscFunctionBegin; 599 ierr = PetscCitationsRegister("@article{MUMPS01,\n author = {P.~R. Amestoy and I.~S. Duff and J.-Y. L'Excellent and J. Koster},\n title = {A fully asynchronous multifrontal solver using distributed dynamic scheduling},\n journal = {SIAM Journal on Matrix Analysis and Applications},\n volume = {23},\n number = {1},\n pages = {15--41},\n year = {2001}\n}\n",&cite1);CHKERRQ(ierr); 600 ierr = PetscCitationsRegister("@article{MUMPS02,\n author = {P.~R. Amestoy and A. Guermouche and J.-Y. L'Excellent and S. Pralet},\n title = {Hybrid scheduling for the parallel solution of linear systems},\n journal = {Parallel Computing},\n volume = {32},\n number = {2},\n pages = {136--156},\n year = {2006}\n}\n",&cite2);CHKERRQ(ierr); 601 mumps->id.nrhs = 1; 602 b_seq = mumps->b_seq; 603 if (mumps->size > 1) { 604 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 605 ierr = VecScatterBegin(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 606 ierr = VecScatterEnd(mumps->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 607 if (!mumps->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 608 } else { /* size == 1 */ 609 ierr = VecCopy(b,x);CHKERRQ(ierr); 610 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 611 } 612 if (!mumps->myid) { /* define rhs on the host */ 613 mumps->id.nrhs = 1; 614 #if defined(PETSC_USE_COMPLEX) 615 #if defined(PETSC_USE_REAL_SINGLE) 616 mumps->id.rhs = (mumps_complex*)array; 617 #else 618 mumps->id.rhs = (mumps_double_complex*)array; 619 #endif 620 #else 621 mumps->id.rhs = array; 622 #endif 623 } 624 625 /* solve phase */ 626 /*-------------*/ 627 mumps->id.job = JOB_SOLVE; 628 PetscMUMPS_c(&mumps->id); 629 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 630 631 if (mumps->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 632 if (mumps->scat_sol && mumps->ICNTL9_pre != mumps->id.ICNTL(9)) { 633 /* when id.ICNTL(9) changes, the contents of lsol_loc may change (not its size, lsol_loc), recreates scat_sol */ 634 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 635 } 636 if (!mumps->scat_sol) { /* create scatter scat_sol */ 637 ierr = ISCreateStride(PETSC_COMM_SELF,mumps->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 638 for (i=0; i<mumps->id.lsol_loc; i++) { 639 mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 640 } 641 ierr = ISCreateGeneral(PETSC_COMM_SELF,mumps->id.lsol_loc,mumps->id.isol_loc,PETSC_COPY_VALUES,&is_petsc);CHKERRQ(ierr); /* to */ 642 ierr = VecScatterCreate(mumps->x_seq,is_iden,x,is_petsc,&mumps->scat_sol);CHKERRQ(ierr); 643 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 644 ierr = ISDestroy(&is_petsc);CHKERRQ(ierr); 645 646 mumps->ICNTL9_pre = mumps->id.ICNTL(9); /* save current value of id.ICNTL(9) */ 647 } 648 649 ierr = VecScatterBegin(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 650 ierr = VecScatterEnd(mumps->scat_sol,mumps->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 651 } 652 PetscFunctionReturn(0); 653 } 654 655 #undef __FUNCT__ 656 #define __FUNCT__ "MatSolveTranspose_MUMPS" 657 PetscErrorCode MatSolveTranspose_MUMPS(Mat A,Vec b,Vec x) 658 { 659 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 660 PetscErrorCode ierr; 661 662 PetscFunctionBegin; 663 mumps->id.ICNTL(9) = 0; 664 ierr = MatSolve_MUMPS(A,b,x);CHKERRQ(ierr); 665 mumps->id.ICNTL(9) = 1; 666 PetscFunctionReturn(0); 667 } 668 669 #undef __FUNCT__ 670 #define __FUNCT__ "MatMatSolve_MUMPS" 671 PetscErrorCode MatMatSolve_MUMPS(Mat A,Mat B,Mat X) 672 { 673 PetscErrorCode ierr; 674 PetscBool flg; 675 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 676 PetscInt i,nrhs,m,M,mx; 677 PetscScalar *array,*bray; 678 679 PetscFunctionBegin; 680 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 681 if (!flg) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 682 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 683 if (!flg) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 684 if (B->rmap->n != X->rmap->n) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONG,"Matrix B and X must have same row distribution"); 685 686 ierr = MatGetLocalSize(B,&m,NULL);CHKERRQ(ierr); 687 ierr = MatGetSize(B,&M,&nrhs);CHKERRQ(ierr); 688 689 if (mumps->size == 1) { 690 /* copy B to X */ 691 ierr = MatDenseGetArray(B,&bray);CHKERRQ(ierr); 692 ierr = MatDenseGetArray(X,&array);CHKERRQ(ierr); 693 for (i=0; i<M*nrhs; i++) array[i] = bray[i]; 694 ierr = MatDenseRestoreArray(B,&bray);CHKERRQ(ierr); 695 696 mumps->id.nrhs = nrhs; 697 mumps->id.lrhs = M; 698 #if defined(PETSC_USE_COMPLEX) 699 #if defined(PETSC_USE_REAL_SINGLE) 700 mumps->id.rhs = (mumps_complex*)array; 701 #else 702 mumps->id.rhs = (mumps_double_complex*)array; 703 #endif 704 #else 705 mumps->id.rhs = array; 706 #endif 707 708 /* solve phase */ 709 /*-------------*/ 710 mumps->id.job = JOB_SOLVE; 711 PetscMUMPS_c(&mumps->id); 712 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 713 714 ierr = MatDenseRestoreArray(X,&array);CHKERRQ(ierr); 715 } else { /************** parallel case ***************/ 716 PetscInt lsol_loc,*isol_loc,*isol_loc_save,*idx,*iidx,*idxx; 717 PetscScalar *sol_loc,*sol_loc_save; 718 IS is_to,is_from; 719 PetscInt k,proc,j; 720 const PetscInt *rstart; 721 722 /* create x_seq to hold local solution */ 723 lsol_loc = nrhs*mumps->id.INFO(23); /* length of sol_loc */ 724 725 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); //save it for MatSovle()!!! 726 sol_loc_save = mumps->id.sol_loc; /* these arrays were allocated for MatSolve_MUMPS() */ 727 isol_loc_save = mumps->id.isol_loc; 728 729 ierr = PetscMalloc2(lsol_loc,&sol_loc,lsol_loc,&isol_loc);CHKERRQ(ierr); 730 #if defined(PETSC_USE_COMPLEX) 731 #if defined(PETSC_USE_REAL_SINGLE) 732 mumps->id.sol_loc = (mumps_complex*)sol_loc; 733 #else 734 mumps->id.sol_loc = (mumps_double_complex*)sol_loc; 735 #endif 736 #else 737 mumps->id.sol_loc = sol_loc; 738 #endif 739 mumps->id.isol_loc = isol_loc; 740 741 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 742 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 743 744 /* copy rhs matrix B into vector bb_mpi */ 745 if (mumps->bb_mpi) { 746 ierr = VecDestroy(&mumps->bb_mpi);CHKERRQ(ierr); 747 ierr = VecDestroy(&mumps->bb_seq);CHKERRQ(ierr); 748 ierr = VecDestroy(&mumps->xx_mpi);CHKERRQ(ierr); 749 ierr = VecScatterDestroy(&mumps->scat_rhss);CHKERRQ(ierr); 750 ierr = VecScatterDestroy(&mumps->scat_sols);CHKERRQ(ierr); 751 } 752 ierr = MatDenseGetArray(B,&bray);CHKERRQ(ierr); 753 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)B),1,nrhs*m,nrhs*M,(const PetscScalar*)bray,&mumps->bb_mpi);CHKERRQ(ierr); 754 ierr = MatDenseRestoreArray(B,&bray);CHKERRQ(ierr); 755 756 ierr = MatDenseGetArray(X,&array);CHKERRQ(ierr); 757 ierr = MatGetLocalSize(B,&mx,NULL);CHKERRQ(ierr); 758 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)X),1,nrhs*mx,nrhs*M,(const PetscScalar*)array,&mumps->xx_mpi);CHKERRQ(ierr); 759 760 /* scatter bb_mpi to bb_seq because MUMPS only supports centralized rhs */ 761 /* idx: maps from k-th index of bb_mpi to (i,j)-th global entry of B; 762 iidx: inverse of idx, will be used by scattering xx_seq -> X */ 763 ierr = PetscMalloc2(nrhs*M,&idx,nrhs*M,&iidx);CHKERRQ(ierr); 764 ierr = MatGetOwnershipRanges(B,&rstart);CHKERRQ(ierr); 765 k = 0; 766 for (proc=0; proc<mumps->size; proc++){ 767 for (j=0; j<nrhs; j++){ 768 for (i=rstart[proc]; i<rstart[proc+1]; i++){ 769 iidx[j*M + i] = k; 770 idx[k++] = j*M + i; 771 } 772 } 773 } 774 775 if (!mumps->myid) { 776 ierr = VecCreateSeq(PETSC_COMM_SELF,nrhs*M,&mumps->bb_seq);CHKERRQ(ierr); 777 ierr = ISCreateGeneral(PETSC_COMM_SELF,nrhs*M,idx,PETSC_COPY_VALUES,&is_to);CHKERRQ(ierr); 778 ierr = ISCreateStride(PETSC_COMM_SELF,nrhs*M,0,1,&is_from);CHKERRQ(ierr); 779 } else { 780 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->bb_seq);CHKERRQ(ierr); 781 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_to);CHKERRQ(ierr); 782 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_from);CHKERRQ(ierr); 783 } 784 ierr = VecScatterCreate(mumps->bb_mpi,is_from,mumps->bb_seq,is_to,&mumps->scat_rhss);CHKERRQ(ierr); 785 ierr = VecScatterBegin(mumps->scat_rhss,mumps->bb_mpi,mumps->bb_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 786 ierr = ISDestroy(&is_to);CHKERRQ(ierr); 787 ierr = ISDestroy(&is_from);CHKERRQ(ierr); 788 ierr = VecScatterEnd(mumps->scat_rhss,mumps->bb_mpi,mumps->bb_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 789 790 if (!mumps->myid) { /* define rhs on the host */ 791 ierr = VecGetArray(mumps->bb_seq,&bray);CHKERRQ(ierr); 792 mumps->id.nrhs = nrhs; 793 mumps->id.lrhs = M; 794 #if defined(PETSC_USE_COMPLEX) 795 #if defined(PETSC_USE_REAL_SINGLE) 796 mumps->id.rhs = (mumps_complex*)bray; 797 #else 798 mumps->id.rhs = (mumps_double_complex*)bray; 799 #endif 800 #else 801 mumps->id.rhs = bray; 802 #endif 803 ierr = VecRestoreArray(mumps->bb_seq,&bray);CHKERRQ(ierr); 804 } 805 806 /* solve phase */ 807 /*-------------*/ 808 mumps->id.job = JOB_SOLVE; 809 PetscMUMPS_c(&mumps->id); 810 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 811 812 /* put mumps distributed solution to petsc vector xx_mpi, which shares local arrays with solution matrix X */ 813 814 /* create scatter scat_sols */ 815 ierr = PetscMalloc1(nrhs*mumps->id.lsol_loc,&idxx);CHKERRQ(ierr); 816 ierr = ISCreateStride(PETSC_COMM_SELF,nrhs*mumps->id.lsol_loc,0,1,&is_from);CHKERRQ(ierr); 817 for (i=0; i<mumps->id.lsol_loc; i++) { 818 mumps->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 819 idxx[i] = iidx[mumps->id.isol_loc[i]]; 820 for (j=1; j<nrhs; j++){ 821 idxx[j*mumps->id.lsol_loc+i] = iidx[mumps->id.isol_loc[i]+j*M]; 822 } 823 } 824 ierr = ISCreateGeneral(PETSC_COMM_SELF,nrhs*mumps->id.lsol_loc,idxx,PETSC_COPY_VALUES,&is_to);CHKERRQ(ierr); 825 ierr = VecScatterCreate(mumps->x_seq,is_from,mumps->xx_mpi,is_to,&mumps->scat_sols);CHKERRQ(ierr); 826 ierr = VecScatterBegin(mumps->scat_sols,mumps->x_seq,mumps->xx_mpi,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 827 ierr = ISDestroy(&is_from);CHKERRQ(ierr); 828 ierr = ISDestroy(&is_to);CHKERRQ(ierr); 829 ierr = VecScatterEnd(mumps->scat_sols,mumps->x_seq,mumps->xx_mpi,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 830 831 ierr = MatDenseRestoreArray(X,&array);CHKERRQ(ierr); 832 //ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 833 //mumps->id.sol_loc = sol_loc_save; 834 //mumps->id.isol_loc = isol_loc_save; 835 //ierr = PetscFree2(sol_loc,isol_loc);CHKERRQ(ierr); 836 ierr = PetscFree2(idx,iidx);CHKERRQ(ierr); 837 ierr = PetscFree(idxx);CHKERRQ(ierr); 838 } 839 PetscFunctionReturn(0); 840 } 841 842 #if !defined(PETSC_USE_COMPLEX) 843 /* 844 input: 845 F: numeric factor 846 output: 847 nneg: total number of negative pivots 848 nzero: 0 849 npos: (global dimension of F) - nneg 850 */ 851 852 #undef __FUNCT__ 853 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 854 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 855 { 856 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 857 PetscErrorCode ierr; 858 PetscMPIInt size; 859 860 PetscFunctionBegin; 861 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);CHKERRQ(ierr); 862 /* 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 */ 863 if (size > 1 && mumps->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",mumps->id.INFOG(13)); 864 865 if (nneg) *nneg = mumps->id.INFOG(12); 866 if (nzero || npos) { 867 if (mumps->id.ICNTL(24) != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"-mat_mumps_icntl_24 must be set as 1 for null pivot row detection"); 868 if (nzero) *nzero = mumps->id.INFOG(28); 869 if (npos) *npos = F->rmap->N - (mumps->id.INFOG(12) + mumps->id.INFOG(28)); 870 } 871 PetscFunctionReturn(0); 872 } 873 #endif /* !defined(PETSC_USE_COMPLEX) */ 874 875 #undef __FUNCT__ 876 #define __FUNCT__ "MatFactorNumeric_MUMPS" 877 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 878 { 879 Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr; 880 PetscErrorCode ierr; 881 Mat F_diag; 882 PetscBool isMPIAIJ; 883 884 PetscFunctionBegin; 885 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 886 887 /* numerical factorization phase */ 888 /*-------------------------------*/ 889 mumps->id.job = JOB_FACTNUMERIC; 890 if (!mumps->id.ICNTL(18)) { /* A is centralized */ 891 if (!mumps->myid) { 892 #if defined(PETSC_USE_COMPLEX) 893 #if defined(PETSC_USE_REAL_SINGLE) 894 mumps->id.a = (mumps_complex*)mumps->val; 895 #else 896 mumps->id.a = (mumps_double_complex*)mumps->val; 897 #endif 898 #else 899 mumps->id.a = mumps->val; 900 #endif 901 } 902 } else { 903 #if defined(PETSC_USE_COMPLEX) 904 #if defined(PETSC_USE_REAL_SINGLE) 905 mumps->id.a_loc = (mumps_complex*)mumps->val; 906 #else 907 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 908 #endif 909 #else 910 mumps->id.a_loc = mumps->val; 911 #endif 912 } 913 PetscMUMPS_c(&mumps->id); 914 if (mumps->id.INFOG(1) < 0) { 915 if (mumps->id.INFO(1) == -13) { 916 if (mumps->id.INFO(2) < 0) { 917 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d megabytes\n",-mumps->id.INFO(2)); 918 } else { 919 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d bytes\n",mumps->id.INFO(2)); 920 } 921 } else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: INFO(1)=%d, INFO(2)=%d\n",mumps->id.INFO(1),mumps->id.INFO(2)); 922 } 923 if (!mumps->myid && mumps->id.ICNTL(16) > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB," mumps->id.ICNTL(16):=%d\n",mumps->id.INFOG(16)); 924 925 (F)->assembled = PETSC_TRUE; 926 mumps->matstruc = SAME_NONZERO_PATTERN; 927 mumps->CleanUpMUMPS = PETSC_TRUE; 928 929 if (mumps->size > 1) { 930 PetscInt lsol_loc; 931 PetscScalar *sol_loc; 932 933 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 934 if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A; 935 else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A; 936 F_diag->assembled = PETSC_TRUE; 937 938 /* distributed solution; Create x_seq=sol_loc for repeated use */ 939 if (mumps->x_seq) { 940 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 941 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 942 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 943 } 944 lsol_loc = mumps->id.INFO(23); /* length of sol_loc */ 945 ierr = PetscMalloc2(lsol_loc,&sol_loc,lsol_loc,&mumps->id.isol_loc);CHKERRQ(ierr); 946 mumps->id.lsol_loc = lsol_loc; 947 #if defined(PETSC_USE_COMPLEX) 948 #if defined(PETSC_USE_REAL_SINGLE) 949 mumps->id.sol_loc = (mumps_complex*)sol_loc; 950 #else 951 mumps->id.sol_loc = (mumps_double_complex*)sol_loc; 952 #endif 953 #else 954 mumps->id.sol_loc = sol_loc; 955 #endif 956 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 957 } 958 PetscFunctionReturn(0); 959 } 960 961 /* Sets MUMPS options from the options database */ 962 #undef __FUNCT__ 963 #define __FUNCT__ "PetscSetMUMPSFromOptions" 964 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 965 { 966 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 967 PetscErrorCode ierr; 968 PetscInt icntl; 969 PetscBool flg; 970 971 PetscFunctionBegin; 972 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 973 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 974 if (flg) mumps->id.ICNTL(1) = icntl; 975 ierr = PetscOptionsInt("-mat_mumps_icntl_2","ICNTL(2): output stream for diagnostic printing, statistics, and warning","None",mumps->id.ICNTL(2),&icntl,&flg);CHKERRQ(ierr); 976 if (flg) mumps->id.ICNTL(2) = icntl; 977 ierr = PetscOptionsInt("-mat_mumps_icntl_3","ICNTL(3): output stream for global information, collected on the host","None",mumps->id.ICNTL(3),&icntl,&flg);CHKERRQ(ierr); 978 if (flg) mumps->id.ICNTL(3) = icntl; 979 980 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 981 if (flg) mumps->id.ICNTL(4) = icntl; 982 if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 983 984 ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): permutes to a zero-free diagonal and/or scale the matrix (0 to 7)","None",mumps->id.ICNTL(6),&icntl,&flg);CHKERRQ(ierr); 985 if (flg) mumps->id.ICNTL(6) = icntl; 986 987 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): computes a symmetric permutation in sequential analysis (0 to 7). 3=Scotch, 4=PORD, 5=Metis","None",mumps->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 988 if (flg) { 989 if (icntl== 1 && mumps->size > 1) 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"); 990 else mumps->id.ICNTL(7) = icntl; 991 } 992 993 ierr = PetscOptionsInt("-mat_mumps_icntl_8","ICNTL(8): scaling strategy (-2 to 8 or 77)","None",mumps->id.ICNTL(8),&mumps->id.ICNTL(8),NULL);CHKERRQ(ierr); 994 /* ierr = PetscOptionsInt("-mat_mumps_icntl_9","ICNTL(9): computes the solution using A or A^T","None",mumps->id.ICNTL(9),&mumps->id.ICNTL(9),NULL);CHKERRQ(ierr); handled by MatSolveTranspose_MUMPS() */ 995 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);CHKERRQ(ierr); 996 ierr = PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): statistics related to an error analysis (via -ksp_view)","None",mumps->id.ICNTL(11),&mumps->id.ICNTL(11),NULL);CHKERRQ(ierr); 997 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): an ordering strategy for symmetric matrices (0 to 3)","None",mumps->id.ICNTL(12),&mumps->id.ICNTL(12),NULL);CHKERRQ(ierr); 998 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): parallelism of the root node (enable ScaLAPACK) and its splitting","None",mumps->id.ICNTL(13),&mumps->id.ICNTL(13),NULL);CHKERRQ(ierr); 999 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): percentage increase in the estimated working space","None",mumps->id.ICNTL(14),&mumps->id.ICNTL(14),NULL);CHKERRQ(ierr); 1000 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): computes the Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);CHKERRQ(ierr); 1001 /* ierr = PetscOptionsInt("-mat_mumps_icntl_20","ICNTL(20): the format (dense or sparse) of the right-hand sides","None",mumps->id.ICNTL(20),&mumps->id.ICNTL(20),NULL);CHKERRQ(ierr); -- sparse rhs is not supported in PETSc API */ 1002 /* ierr = PetscOptionsInt("-mat_mumps_icntl_21","ICNTL(21): the distribution (centralized or distributed) of the solution vectors","None",mumps->id.ICNTL(21),&mumps->id.ICNTL(21),NULL);CHKERRQ(ierr); we only use distributed solution vector */ 1003 1004 ierr = PetscOptionsInt("-mat_mumps_icntl_22","ICNTL(22): in-core/out-of-core factorization and solve (0 or 1)","None",mumps->id.ICNTL(22),&mumps->id.ICNTL(22),NULL);CHKERRQ(ierr); 1005 ierr = PetscOptionsInt("-mat_mumps_icntl_23","ICNTL(23): max size of the working memory (MB) that can allocate per processor","None",mumps->id.ICNTL(23),&mumps->id.ICNTL(23),NULL);CHKERRQ(ierr); 1006 ierr = PetscOptionsInt("-mat_mumps_icntl_24","ICNTL(24): detection of null pivot rows (0 or 1)","None",mumps->id.ICNTL(24),&mumps->id.ICNTL(24),NULL);CHKERRQ(ierr); 1007 if (mumps->id.ICNTL(24)) { 1008 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 1009 } 1010 1011 ierr = PetscOptionsInt("-mat_mumps_icntl_25","ICNTL(25): compute a solution of a deficient matrix and a null space basis","None",mumps->id.ICNTL(25),&mumps->id.ICNTL(25),NULL);CHKERRQ(ierr); 1012 ierr = PetscOptionsInt("-mat_mumps_icntl_26","ICNTL(26): drives the solution phase if a Schur complement matrix","None",mumps->id.ICNTL(26),&mumps->id.ICNTL(26),NULL);CHKERRQ(ierr); 1013 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): the blocking size for multiple right-hand sides","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);CHKERRQ(ierr); 1014 ierr = PetscOptionsInt("-mat_mumps_icntl_28","ICNTL(28): use 1 for sequential analysis and ictnl(7) ordering, or 2 for parallel analysis and ictnl(29) ordering","None",mumps->id.ICNTL(28),&mumps->id.ICNTL(28),NULL);CHKERRQ(ierr); 1015 ierr = PetscOptionsInt("-mat_mumps_icntl_29","ICNTL(29): parallel ordering 1 = ptscotch, 2 = parmetis","None",mumps->id.ICNTL(29),&mumps->id.ICNTL(29),NULL);CHKERRQ(ierr); 1016 ierr = PetscOptionsInt("-mat_mumps_icntl_30","ICNTL(30): compute user-specified set of entries in inv(A)","None",mumps->id.ICNTL(30),&mumps->id.ICNTL(30),NULL);CHKERRQ(ierr); 1017 ierr = PetscOptionsInt("-mat_mumps_icntl_31","ICNTL(31): indicates which factors may be discarded during factorization","None",mumps->id.ICNTL(31),&mumps->id.ICNTL(31),NULL);CHKERRQ(ierr); 1018 /* ierr = PetscOptionsInt("-mat_mumps_icntl_32","ICNTL(32): performs the forward elemination of the right-hand sides during factorization","None",mumps->id.ICNTL(32),&mumps->id.ICNTL(32),NULL);CHKERRQ(ierr); -- not supported by PETSc API */ 1019 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);CHKERRQ(ierr); 1020 1021 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);CHKERRQ(ierr); 1022 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);CHKERRQ(ierr); 1023 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);CHKERRQ(ierr); 1024 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);CHKERRQ(ierr); 1025 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);CHKERRQ(ierr); 1026 1027 ierr = PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL); 1028 PetscOptionsEnd(); 1029 PetscFunctionReturn(0); 1030 } 1031 1032 #undef __FUNCT__ 1033 #define __FUNCT__ "PetscInitializeMUMPS" 1034 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps) 1035 { 1036 PetscErrorCode ierr; 1037 1038 PetscFunctionBegin; 1039 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid); 1040 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);CHKERRQ(ierr); 1041 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));CHKERRQ(ierr); 1042 1043 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 1044 1045 mumps->id.job = JOB_INIT; 1046 mumps->id.par = 1; /* host participates factorizaton and solve */ 1047 mumps->id.sym = mumps->sym; 1048 PetscMUMPS_c(&mumps->id); 1049 1050 mumps->CleanUpMUMPS = PETSC_FALSE; 1051 mumps->scat_rhs = NULL; 1052 mumps->scat_sol = NULL; 1053 1054 /* set PETSc-MUMPS default options - override MUMPS default */ 1055 mumps->id.ICNTL(3) = 0; 1056 mumps->id.ICNTL(4) = 0; 1057 if (mumps->size == 1) { 1058 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 1059 } else { 1060 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 1061 mumps->id.ICNTL(20) = 0; /* rhs is in dense format */ 1062 mumps->id.ICNTL(21) = 1; /* distributed solution */ 1063 } 1064 PetscFunctionReturn(0); 1065 } 1066 1067 /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */ 1068 #undef __FUNCT__ 1069 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 1070 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1071 { 1072 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1073 PetscErrorCode ierr; 1074 Vec b; 1075 IS is_iden; 1076 const PetscInt M = A->rmap->N; 1077 1078 PetscFunctionBegin; 1079 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1080 1081 /* Set MUMPS options from the options database */ 1082 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1083 1084 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1085 1086 /* analysis phase */ 1087 /*----------------*/ 1088 mumps->id.job = JOB_FACTSYMBOLIC; 1089 mumps->id.n = M; 1090 switch (mumps->id.ICNTL(18)) { 1091 case 0: /* centralized assembled matrix input */ 1092 if (!mumps->myid) { 1093 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1094 if (mumps->id.ICNTL(6)>1) { 1095 #if defined(PETSC_USE_COMPLEX) 1096 #if defined(PETSC_USE_REAL_SINGLE) 1097 mumps->id.a = (mumps_complex*)mumps->val; 1098 #else 1099 mumps->id.a = (mumps_double_complex*)mumps->val; 1100 #endif 1101 #else 1102 mumps->id.a = mumps->val; 1103 #endif 1104 } 1105 if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */ 1106 /* 1107 PetscBool flag; 1108 ierr = ISEqual(r,c,&flag);CHKERRQ(ierr); 1109 if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm"); 1110 ierr = ISView(r,PETSC_VIEWER_STDOUT_SELF); 1111 */ 1112 if (!mumps->myid) { 1113 const PetscInt *idx; 1114 PetscInt i,*perm_in; 1115 1116 ierr = PetscMalloc1(M,&perm_in);CHKERRQ(ierr); 1117 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 1118 1119 mumps->id.perm_in = perm_in; 1120 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 1121 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 1122 } 1123 } 1124 } 1125 break; 1126 case 3: /* distributed assembled matrix input (size>1) */ 1127 mumps->id.nz_loc = mumps->nz; 1128 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1129 if (mumps->id.ICNTL(6)>1) { 1130 #if defined(PETSC_USE_COMPLEX) 1131 #if defined(PETSC_USE_REAL_SINGLE) 1132 mumps->id.a_loc = (mumps_complex*)mumps->val; 1133 #else 1134 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1135 #endif 1136 #else 1137 mumps->id.a_loc = mumps->val; 1138 #endif 1139 } 1140 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1141 if (!mumps->myid) { 1142 ierr = VecCreateSeq(PETSC_COMM_SELF,A->rmap->N,&mumps->b_seq);CHKERRQ(ierr); 1143 ierr = ISCreateStride(PETSC_COMM_SELF,A->rmap->N,0,1,&is_iden);CHKERRQ(ierr); 1144 } else { 1145 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1146 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1147 } 1148 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1149 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1150 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1151 ierr = VecDestroy(&b);CHKERRQ(ierr); 1152 break; 1153 } 1154 PetscMUMPS_c(&mumps->id); 1155 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 1156 1157 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1158 F->ops->solve = MatSolve_MUMPS; 1159 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1160 F->ops->matsolve = MatMatSolve_MUMPS; 1161 PetscFunctionReturn(0); 1162 } 1163 1164 /* Note the Petsc r and c permutations are ignored */ 1165 #undef __FUNCT__ 1166 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 1167 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1168 { 1169 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1170 PetscErrorCode ierr; 1171 Vec b; 1172 IS is_iden; 1173 const PetscInt M = A->rmap->N; 1174 1175 PetscFunctionBegin; 1176 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1177 1178 /* Set MUMPS options from the options database */ 1179 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1180 1181 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1182 1183 /* analysis phase */ 1184 /*----------------*/ 1185 mumps->id.job = JOB_FACTSYMBOLIC; 1186 mumps->id.n = M; 1187 switch (mumps->id.ICNTL(18)) { 1188 case 0: /* centralized assembled matrix input */ 1189 if (!mumps->myid) { 1190 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1191 if (mumps->id.ICNTL(6)>1) { 1192 #if defined(PETSC_USE_COMPLEX) 1193 #if defined(PETSC_USE_REAL_SINGLE) 1194 mumps->id.a = (mumps_complex*)mumps->val; 1195 #else 1196 mumps->id.a = (mumps_double_complex*)mumps->val; 1197 #endif 1198 #else 1199 mumps->id.a = mumps->val; 1200 #endif 1201 } 1202 } 1203 break; 1204 case 3: /* distributed assembled matrix input (size>1) */ 1205 mumps->id.nz_loc = mumps->nz; 1206 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1207 if (mumps->id.ICNTL(6)>1) { 1208 #if defined(PETSC_USE_COMPLEX) 1209 #if defined(PETSC_USE_REAL_SINGLE) 1210 mumps->id.a_loc = (mumps_complex*)mumps->val; 1211 #else 1212 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1213 #endif 1214 #else 1215 mumps->id.a_loc = mumps->val; 1216 #endif 1217 } 1218 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1219 if (!mumps->myid) { 1220 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1221 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1222 } else { 1223 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1224 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1225 } 1226 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1227 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1228 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1229 ierr = VecDestroy(&b);CHKERRQ(ierr); 1230 break; 1231 } 1232 PetscMUMPS_c(&mumps->id); 1233 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 1234 1235 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1236 F->ops->solve = MatSolve_MUMPS; 1237 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1238 PetscFunctionReturn(0); 1239 } 1240 1241 /* Note the Petsc r permutation and factor info are ignored */ 1242 #undef __FUNCT__ 1243 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 1244 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 1245 { 1246 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1247 PetscErrorCode ierr; 1248 Vec b; 1249 IS is_iden; 1250 const PetscInt M = A->rmap->N; 1251 1252 PetscFunctionBegin; 1253 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1254 1255 /* Set MUMPS options from the options database */ 1256 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1257 1258 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1259 1260 /* analysis phase */ 1261 /*----------------*/ 1262 mumps->id.job = JOB_FACTSYMBOLIC; 1263 mumps->id.n = M; 1264 switch (mumps->id.ICNTL(18)) { 1265 case 0: /* centralized assembled matrix input */ 1266 if (!mumps->myid) { 1267 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1268 if (mumps->id.ICNTL(6)>1) { 1269 #if defined(PETSC_USE_COMPLEX) 1270 #if defined(PETSC_USE_REAL_SINGLE) 1271 mumps->id.a = (mumps_complex*)mumps->val; 1272 #else 1273 mumps->id.a = (mumps_double_complex*)mumps->val; 1274 #endif 1275 #else 1276 mumps->id.a = mumps->val; 1277 #endif 1278 } 1279 } 1280 break; 1281 case 3: /* distributed assembled matrix input (size>1) */ 1282 mumps->id.nz_loc = mumps->nz; 1283 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1284 if (mumps->id.ICNTL(6)>1) { 1285 #if defined(PETSC_USE_COMPLEX) 1286 #if defined(PETSC_USE_REAL_SINGLE) 1287 mumps->id.a_loc = (mumps_complex*)mumps->val; 1288 #else 1289 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1290 #endif 1291 #else 1292 mumps->id.a_loc = mumps->val; 1293 #endif 1294 } 1295 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1296 if (!mumps->myid) { 1297 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1298 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1299 } else { 1300 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1301 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1302 } 1303 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1304 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1305 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1306 ierr = VecDestroy(&b);CHKERRQ(ierr); 1307 break; 1308 } 1309 PetscMUMPS_c(&mumps->id); 1310 if (mumps->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",mumps->id.INFOG(1)); 1311 1312 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1313 F->ops->solve = MatSolve_MUMPS; 1314 F->ops->solvetranspose = MatSolve_MUMPS; 1315 F->ops->matsolve = MatMatSolve_MUMPS; 1316 #if defined(PETSC_USE_COMPLEX) 1317 F->ops->getinertia = NULL; 1318 #else 1319 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1320 #endif 1321 PetscFunctionReturn(0); 1322 } 1323 1324 //update!!! 1325 #undef __FUNCT__ 1326 #define __FUNCT__ "MatView_MUMPS" 1327 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1328 { 1329 PetscErrorCode ierr; 1330 PetscBool iascii; 1331 PetscViewerFormat format; 1332 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1333 1334 PetscFunctionBegin; 1335 /* check if matrix is mumps type */ 1336 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1337 1338 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1339 if (iascii) { 1340 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1341 if (format == PETSC_VIEWER_ASCII_INFO) { 1342 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1343 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1344 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1345 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));CHKERRQ(ierr); 1346 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));CHKERRQ(ierr); 1347 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));CHKERRQ(ierr); 1348 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));CHKERRQ(ierr); 1349 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));CHKERRQ(ierr); 1350 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));CHKERRQ(ierr); 1351 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));CHKERRQ(ierr); 1352 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));CHKERRQ(ierr); 1353 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));CHKERRQ(ierr); 1354 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));CHKERRQ(ierr); 1355 if (mumps->id.ICNTL(11)>0) { 1356 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));CHKERRQ(ierr); 1357 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));CHKERRQ(ierr); 1358 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));CHKERRQ(ierr); 1359 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));CHKERRQ(ierr); 1360 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));CHKERRQ(ierr); 1361 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));CHKERRQ(ierr); 1362 } 1363 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));CHKERRQ(ierr); 1364 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));CHKERRQ(ierr); 1365 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));CHKERRQ(ierr); 1366 /* ICNTL(15-17) not used */ 1367 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));CHKERRQ(ierr); 1368 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));CHKERRQ(ierr); 1369 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));CHKERRQ(ierr); 1370 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));CHKERRQ(ierr); 1371 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));CHKERRQ(ierr); 1372 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));CHKERRQ(ierr); 1373 1374 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));CHKERRQ(ierr); 1375 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));CHKERRQ(ierr); 1376 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));CHKERRQ(ierr); 1377 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));CHKERRQ(ierr); 1378 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));CHKERRQ(ierr); 1379 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));CHKERRQ(ierr); 1380 1381 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));CHKERRQ(ierr); 1382 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));CHKERRQ(ierr); 1383 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));CHKERRQ(ierr); 1384 1385 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));CHKERRQ(ierr); 1386 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));CHKERRQ(ierr); 1387 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));CHKERRQ(ierr); 1388 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));CHKERRQ(ierr); 1389 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));CHKERRQ(ierr); 1390 1391 /* infomation local to each processor */ 1392 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1393 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1394 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));CHKERRQ(ierr); 1395 ierr = PetscViewerFlush(viewer); 1396 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1397 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));CHKERRQ(ierr); 1398 ierr = PetscViewerFlush(viewer); 1399 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1400 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));CHKERRQ(ierr); 1401 ierr = PetscViewerFlush(viewer); 1402 1403 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1404 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));CHKERRQ(ierr); 1405 ierr = PetscViewerFlush(viewer); 1406 1407 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1408 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));CHKERRQ(ierr); 1409 ierr = PetscViewerFlush(viewer); 1410 1411 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1412 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));CHKERRQ(ierr); 1413 ierr = PetscViewerFlush(viewer); 1414 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1415 1416 if (!mumps->myid) { /* information from the host */ 1417 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));CHKERRQ(ierr); 1418 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));CHKERRQ(ierr); 1419 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));CHKERRQ(ierr); 1420 ierr = PetscViewerASCIIPrintf(viewer," (RINFOG(12) RINFOG(13))*2^INFOG(34) (determinant): (%g,%g)*(2^%d)\n",mumps->id.RINFOG(12),mumps->id.RINFOG(13),mumps->id.INFOG(34));CHKERRQ(ierr); 1421 1422 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));CHKERRQ(ierr); 1423 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));CHKERRQ(ierr); 1424 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));CHKERRQ(ierr); 1425 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));CHKERRQ(ierr); 1426 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));CHKERRQ(ierr); 1427 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));CHKERRQ(ierr); 1428 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));CHKERRQ(ierr); 1429 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));CHKERRQ(ierr); 1430 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));CHKERRQ(ierr); 1431 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));CHKERRQ(ierr); 1432 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));CHKERRQ(ierr); 1433 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));CHKERRQ(ierr); 1434 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));CHKERRQ(ierr); 1435 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",mumps->id.INFOG(16));CHKERRQ(ierr); 1436 ierr = PetscViewerASCIIPrintf(viewer," INFOG(17) (estimated size of all MUMPS internal data for factorization after analysis: sum over all processors): %d \n",mumps->id.INFOG(17));CHKERRQ(ierr); 1437 ierr = PetscViewerASCIIPrintf(viewer," INFOG(18) (size of all MUMPS internal data allocated during factorization: value on the most memory consuming processor): %d \n",mumps->id.INFOG(18));CHKERRQ(ierr); 1438 ierr = PetscViewerASCIIPrintf(viewer," INFOG(19) (size of all MUMPS internal data allocated during factorization: sum over all processors): %d \n",mumps->id.INFOG(19));CHKERRQ(ierr); 1439 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));CHKERRQ(ierr); 1440 ierr = PetscViewerASCIIPrintf(viewer," INFOG(21) (size in MB of memory effectively used during factorization - value on the most memory consuming processor): %d \n",mumps->id.INFOG(21));CHKERRQ(ierr); 1441 ierr = PetscViewerASCIIPrintf(viewer," INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): %d \n",mumps->id.INFOG(22));CHKERRQ(ierr); 1442 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));CHKERRQ(ierr); 1443 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));CHKERRQ(ierr); 1444 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));CHKERRQ(ierr); 1445 ierr = PetscViewerASCIIPrintf(viewer," INFOG(28) (after factorization: number of null pivots encountered): %d\n",mumps->id.INFOG(28));CHKERRQ(ierr); 1446 ierr = PetscViewerASCIIPrintf(viewer," INFOG(29) (after factorization: effective number of entries in the factors (sum over all processors)): %d\n",mumps->id.INFOG(29));CHKERRQ(ierr); 1447 ierr = PetscViewerASCIIPrintf(viewer," INFOG(30, 31) (after solution: size in Mbytes of memory used during solution phase): %d, %d\n",mumps->id.INFOG(30),mumps->id.INFOG(31));CHKERRQ(ierr); 1448 ierr = PetscViewerASCIIPrintf(viewer," INFOG(32) (after analysis: type of analysis done): %d\n",mumps->id.INFOG(32));CHKERRQ(ierr); 1449 ierr = PetscViewerASCIIPrintf(viewer," INFOG(33) (value used for ICNTL(8)): %d\n",mumps->id.INFOG(33));CHKERRQ(ierr); 1450 ierr = PetscViewerASCIIPrintf(viewer," INFOG(34) (exponent of the determinant if determinant is requested): %d\n",mumps->id.INFOG(34));CHKERRQ(ierr); 1451 } 1452 } 1453 } 1454 PetscFunctionReturn(0); 1455 } 1456 1457 #undef __FUNCT__ 1458 #define __FUNCT__ "MatGetInfo_MUMPS" 1459 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1460 { 1461 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1462 1463 PetscFunctionBegin; 1464 info->block_size = 1.0; 1465 info->nz_allocated = mumps->id.INFOG(20); 1466 info->nz_used = mumps->id.INFOG(20); 1467 info->nz_unneeded = 0.0; 1468 info->assemblies = 0.0; 1469 info->mallocs = 0.0; 1470 info->memory = 0.0; 1471 info->fill_ratio_given = 0; 1472 info->fill_ratio_needed = 0; 1473 info->factor_mallocs = 0; 1474 PetscFunctionReturn(0); 1475 } 1476 1477 /* -------------------------------------------------------------------------------------------*/ 1478 #undef __FUNCT__ 1479 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1480 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1481 { 1482 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1483 1484 PetscFunctionBegin; 1485 mumps->id.ICNTL(icntl) = ival; 1486 PetscFunctionReturn(0); 1487 } 1488 1489 #undef __FUNCT__ 1490 #define __FUNCT__ "MatMumpsGetIcntl_MUMPS" 1491 PetscErrorCode MatMumpsGetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt *ival) 1492 { 1493 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1494 1495 PetscFunctionBegin; 1496 *ival = mumps->id.ICNTL(icntl); 1497 PetscFunctionReturn(0); 1498 } 1499 1500 #undef __FUNCT__ 1501 #define __FUNCT__ "MatMumpsSetIcntl" 1502 /*@ 1503 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1504 1505 Logically Collective on Mat 1506 1507 Input Parameters: 1508 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1509 . icntl - index of MUMPS parameter array ICNTL() 1510 - ival - value of MUMPS ICNTL(icntl) 1511 1512 Options Database: 1513 . -mat_mumps_icntl_<icntl> <ival> 1514 1515 Level: beginner 1516 1517 References: MUMPS Users' Guide 1518 1519 .seealso: MatGetFactor() 1520 @*/ 1521 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1522 { 1523 PetscErrorCode ierr; 1524 1525 PetscFunctionBegin; 1526 PetscValidLogicalCollectiveInt(F,icntl,2); 1527 PetscValidLogicalCollectiveInt(F,ival,3); 1528 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1529 PetscFunctionReturn(0); 1530 } 1531 1532 #undef __FUNCT__ 1533 #define __FUNCT__ "MatMumpsGetIcntl" 1534 /*@ 1535 MatMumpsGetIcntl - Get MUMPS parameter ICNTL() 1536 1537 Logically Collective on Mat 1538 1539 Input Parameters: 1540 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1541 - icntl - index of MUMPS parameter array ICNTL() 1542 1543 Output Parameter: 1544 . ival - value of MUMPS ICNTL(icntl) 1545 1546 Level: beginner 1547 1548 References: MUMPS Users' Guide 1549 1550 .seealso: MatGetFactor() 1551 @*/ 1552 PetscErrorCode MatMumpsGetIcntl(Mat F,PetscInt icntl,PetscInt *ival) 1553 { 1554 PetscErrorCode ierr; 1555 1556 PetscFunctionBegin; 1557 PetscValidLogicalCollectiveInt(F,icntl,2); 1558 PetscValidIntPointer(ival,3); 1559 ierr = PetscTryMethod(F,"MatMumpsGetIcntl_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1560 PetscFunctionReturn(0); 1561 } 1562 1563 /* -------------------------------------------------------------------------------------------*/ 1564 #undef __FUNCT__ 1565 #define __FUNCT__ "MatMumpsSetCntl_MUMPS" 1566 PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val) 1567 { 1568 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1569 1570 PetscFunctionBegin; 1571 mumps->id.CNTL(icntl) = val; 1572 PetscFunctionReturn(0); 1573 } 1574 1575 #undef __FUNCT__ 1576 #define __FUNCT__ "MatMumpsGetCntl_MUMPS" 1577 PetscErrorCode MatMumpsGetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal *val) 1578 { 1579 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1580 1581 PetscFunctionBegin; 1582 *val = mumps->id.CNTL(icntl); 1583 PetscFunctionReturn(0); 1584 } 1585 1586 #undef __FUNCT__ 1587 #define __FUNCT__ "MatMumpsSetCntl" 1588 /*@ 1589 MatMumpsSetCntl - Set MUMPS parameter CNTL() 1590 1591 Logically Collective on Mat 1592 1593 Input Parameters: 1594 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1595 . icntl - index of MUMPS parameter array CNTL() 1596 - val - value of MUMPS CNTL(icntl) 1597 1598 Options Database: 1599 . -mat_mumps_cntl_<icntl> <val> 1600 1601 Level: beginner 1602 1603 References: MUMPS Users' Guide 1604 1605 .seealso: MatGetFactor() 1606 @*/ 1607 PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val) 1608 { 1609 PetscErrorCode ierr; 1610 1611 PetscFunctionBegin; 1612 PetscValidLogicalCollectiveInt(F,icntl,2); 1613 PetscValidLogicalCollectiveReal(F,val,3); 1614 ierr = PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));CHKERRQ(ierr); 1615 PetscFunctionReturn(0); 1616 } 1617 1618 #undef __FUNCT__ 1619 #define __FUNCT__ "MatMumpsGetCntl" 1620 /*@ 1621 MatMumpsGetCntl - Get MUMPS parameter CNTL() 1622 1623 Logically Collective on Mat 1624 1625 Input Parameters: 1626 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1627 - icntl - index of MUMPS parameter array CNTL() 1628 1629 Output Parameter: 1630 . val - value of MUMPS CNTL(icntl) 1631 1632 Level: beginner 1633 1634 References: MUMPS Users' Guide 1635 1636 .seealso: MatGetFactor() 1637 @*/ 1638 PetscErrorCode MatMumpsGetCntl(Mat F,PetscInt icntl,PetscReal *val) 1639 { 1640 PetscErrorCode ierr; 1641 1642 PetscFunctionBegin; 1643 PetscValidLogicalCollectiveInt(F,icntl,2); 1644 PetscValidRealPointer(val,3); 1645 ierr = PetscTryMethod(F,"MatMumpsGetCntl_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1646 PetscFunctionReturn(0); 1647 } 1648 1649 #undef __FUNCT__ 1650 #define __FUNCT__ "MatMumpsGetInfo_MUMPS" 1651 PetscErrorCode MatMumpsGetInfo_MUMPS(Mat F,PetscInt icntl,PetscInt *info) 1652 { 1653 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1654 1655 PetscFunctionBegin; 1656 *info = mumps->id.INFO(icntl); 1657 PetscFunctionReturn(0); 1658 } 1659 1660 #undef __FUNCT__ 1661 #define __FUNCT__ "MatMumpsGetInfog_MUMPS" 1662 PetscErrorCode MatMumpsGetInfog_MUMPS(Mat F,PetscInt icntl,PetscInt *infog) 1663 { 1664 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1665 1666 PetscFunctionBegin; 1667 *infog = mumps->id.INFOG(icntl); 1668 PetscFunctionReturn(0); 1669 } 1670 1671 #undef __FUNCT__ 1672 #define __FUNCT__ "MatMumpsGetRinfo_MUMPS" 1673 PetscErrorCode MatMumpsGetRinfo_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfo) 1674 { 1675 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1676 1677 PetscFunctionBegin; 1678 *rinfo = mumps->id.RINFO(icntl); 1679 PetscFunctionReturn(0); 1680 } 1681 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "MatMumpsGetRinfog_MUMPS" 1684 PetscErrorCode MatMumpsGetRinfog_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfog) 1685 { 1686 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1687 1688 PetscFunctionBegin; 1689 *rinfog = mumps->id.RINFOG(icntl); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 #undef __FUNCT__ 1694 #define __FUNCT__ "MatMumpsGetInfo" 1695 /*@ 1696 MatMumpsGetInfo - Get MUMPS parameter INFO() 1697 1698 Logically Collective on Mat 1699 1700 Input Parameters: 1701 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1702 - icntl - index of MUMPS parameter array INFO() 1703 1704 Output Parameter: 1705 . ival - value of MUMPS INFO(icntl) 1706 1707 Level: beginner 1708 1709 References: MUMPS Users' Guide 1710 1711 .seealso: MatGetFactor() 1712 @*/ 1713 PetscErrorCode MatMumpsGetInfo(Mat F,PetscInt icntl,PetscInt *ival) 1714 { 1715 PetscErrorCode ierr; 1716 1717 PetscFunctionBegin; 1718 PetscValidIntPointer(ival,3); 1719 ierr = PetscTryMethod(F,"MatMumpsGetInfo_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1720 PetscFunctionReturn(0); 1721 } 1722 1723 #undef __FUNCT__ 1724 #define __FUNCT__ "MatMumpsGetInfog" 1725 /*@ 1726 MatMumpsGetInfog - Get MUMPS parameter INFOG() 1727 1728 Logically Collective on Mat 1729 1730 Input Parameters: 1731 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1732 - icntl - index of MUMPS parameter array INFOG() 1733 1734 Output Parameter: 1735 . ival - value of MUMPS INFOG(icntl) 1736 1737 Level: beginner 1738 1739 References: MUMPS Users' Guide 1740 1741 .seealso: MatGetFactor() 1742 @*/ 1743 PetscErrorCode MatMumpsGetInfog(Mat F,PetscInt icntl,PetscInt *ival) 1744 { 1745 PetscErrorCode ierr; 1746 1747 PetscFunctionBegin; 1748 PetscValidIntPointer(ival,3); 1749 ierr = PetscTryMethod(F,"MatMumpsGetInfog_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "MatMumpsGetRinfo" 1755 /*@ 1756 MatMumpsGetRinfo - Get MUMPS parameter RINFO() 1757 1758 Logically Collective on Mat 1759 1760 Input Parameters: 1761 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1762 - icntl - index of MUMPS parameter array RINFO() 1763 1764 Output Parameter: 1765 . val - value of MUMPS RINFO(icntl) 1766 1767 Level: beginner 1768 1769 References: MUMPS Users' Guide 1770 1771 .seealso: MatGetFactor() 1772 @*/ 1773 PetscErrorCode MatMumpsGetRinfo(Mat F,PetscInt icntl,PetscReal *val) 1774 { 1775 PetscErrorCode ierr; 1776 1777 PetscFunctionBegin; 1778 PetscValidRealPointer(val,3); 1779 ierr = PetscTryMethod(F,"MatMumpsGetRinfo_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1780 PetscFunctionReturn(0); 1781 } 1782 1783 #undef __FUNCT__ 1784 #define __FUNCT__ "MatMumpsGetRinfog" 1785 /*@ 1786 MatMumpsGetRinfog - Get MUMPS parameter RINFOG() 1787 1788 Logically Collective on Mat 1789 1790 Input Parameters: 1791 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1792 - icntl - index of MUMPS parameter array RINFOG() 1793 1794 Output Parameter: 1795 . val - value of MUMPS RINFOG(icntl) 1796 1797 Level: beginner 1798 1799 References: MUMPS Users' Guide 1800 1801 .seealso: MatGetFactor() 1802 @*/ 1803 PetscErrorCode MatMumpsGetRinfog(Mat F,PetscInt icntl,PetscReal *val) 1804 { 1805 PetscErrorCode ierr; 1806 1807 PetscFunctionBegin; 1808 PetscValidRealPointer(val,3); 1809 ierr = PetscTryMethod(F,"MatMumpsGetRinfog_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1810 PetscFunctionReturn(0); 1811 } 1812 1813 /*MC 1814 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1815 distributed and sequential matrices via the external package MUMPS. 1816 1817 Works with MATAIJ and MATSBAIJ matrices 1818 1819 Options Database Keys: 1820 + -mat_mumps_icntl_1 <6>: ICNTL(1): output stream for error messages (None) 1821 . -mat_mumps_icntl_2 <0>: ICNTL(2): output stream for diagnostic printing, statistics, and warning (None) 1822 . -mat_mumps_icntl_3 <0>: ICNTL(3): output stream for global information, collected on the host (None) 1823 . -mat_mumps_icntl_4 <0>: ICNTL(4): level of printing (0 to 4) (None) 1824 . -mat_mumps_icntl_6 <7>: ICNTL(6): permutes to a zero-free diagonal and/or scale the matrix (0 to 7) (None) 1825 . -mat_mumps_icntl_7 <7>: ICNTL(7): computes a symmetric permutation in sequential analysis (0 to 7). 3=Scotch, 4=PORD, 5=Metis (None) 1826 . -mat_mumps_icntl_8 <77>: ICNTL(8): scaling strategy (-2 to 8 or 77) (None) 1827 . -mat_mumps_icntl_10 <0>: ICNTL(10): max num of refinements (None) 1828 . -mat_mumps_icntl_11 <0>: ICNTL(11): statistics related to an error analysis (via -ksp_view) (None) 1829 . -mat_mumps_icntl_12 <1>: ICNTL(12): an ordering strategy for symmetric matrices (0 to 3) (None) 1830 . -mat_mumps_icntl_13 <0>: ICNTL(13): parallelism of the root node (enable ScaLAPACK) and its splitting (None) 1831 . -mat_mumps_icntl_14 <20>: ICNTL(14): percentage increase in the estimated working space (None) 1832 . -mat_mumps_icntl_19 <0>: ICNTL(19): computes the Schur complement (None) 1833 . -mat_mumps_icntl_22 <0>: ICNTL(22): in-core/out-of-core factorization and solve (0 or 1) (None) 1834 . -mat_mumps_icntl_23 <0>: ICNTL(23): max size of the working memory (MB) that can allocate per processor (None) 1835 . -mat_mumps_icntl_24 <0>: ICNTL(24): detection of null pivot rows (0 or 1) (None) 1836 . -mat_mumps_icntl_25 <0>: ICNTL(25): compute a solution of a deficient matrix and a null space basis (None) 1837 . -mat_mumps_icntl_26 <0>: ICNTL(26): drives the solution phase if a Schur complement matrix (None) 1838 . -mat_mumps_icntl_28 <1>: ICNTL(28): use 1 for sequential analysis and ictnl(7) ordering, or 2 for parallel analysis and ictnl(29) ordering (None) 1839 . -mat_mumps_icntl_29 <0>: ICNTL(29): parallel ordering 1 = ptscotch, 2 = parmetis (None) 1840 . -mat_mumps_icntl_30 <0>: ICNTL(30): compute user-specified set of entries in inv(A) (None) 1841 . -mat_mumps_icntl_31 <0>: ICNTL(31): indicates which factors may be discarded during factorization (None) 1842 . -mat_mumps_icntl_33 <0>: ICNTL(33): compute determinant (None) 1843 . -mat_mumps_cntl_1 <0.01>: CNTL(1): relative pivoting threshold (None) 1844 . -mat_mumps_cntl_2 <1.49012e-08>: CNTL(2): stopping criterion of refinement (None) 1845 . -mat_mumps_cntl_3 <0>: CNTL(3): absolute pivoting threshold (None) 1846 . -mat_mumps_cntl_4 <-1>: CNTL(4): value for static pivoting (None) 1847 - -mat_mumps_cntl_5 <0>: CNTL(5): fixation for null pivots (None) 1848 1849 Level: beginner 1850 1851 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1852 1853 M*/ 1854 1855 #undef __FUNCT__ 1856 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1857 static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1858 { 1859 PetscFunctionBegin; 1860 *type = MATSOLVERMUMPS; 1861 PetscFunctionReturn(0); 1862 } 1863 1864 /* MatGetFactor for Seq and MPI AIJ matrices */ 1865 #undef __FUNCT__ 1866 #define __FUNCT__ "MatGetFactor_aij_mumps" 1867 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1868 { 1869 Mat B; 1870 PetscErrorCode ierr; 1871 Mat_MUMPS *mumps; 1872 PetscBool isSeqAIJ; 1873 1874 PetscFunctionBegin; 1875 /* Create the factorization matrix */ 1876 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1877 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1878 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1879 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1880 if (isSeqAIJ) { 1881 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1882 } else { 1883 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1884 } 1885 1886 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1887 1888 B->ops->view = MatView_MUMPS; 1889 B->ops->getinfo = MatGetInfo_MUMPS; 1890 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1891 1892 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1893 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1894 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1895 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1896 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1897 1898 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1899 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1900 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1901 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1902 if (ftype == MAT_FACTOR_LU) { 1903 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1904 B->factortype = MAT_FACTOR_LU; 1905 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1906 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1907 mumps->sym = 0; 1908 } else { 1909 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1910 B->factortype = MAT_FACTOR_CHOLESKY; 1911 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1912 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1913 if (A->spd_set && A->spd) mumps->sym = 1; 1914 else mumps->sym = 2; 1915 } 1916 1917 mumps->isAIJ = PETSC_TRUE; 1918 mumps->Destroy = B->ops->destroy; 1919 B->ops->destroy = MatDestroy_MUMPS; 1920 B->spptr = (void*)mumps; 1921 1922 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1923 1924 *F = B; 1925 PetscFunctionReturn(0); 1926 } 1927 1928 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1929 #undef __FUNCT__ 1930 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1931 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1932 { 1933 Mat B; 1934 PetscErrorCode ierr; 1935 Mat_MUMPS *mumps; 1936 PetscBool isSeqSBAIJ; 1937 1938 PetscFunctionBegin; 1939 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1940 if (A->rmap->bs > 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with block size > 1 with MUMPS Cholesky, use AIJ matrix instead"); 1941 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1942 /* Create the factorization matrix */ 1943 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1944 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1945 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1946 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1947 if (isSeqSBAIJ) { 1948 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 1949 1950 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1951 } else { 1952 ierr = MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);CHKERRQ(ierr); 1953 1954 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1955 } 1956 1957 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1958 B->ops->view = MatView_MUMPS; 1959 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1960 1961 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1962 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1963 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1964 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1965 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1966 1967 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1968 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1969 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1970 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1971 1972 B->factortype = MAT_FACTOR_CHOLESKY; 1973 if (A->spd_set && A->spd) mumps->sym = 1; 1974 else mumps->sym = 2; 1975 1976 mumps->isAIJ = PETSC_FALSE; 1977 mumps->Destroy = B->ops->destroy; 1978 B->ops->destroy = MatDestroy_MUMPS; 1979 B->spptr = (void*)mumps; 1980 1981 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1982 1983 *F = B; 1984 PetscFunctionReturn(0); 1985 } 1986 1987 #undef __FUNCT__ 1988 #define __FUNCT__ "MatGetFactor_baij_mumps" 1989 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 1990 { 1991 Mat B; 1992 PetscErrorCode ierr; 1993 Mat_MUMPS *mumps; 1994 PetscBool isSeqBAIJ; 1995 1996 PetscFunctionBegin; 1997 /* Create the factorization matrix */ 1998 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1999 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 2000 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 2001 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2002 if (isSeqBAIJ) { 2003 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);CHKERRQ(ierr); 2004 } else { 2005 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 2006 } 2007 2008 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 2009 if (ftype == MAT_FACTOR_LU) { 2010 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 2011 B->factortype = MAT_FACTOR_LU; 2012 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 2013 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 2014 mumps->sym = 0; 2015 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 2016 2017 B->ops->view = MatView_MUMPS; 2018 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 2019 2020 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 2021 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 2022 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 2023 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 2024 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 2025 2026 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 2027 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 2028 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 2029 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 2030 2031 mumps->isAIJ = PETSC_TRUE; 2032 mumps->Destroy = B->ops->destroy; 2033 B->ops->destroy = MatDestroy_MUMPS; 2034 B->spptr = (void*)mumps; 2035 2036 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 2037 2038 *F = B; 2039 PetscFunctionReturn(0); 2040 } 2041 2042 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 2043 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat,MatFactorType,Mat*); 2044 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 2045 2046 #undef __FUNCT__ 2047 #define __FUNCT__ "MatSolverPackageRegister_MUMPS" 2048 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MUMPS(void) 2049 { 2050 PetscErrorCode ierr; 2051 2052 PetscFunctionBegin; 2053 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2054 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2055 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2056 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2057 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPISBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2058 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2059 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2060 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2061 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2062 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2063 PetscFunctionReturn(0); 2064 } 2065 2066