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