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