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; 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,NULL,&nrhs);CHKERRQ(ierr); 691 mumps->id.nrhs = nrhs; 692 ierr = MatCopy(B,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 693 ierr = MatDenseGetArray(X,&array);CHKERRQ(ierr); 694 #if defined(PETSC_USE_COMPLEX) 695 #if defined(PETSC_USE_REAL_SINGLE) 696 mumps->id.rhs = (mumps_complex*)array; 697 #else 698 mumps->id.rhs = (mumps_double_complex*)array; 699 #endif 700 #else 701 mumps->id.rhs = array; 702 #endif 703 ierr = MatDenseRestoreArray(X,&array);CHKERRQ(ierr); 704 705 /* solve phase */ 706 /*-------------*/ 707 mumps->id.job = JOB_SOLVE; 708 PetscMUMPS_c(&mumps->id); 709 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)); 710 PetscFunctionReturn(0); 711 } 712 713 #if !defined(PETSC_USE_COMPLEX) 714 /* 715 input: 716 F: numeric factor 717 output: 718 nneg: total number of negative pivots 719 nzero: 0 720 npos: (global dimension of F) - nneg 721 */ 722 723 #undef __FUNCT__ 724 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 725 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 726 { 727 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 728 PetscErrorCode ierr; 729 PetscMPIInt size; 730 731 PetscFunctionBegin; 732 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)F),&size);CHKERRQ(ierr); 733 /* 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 */ 734 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)); 735 736 if (nneg) *nneg = mumps->id.INFOG(12); 737 if (nzero || npos) { 738 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"); 739 if (nzero) *nzero = mumps->id.INFOG(28); 740 if (npos) *npos = F->rmap->N - (mumps->id.INFOG(12) + mumps->id.INFOG(28)); 741 } 742 PetscFunctionReturn(0); 743 } 744 #endif /* !defined(PETSC_USE_COMPLEX) */ 745 746 #undef __FUNCT__ 747 #define __FUNCT__ "MatFactorNumeric_MUMPS" 748 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 749 { 750 Mat_MUMPS *mumps =(Mat_MUMPS*)(F)->spptr; 751 PetscErrorCode ierr; 752 Mat F_diag; 753 PetscBool isMPIAIJ; 754 755 PetscFunctionBegin; 756 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_REUSE_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 757 758 /* numerical factorization phase */ 759 /*-------------------------------*/ 760 mumps->id.job = JOB_FACTNUMERIC; 761 if (!mumps->id.ICNTL(18)) { 762 if (!mumps->myid) { 763 #if defined(PETSC_USE_COMPLEX) 764 #if defined(PETSC_USE_REAL_SINGLE) 765 mumps->id.a = (mumps_complex*)mumps->val; 766 #else 767 mumps->id.a = (mumps_double_complex*)mumps->val; 768 #endif 769 #else 770 mumps->id.a = mumps->val; 771 #endif 772 } 773 } else { 774 #if defined(PETSC_USE_COMPLEX) 775 #if defined(PETSC_USE_REAL_SINGLE) 776 mumps->id.a_loc = (mumps_complex*)mumps->val; 777 #else 778 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 779 #endif 780 #else 781 mumps->id.a_loc = mumps->val; 782 #endif 783 } 784 PetscMUMPS_c(&mumps->id); 785 if (mumps->id.INFOG(1) < 0) { 786 if (mumps->id.INFO(1) == -13) { 787 if (mumps->id.INFO(2) < 0) { 788 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)); 789 } else { 790 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)); 791 } 792 } 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)); 793 } 794 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)); 795 796 (F)->assembled = PETSC_TRUE; 797 mumps->matstruc = SAME_NONZERO_PATTERN; 798 mumps->CleanUpMUMPS = PETSC_TRUE; 799 800 if (mumps->size > 1) { 801 PetscInt lsol_loc; 802 PetscScalar *sol_loc; 803 804 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 805 if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A; 806 else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A; 807 F_diag->assembled = PETSC_TRUE; 808 809 /* distributed solution; Create x_seq=sol_loc for repeated use */ 810 if (mumps->x_seq) { 811 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 812 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 813 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 814 } 815 lsol_loc = mumps->id.INFO(23); /* length of sol_loc */ 816 ierr = PetscMalloc2(lsol_loc,&sol_loc,lsol_loc,&mumps->id.isol_loc);CHKERRQ(ierr); 817 mumps->id.lsol_loc = lsol_loc; 818 #if defined(PETSC_USE_COMPLEX) 819 #if defined(PETSC_USE_REAL_SINGLE) 820 mumps->id.sol_loc = (mumps_complex*)sol_loc; 821 #else 822 mumps->id.sol_loc = (mumps_double_complex*)sol_loc; 823 #endif 824 #else 825 mumps->id.sol_loc = sol_loc; 826 #endif 827 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 828 } 829 PetscFunctionReturn(0); 830 } 831 832 /* Sets MUMPS options from the options database */ 833 #undef __FUNCT__ 834 #define __FUNCT__ "PetscSetMUMPSFromOptions" 835 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 836 { 837 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 838 PetscErrorCode ierr; 839 PetscInt icntl; 840 PetscBool flg; 841 842 PetscFunctionBegin; 843 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 844 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 845 if (flg) mumps->id.ICNTL(1) = icntl; 846 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); 847 if (flg) mumps->id.ICNTL(2) = icntl; 848 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); 849 if (flg) mumps->id.ICNTL(3) = icntl; 850 851 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 852 if (flg) mumps->id.ICNTL(4) = icntl; 853 if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 854 855 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); 856 if (flg) mumps->id.ICNTL(6) = icntl; 857 858 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); 859 if (flg) { 860 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"); 861 else mumps->id.ICNTL(7) = icntl; 862 } 863 864 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); 865 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);CHKERRQ(ierr); 866 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); 867 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); 868 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); 869 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); 870 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);CHKERRQ(ierr); 871 872 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); 873 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); 874 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); 875 if (mumps->id.ICNTL(24)) { 876 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 877 } 878 879 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); 880 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); 881 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",mumps->id.ICNTL(27),&mumps->id.ICNTL(27),NULL);CHKERRQ(ierr); 882 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); 883 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); 884 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); 885 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); 886 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);CHKERRQ(ierr); 887 888 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);CHKERRQ(ierr); 889 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);CHKERRQ(ierr); 890 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);CHKERRQ(ierr); 891 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);CHKERRQ(ierr); 892 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);CHKERRQ(ierr); 893 894 ierr = PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL); 895 PetscOptionsEnd(); 896 PetscFunctionReturn(0); 897 } 898 899 #undef __FUNCT__ 900 #define __FUNCT__ "PetscInitializeMUMPS" 901 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps) 902 { 903 PetscErrorCode ierr; 904 905 PetscFunctionBegin; 906 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid); 907 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);CHKERRQ(ierr); 908 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));CHKERRQ(ierr); 909 910 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 911 912 mumps->id.job = JOB_INIT; 913 mumps->id.par = 1; /* host participates factorizaton and solve */ 914 mumps->id.sym = mumps->sym; 915 PetscMUMPS_c(&mumps->id); 916 917 mumps->CleanUpMUMPS = PETSC_FALSE; 918 mumps->scat_rhs = NULL; 919 mumps->scat_sol = NULL; 920 921 /* set PETSc-MUMPS default options - override MUMPS default */ 922 mumps->id.ICNTL(3) = 0; 923 mumps->id.ICNTL(4) = 0; 924 if (mumps->size == 1) { 925 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 926 } else { 927 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 928 mumps->id.ICNTL(21) = 1; /* distributed solution */ 929 } 930 931 /* schur */ 932 mumps->id.size_schur = 0; 933 mumps->id.listvar_schur = NULL; 934 mumps->id.schur = NULL; 935 PetscFunctionReturn(0); 936 } 937 938 /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */ 939 #undef __FUNCT__ 940 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 941 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 942 { 943 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 944 PetscErrorCode ierr; 945 Vec b; 946 IS is_iden; 947 const PetscInt M = A->rmap->N; 948 949 PetscFunctionBegin; 950 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 951 952 /* Set MUMPS options from the options database */ 953 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 954 955 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 956 957 /* analysis phase */ 958 /*----------------*/ 959 mumps->id.job = JOB_FACTSYMBOLIC; 960 mumps->id.n = M; 961 switch (mumps->id.ICNTL(18)) { 962 case 0: /* centralized assembled matrix input */ 963 if (!mumps->myid) { 964 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 965 if (mumps->id.ICNTL(6)>1) { 966 #if defined(PETSC_USE_COMPLEX) 967 #if defined(PETSC_USE_REAL_SINGLE) 968 mumps->id.a = (mumps_complex*)mumps->val; 969 #else 970 mumps->id.a = (mumps_double_complex*)mumps->val; 971 #endif 972 #else 973 mumps->id.a = mumps->val; 974 #endif 975 } 976 if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */ 977 /* 978 PetscBool flag; 979 ierr = ISEqual(r,c,&flag);CHKERRQ(ierr); 980 if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm"); 981 ierr = ISView(r,PETSC_VIEWER_STDOUT_SELF); 982 */ 983 if (!mumps->myid) { 984 const PetscInt *idx; 985 PetscInt i,*perm_in; 986 987 ierr = PetscMalloc1(M,&perm_in);CHKERRQ(ierr); 988 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 989 990 mumps->id.perm_in = perm_in; 991 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 992 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 993 } 994 } 995 } 996 break; 997 case 3: /* distributed assembled matrix input (size>1) */ 998 mumps->id.nz_loc = mumps->nz; 999 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1000 if (mumps->id.ICNTL(6)>1) { 1001 #if defined(PETSC_USE_COMPLEX) 1002 #if defined(PETSC_USE_REAL_SINGLE) 1003 mumps->id.a_loc = (mumps_complex*)mumps->val; 1004 #else 1005 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1006 #endif 1007 #else 1008 mumps->id.a_loc = mumps->val; 1009 #endif 1010 } 1011 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1012 if (!mumps->myid) { 1013 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1014 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1015 } else { 1016 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1017 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1018 } 1019 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1020 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1021 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1022 ierr = VecDestroy(&b);CHKERRQ(ierr); 1023 break; 1024 } 1025 PetscMUMPS_c(&mumps->id); 1026 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)); 1027 1028 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1029 F->ops->solve = MatSolve_MUMPS; 1030 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1031 if (mumps->size > 1) { 1032 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 1033 } else { 1034 F->ops->matsolve = MatMatSolve_MUMPS; /* use MatMatSolve_MUMPS() for sequential matrices */ 1035 } 1036 PetscFunctionReturn(0); 1037 } 1038 1039 /* Note the Petsc r and c permutations are ignored */ 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 1042 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1043 { 1044 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1045 PetscErrorCode ierr; 1046 Vec b; 1047 IS is_iden; 1048 const PetscInt M = A->rmap->N; 1049 1050 PetscFunctionBegin; 1051 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1052 1053 /* Set MUMPS options from the options database */ 1054 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1055 1056 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1057 1058 /* analysis phase */ 1059 /*----------------*/ 1060 mumps->id.job = JOB_FACTSYMBOLIC; 1061 mumps->id.n = M; 1062 switch (mumps->id.ICNTL(18)) { 1063 case 0: /* centralized assembled matrix input */ 1064 if (!mumps->myid) { 1065 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1066 if (mumps->id.ICNTL(6)>1) { 1067 #if defined(PETSC_USE_COMPLEX) 1068 #if defined(PETSC_USE_REAL_SINGLE) 1069 mumps->id.a = (mumps_complex*)mumps->val; 1070 #else 1071 mumps->id.a = (mumps_double_complex*)mumps->val; 1072 #endif 1073 #else 1074 mumps->id.a = mumps->val; 1075 #endif 1076 } 1077 } 1078 break; 1079 case 3: /* distributed assembled matrix input (size>1) */ 1080 mumps->id.nz_loc = mumps->nz; 1081 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1082 if (mumps->id.ICNTL(6)>1) { 1083 #if defined(PETSC_USE_COMPLEX) 1084 #if defined(PETSC_USE_REAL_SINGLE) 1085 mumps->id.a_loc = (mumps_complex*)mumps->val; 1086 #else 1087 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1088 #endif 1089 #else 1090 mumps->id.a_loc = mumps->val; 1091 #endif 1092 } 1093 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1094 if (!mumps->myid) { 1095 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1096 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1097 } else { 1098 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1099 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1100 } 1101 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1102 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1103 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1104 ierr = VecDestroy(&b);CHKERRQ(ierr); 1105 break; 1106 } 1107 PetscMUMPS_c(&mumps->id); 1108 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)); 1109 1110 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1111 F->ops->solve = MatSolve_MUMPS; 1112 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1113 PetscFunctionReturn(0); 1114 } 1115 1116 /* Note the Petsc r permutation and factor info are ignored */ 1117 #undef __FUNCT__ 1118 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 1119 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 1120 { 1121 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1122 PetscErrorCode ierr; 1123 Vec b; 1124 IS is_iden; 1125 const PetscInt M = A->rmap->N; 1126 1127 PetscFunctionBegin; 1128 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1129 1130 /* Set MUMPS options from the options database */ 1131 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1132 1133 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1134 1135 /* analysis phase */ 1136 /*----------------*/ 1137 mumps->id.job = JOB_FACTSYMBOLIC; 1138 mumps->id.n = M; 1139 switch (mumps->id.ICNTL(18)) { 1140 case 0: /* centralized assembled matrix input */ 1141 if (!mumps->myid) { 1142 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1143 if (mumps->id.ICNTL(6)>1) { 1144 #if defined(PETSC_USE_COMPLEX) 1145 #if defined(PETSC_USE_REAL_SINGLE) 1146 mumps->id.a = (mumps_complex*)mumps->val; 1147 #else 1148 mumps->id.a = (mumps_double_complex*)mumps->val; 1149 #endif 1150 #else 1151 mumps->id.a = mumps->val; 1152 #endif 1153 } 1154 } 1155 break; 1156 case 3: /* distributed assembled matrix input (size>1) */ 1157 mumps->id.nz_loc = mumps->nz; 1158 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1159 if (mumps->id.ICNTL(6)>1) { 1160 #if defined(PETSC_USE_COMPLEX) 1161 #if defined(PETSC_USE_REAL_SINGLE) 1162 mumps->id.a_loc = (mumps_complex*)mumps->val; 1163 #else 1164 mumps->id.a_loc = (mumps_double_complex*)mumps->val; 1165 #endif 1166 #else 1167 mumps->id.a_loc = mumps->val; 1168 #endif 1169 } 1170 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1171 if (!mumps->myid) { 1172 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1173 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1174 } else { 1175 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1176 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1177 } 1178 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1179 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1180 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1181 ierr = VecDestroy(&b);CHKERRQ(ierr); 1182 break; 1183 } 1184 PetscMUMPS_c(&mumps->id); 1185 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)); 1186 1187 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1188 F->ops->solve = MatSolve_MUMPS; 1189 F->ops->solvetranspose = MatSolve_MUMPS; 1190 F->ops->matsolve = 0; /* use MatMatSolve_Basic() until mumps supports distributed rhs */ 1191 #if !defined(PETSC_USE_COMPLEX) 1192 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1193 #else 1194 F->ops->getinertia = NULL; 1195 #endif 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "MatView_MUMPS" 1201 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1202 { 1203 PetscErrorCode ierr; 1204 PetscBool iascii; 1205 PetscViewerFormat format; 1206 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1207 1208 PetscFunctionBegin; 1209 /* check if matrix is mumps type */ 1210 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1211 1212 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1213 if (iascii) { 1214 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1215 if (format == PETSC_VIEWER_ASCII_INFO) { 1216 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1217 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1218 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1219 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));CHKERRQ(ierr); 1220 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));CHKERRQ(ierr); 1221 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));CHKERRQ(ierr); 1222 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));CHKERRQ(ierr); 1223 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));CHKERRQ(ierr); 1224 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));CHKERRQ(ierr); 1225 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));CHKERRQ(ierr); 1226 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));CHKERRQ(ierr); 1227 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));CHKERRQ(ierr); 1228 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));CHKERRQ(ierr); 1229 if (mumps->id.ICNTL(11)>0) { 1230 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));CHKERRQ(ierr); 1231 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));CHKERRQ(ierr); 1232 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));CHKERRQ(ierr); 1233 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));CHKERRQ(ierr); 1234 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));CHKERRQ(ierr); 1235 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));CHKERRQ(ierr); 1236 } 1237 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));CHKERRQ(ierr); 1238 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));CHKERRQ(ierr); 1239 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));CHKERRQ(ierr); 1240 /* ICNTL(15-17) not used */ 1241 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));CHKERRQ(ierr); 1242 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));CHKERRQ(ierr); 1243 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));CHKERRQ(ierr); 1244 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (somumpstion struct): %d \n",mumps->id.ICNTL(21));CHKERRQ(ierr); 1245 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));CHKERRQ(ierr); 1246 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));CHKERRQ(ierr); 1247 1248 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));CHKERRQ(ierr); 1249 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));CHKERRQ(ierr); 1250 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));CHKERRQ(ierr); 1251 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));CHKERRQ(ierr); 1252 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));CHKERRQ(ierr); 1253 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));CHKERRQ(ierr); 1254 1255 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));CHKERRQ(ierr); 1256 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));CHKERRQ(ierr); 1257 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));CHKERRQ(ierr); 1258 1259 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));CHKERRQ(ierr); 1260 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));CHKERRQ(ierr); 1261 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absomumpste pivoting threshold): %g \n",mumps->id.CNTL(3));CHKERRQ(ierr); 1262 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (vamumpse of static pivoting): %g \n",mumps->id.CNTL(4));CHKERRQ(ierr); 1263 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));CHKERRQ(ierr); 1264 1265 /* infomation local to each processor */ 1266 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1267 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1268 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));CHKERRQ(ierr); 1269 ierr = PetscViewerFlush(viewer); 1270 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1271 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));CHKERRQ(ierr); 1272 ierr = PetscViewerFlush(viewer); 1273 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1274 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));CHKERRQ(ierr); 1275 ierr = PetscViewerFlush(viewer); 1276 1277 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1278 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));CHKERRQ(ierr); 1279 ierr = PetscViewerFlush(viewer); 1280 1281 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1282 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));CHKERRQ(ierr); 1283 ierr = PetscViewerFlush(viewer); 1284 1285 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1286 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));CHKERRQ(ierr); 1287 ierr = PetscViewerFlush(viewer); 1288 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1289 1290 if (!mumps->myid) { /* information from the host */ 1291 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));CHKERRQ(ierr); 1292 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));CHKERRQ(ierr); 1293 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));CHKERRQ(ierr); 1294 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); 1295 1296 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));CHKERRQ(ierr); 1297 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));CHKERRQ(ierr); 1298 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));CHKERRQ(ierr); 1299 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));CHKERRQ(ierr); 1300 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));CHKERRQ(ierr); 1301 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));CHKERRQ(ierr); 1302 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));CHKERRQ(ierr); 1303 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));CHKERRQ(ierr); 1304 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));CHKERRQ(ierr); 1305 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));CHKERRQ(ierr); 1306 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));CHKERRQ(ierr); 1307 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));CHKERRQ(ierr); 1308 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));CHKERRQ(ierr); 1309 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); 1310 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); 1311 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); 1312 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); 1313 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));CHKERRQ(ierr); 1314 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); 1315 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); 1316 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));CHKERRQ(ierr); 1317 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));CHKERRQ(ierr); 1318 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));CHKERRQ(ierr); 1319 ierr = PetscViewerASCIIPrintf(viewer," INFOG(28) (after factorization: number of null pivots encountered): %d\n",mumps->id.INFOG(28));CHKERRQ(ierr); 1320 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); 1321 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); 1322 ierr = PetscViewerASCIIPrintf(viewer," INFOG(32) (after analysis: type of analysis done): %d\n",mumps->id.INFOG(32));CHKERRQ(ierr); 1323 ierr = PetscViewerASCIIPrintf(viewer," INFOG(33) (value used for ICNTL(8)): %d\n",mumps->id.INFOG(33));CHKERRQ(ierr); 1324 ierr = PetscViewerASCIIPrintf(viewer," INFOG(34) (exponent of the determinant if determinant is requested): %d\n",mumps->id.INFOG(34));CHKERRQ(ierr); 1325 } 1326 } 1327 } 1328 PetscFunctionReturn(0); 1329 } 1330 1331 #undef __FUNCT__ 1332 #define __FUNCT__ "MatGetInfo_MUMPS" 1333 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1334 { 1335 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1336 1337 PetscFunctionBegin; 1338 info->block_size = 1.0; 1339 info->nz_allocated = mumps->id.INFOG(20); 1340 info->nz_used = mumps->id.INFOG(20); 1341 info->nz_unneeded = 0.0; 1342 info->assemblies = 0.0; 1343 info->mallocs = 0.0; 1344 info->memory = 0.0; 1345 info->fill_ratio_given = 0; 1346 info->fill_ratio_needed = 0; 1347 info->factor_mallocs = 0; 1348 PetscFunctionReturn(0); 1349 } 1350 1351 /* -------------------------------------------------------------------------------------------*/ 1352 #undef __FUNCT__ 1353 #define __FUNCT__ "MatMumpsSetSchurIndices_MUMPS" 1354 PetscErrorCode MatMumpsSetSchurIndices_MUMPS(Mat F,PetscInt size,PetscInt idxs[]) 1355 { 1356 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1357 PetscErrorCode ierr; 1358 1359 PetscFunctionBegin; 1360 if (mumps->size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MUMPS parallel computation of Schur complements not yet supported from PETSc\n"); 1361 if (mumps->id.size_schur != size) { 1362 ierr = PetscFree2(mumps->id.listvar_schur,mumps->id.schur);CHKERRQ(ierr); 1363 mumps->id.size_schur = size; 1364 mumps->id.schur_lld = size; 1365 ierr = PetscMalloc2(size,&mumps->id.listvar_schur,size*size,&mumps->id.schur);CHKERRQ(ierr); 1366 } 1367 ierr = PetscMemcpy(mumps->id.listvar_schur,idxs,size*sizeof(PetscInt));CHKERRQ(ierr); 1368 if (F->factortype == MAT_FACTOR_LU) { 1369 mumps->id.ICNTL(19) = 3; /* return full matrix */ 1370 } else { 1371 mumps->id.ICNTL(19) = 2; /* return lower triangular part */ 1372 } 1373 PetscFunctionReturn(0); 1374 } 1375 1376 #undef __FUNCT__ 1377 #define __FUNCT__ "MatMumpsSetSchurIndices" 1378 /*@ 1379 MatMumpsSetSchurIndices - Set indices defining the Schur complement that MUMPS will compute during the factorization steps 1380 1381 Logically Collective on Mat 1382 1383 Input Parameters: 1384 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1385 . size - size of the Schur complement indices 1386 - idxs[] - array of Schur complement indices 1387 1388 Notes: 1389 The user has to free the array idxs[] since it is copied by the routine. 1390 Currently implemented for sequential matrices 1391 1392 Level: advanced 1393 1394 References: MUMPS Users' Guide 1395 1396 .seealso: MatGetFactor() 1397 @*/ 1398 PetscErrorCode MatMumpsSetSchurIndices(Mat F,PetscInt size,PetscInt idxs[]) 1399 { 1400 PetscErrorCode ierr; 1401 1402 PetscFunctionBegin; 1403 PetscValidIntPointer(idxs,3); 1404 ierr = PetscTryMethod(F,"MatMumpsSetSchurIndices_C",(Mat,PetscInt,PetscInt[]),(F,size,idxs));CHKERRQ(ierr); 1405 PetscFunctionReturn(0); 1406 } 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 PetscErrorCode ierr; 1417 1418 PetscFunctionBegin; 1419 if (mumps->id.job != JOB_FACTNUMERIC) { 1420 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_WRONG,"Numerical factorization phase not yet performed! You should call MatFactorSymbolic/Numeric before"); 1421 } else if (mumps->id.size_schur == 0) { 1422 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_WRONG,"Schur indices not set! You should call MatMumpsSetSchurIndices before"); 1423 } 1424 1425 ierr = MatCreate(PetscObjectComm((PetscObject)F),&St);CHKERRQ(ierr); 1426 ierr = MatSetSizes(St,PETSC_DECIDE,PETSC_DECIDE,mumps->id.size_schur,mumps->id.size_schur);CHKERRQ(ierr); 1427 ierr = MatSetType(St,MATDENSE);CHKERRQ(ierr); 1428 ierr = MatSetUp(St);CHKERRQ(ierr); 1429 ierr = MatDenseGetArray(St,&array);CHKERRQ(ierr); 1430 if (mumps->sym == 0) { /* MUMPS always return a full matrix */ 1431 if (mumps->id.ICNTL(19) == 1) { /* stored by rows */ 1432 PetscInt i,j,N=mumps->id.size_schur; 1433 for (i=0;i<N;i++) { 1434 for (j=0;j<N;j++) { 1435 PetscScalar val = mumps->id.schur[i*N+j]; 1436 array[j*N+i] = val; 1437 } 1438 } 1439 } else { /* stored by columns */ 1440 ierr = PetscMemcpy(array,mumps->id.schur,mumps->id.size_schur*mumps->id.size_schur*sizeof(PetscScalar));CHKERRQ(ierr); 1441 } 1442 } else { /* either full or lower-triangular (not packed) */ 1443 if (mumps->id.ICNTL(19) == 2) { /* lower triangular stored by columns */ 1444 PetscInt i,j,N=mumps->id.size_schur; 1445 for (i=0;i<N;i++) { 1446 for (j=i;j<N;j++) { 1447 PetscScalar val = mumps->id.schur[i*N+j]; 1448 array[i*N+j] = val; 1449 } 1450 for (j=i;j<N;j++) { 1451 PetscScalar val = mumps->id.schur[i*N+j]; 1452 array[j*N+i] = val; 1453 } 1454 } 1455 } else if (mumps->id.ICNTL(19) == 3) { /* full matrix */ 1456 ierr = PetscMemcpy(array,mumps->id.schur,mumps->id.size_schur*mumps->id.size_schur*sizeof(PetscScalar));CHKERRQ(ierr); 1457 } else { /* ICNTL(19) == 1 lower triangular stored by rows */ 1458 PetscInt i,j,N=mumps->id.size_schur; 1459 for (i=0;i<N;i++) { 1460 for (j=0;j<i+1;j++) { 1461 PetscScalar val = mumps->id.schur[i*N+j]; 1462 array[i*N+j] = val; 1463 } 1464 for (j=0;j<i+1;j++) { 1465 PetscScalar val = mumps->id.schur[i*N+j]; 1466 array[j*N+i] = val; 1467 } 1468 } 1469 } 1470 } 1471 ierr = MatDenseRestoreArray(St,&array);CHKERRQ(ierr); 1472 *S = St; 1473 PetscFunctionReturn(0); 1474 } 1475 1476 #undef __FUNCT__ 1477 #define __FUNCT__ "MatMumpsGetSchurComplement" 1478 /*@ 1479 MatMumpsGetSchurComplement - Get Schur complement matrix computed by MUMPS during the factorization step 1480 1481 Logically Collective on Mat 1482 1483 Input Parameters: 1484 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1485 . *S - location where to return the Schur complement (MATDENSE) 1486 1487 Notes: 1488 Currently implemented for sequential matrices 1489 1490 Level: advanced 1491 1492 References: MUMPS Users' Guide 1493 1494 .seealso: MatGetFactor() 1495 @*/ 1496 PetscErrorCode MatMumpsGetSchurComplement(Mat F,Mat* S) 1497 { 1498 PetscErrorCode ierr; 1499 1500 PetscFunctionBegin; 1501 ierr = PetscTryMethod(F,"MatMumpsGetSchurComplement_C",(Mat,Mat*),(F,S));CHKERRQ(ierr); 1502 PetscFunctionReturn(0); 1503 } 1504 1505 /* -------------------------------------------------------------------------------------------*/ 1506 #undef __FUNCT__ 1507 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 1508 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 1509 { 1510 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1511 1512 PetscFunctionBegin; 1513 mumps->id.ICNTL(icntl) = ival; 1514 PetscFunctionReturn(0); 1515 } 1516 1517 #undef __FUNCT__ 1518 #define __FUNCT__ "MatMumpsGetIcntl_MUMPS" 1519 PetscErrorCode MatMumpsGetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt *ival) 1520 { 1521 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1522 1523 PetscFunctionBegin; 1524 *ival = mumps->id.ICNTL(icntl); 1525 PetscFunctionReturn(0); 1526 } 1527 1528 #undef __FUNCT__ 1529 #define __FUNCT__ "MatMumpsSetIcntl" 1530 /*@ 1531 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1532 1533 Logically Collective on Mat 1534 1535 Input Parameters: 1536 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1537 . icntl - index of MUMPS parameter array ICNTL() 1538 - ival - value of MUMPS ICNTL(icntl) 1539 1540 Options Database: 1541 . -mat_mumps_icntl_<icntl> <ival> 1542 1543 Level: beginner 1544 1545 References: MUMPS Users' Guide 1546 1547 .seealso: MatGetFactor() 1548 @*/ 1549 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 1550 { 1551 PetscErrorCode ierr; 1552 1553 PetscFunctionBegin; 1554 PetscValidLogicalCollectiveInt(F,icntl,2); 1555 PetscValidLogicalCollectiveInt(F,ival,3); 1556 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 1557 PetscFunctionReturn(0); 1558 } 1559 1560 #undef __FUNCT__ 1561 #define __FUNCT__ "MatMumpsGetIcntl" 1562 /*@ 1563 MatMumpsGetIcntl - Get MUMPS parameter ICNTL() 1564 1565 Logically Collective on Mat 1566 1567 Input Parameters: 1568 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1569 - icntl - index of MUMPS parameter array ICNTL() 1570 1571 Output Parameter: 1572 . ival - value of MUMPS ICNTL(icntl) 1573 1574 Level: beginner 1575 1576 References: MUMPS Users' Guide 1577 1578 .seealso: MatGetFactor() 1579 @*/ 1580 PetscErrorCode MatMumpsGetIcntl(Mat F,PetscInt icntl,PetscInt *ival) 1581 { 1582 PetscErrorCode ierr; 1583 1584 PetscFunctionBegin; 1585 PetscValidLogicalCollectiveInt(F,icntl,2); 1586 PetscValidIntPointer(ival,3); 1587 ierr = PetscTryMethod(F,"MatMumpsGetIcntl_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1588 PetscFunctionReturn(0); 1589 } 1590 1591 /* -------------------------------------------------------------------------------------------*/ 1592 #undef __FUNCT__ 1593 #define __FUNCT__ "MatMumpsSetCntl_MUMPS" 1594 PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val) 1595 { 1596 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1597 1598 PetscFunctionBegin; 1599 mumps->id.CNTL(icntl) = val; 1600 PetscFunctionReturn(0); 1601 } 1602 1603 #undef __FUNCT__ 1604 #define __FUNCT__ "MatMumpsGetCntl_MUMPS" 1605 PetscErrorCode MatMumpsGetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal *val) 1606 { 1607 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1608 1609 PetscFunctionBegin; 1610 *val = mumps->id.CNTL(icntl); 1611 PetscFunctionReturn(0); 1612 } 1613 1614 #undef __FUNCT__ 1615 #define __FUNCT__ "MatMumpsSetCntl" 1616 /*@ 1617 MatMumpsSetCntl - Set MUMPS parameter CNTL() 1618 1619 Logically Collective on Mat 1620 1621 Input Parameters: 1622 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1623 . icntl - index of MUMPS parameter array CNTL() 1624 - val - value of MUMPS CNTL(icntl) 1625 1626 Options Database: 1627 . -mat_mumps_cntl_<icntl> <val> 1628 1629 Level: beginner 1630 1631 References: MUMPS Users' Guide 1632 1633 .seealso: MatGetFactor() 1634 @*/ 1635 PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val) 1636 { 1637 PetscErrorCode ierr; 1638 1639 PetscFunctionBegin; 1640 PetscValidLogicalCollectiveInt(F,icntl,2); 1641 PetscValidLogicalCollectiveReal(F,val,3); 1642 ierr = PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));CHKERRQ(ierr); 1643 PetscFunctionReturn(0); 1644 } 1645 1646 #undef __FUNCT__ 1647 #define __FUNCT__ "MatMumpsGetCntl" 1648 /*@ 1649 MatMumpsGetCntl - Get MUMPS parameter CNTL() 1650 1651 Logically Collective on Mat 1652 1653 Input Parameters: 1654 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1655 - icntl - index of MUMPS parameter array CNTL() 1656 1657 Output Parameter: 1658 . val - value of MUMPS CNTL(icntl) 1659 1660 Level: beginner 1661 1662 References: MUMPS Users' Guide 1663 1664 .seealso: MatGetFactor() 1665 @*/ 1666 PetscErrorCode MatMumpsGetCntl(Mat F,PetscInt icntl,PetscReal *val) 1667 { 1668 PetscErrorCode ierr; 1669 1670 PetscFunctionBegin; 1671 PetscValidLogicalCollectiveInt(F,icntl,2); 1672 PetscValidRealPointer(val,3); 1673 ierr = PetscTryMethod(F,"MatMumpsGetCntl_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1674 PetscFunctionReturn(0); 1675 } 1676 1677 #undef __FUNCT__ 1678 #define __FUNCT__ "MatMumpsGetInfo_MUMPS" 1679 PetscErrorCode MatMumpsGetInfo_MUMPS(Mat F,PetscInt icntl,PetscInt *info) 1680 { 1681 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1682 1683 PetscFunctionBegin; 1684 *info = mumps->id.INFO(icntl); 1685 PetscFunctionReturn(0); 1686 } 1687 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "MatMumpsGetInfog_MUMPS" 1690 PetscErrorCode MatMumpsGetInfog_MUMPS(Mat F,PetscInt icntl,PetscInt *infog) 1691 { 1692 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1693 1694 PetscFunctionBegin; 1695 *infog = mumps->id.INFOG(icntl); 1696 PetscFunctionReturn(0); 1697 } 1698 1699 #undef __FUNCT__ 1700 #define __FUNCT__ "MatMumpsGetRinfo_MUMPS" 1701 PetscErrorCode MatMumpsGetRinfo_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfo) 1702 { 1703 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1704 1705 PetscFunctionBegin; 1706 *rinfo = mumps->id.RINFO(icntl); 1707 PetscFunctionReturn(0); 1708 } 1709 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "MatMumpsGetRinfog_MUMPS" 1712 PetscErrorCode MatMumpsGetRinfog_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfog) 1713 { 1714 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1715 1716 PetscFunctionBegin; 1717 *rinfog = mumps->id.RINFOG(icntl); 1718 PetscFunctionReturn(0); 1719 } 1720 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "MatMumpsGetInfo" 1723 /*@ 1724 MatMumpsGetInfo - Get MUMPS parameter INFO() 1725 1726 Logically Collective on Mat 1727 1728 Input Parameters: 1729 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1730 - icntl - index of MUMPS parameter array INFO() 1731 1732 Output Parameter: 1733 . ival - value of MUMPS INFO(icntl) 1734 1735 Level: beginner 1736 1737 References: MUMPS Users' Guide 1738 1739 .seealso: MatGetFactor() 1740 @*/ 1741 PetscErrorCode MatMumpsGetInfo(Mat F,PetscInt icntl,PetscInt *ival) 1742 { 1743 PetscErrorCode ierr; 1744 1745 PetscFunctionBegin; 1746 PetscValidIntPointer(ival,3); 1747 ierr = PetscTryMethod(F,"MatMumpsGetInfo_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1748 PetscFunctionReturn(0); 1749 } 1750 1751 #undef __FUNCT__ 1752 #define __FUNCT__ "MatMumpsGetInfog" 1753 /*@ 1754 MatMumpsGetInfog - Get MUMPS parameter INFOG() 1755 1756 Logically Collective on Mat 1757 1758 Input Parameters: 1759 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1760 - icntl - index of MUMPS parameter array INFOG() 1761 1762 Output Parameter: 1763 . ival - value of MUMPS INFOG(icntl) 1764 1765 Level: beginner 1766 1767 References: MUMPS Users' Guide 1768 1769 .seealso: MatGetFactor() 1770 @*/ 1771 PetscErrorCode MatMumpsGetInfog(Mat F,PetscInt icntl,PetscInt *ival) 1772 { 1773 PetscErrorCode ierr; 1774 1775 PetscFunctionBegin; 1776 PetscValidIntPointer(ival,3); 1777 ierr = PetscTryMethod(F,"MatMumpsGetInfog_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 1778 PetscFunctionReturn(0); 1779 } 1780 1781 #undef __FUNCT__ 1782 #define __FUNCT__ "MatMumpsGetRinfo" 1783 /*@ 1784 MatMumpsGetRinfo - Get MUMPS parameter RINFO() 1785 1786 Logically Collective on Mat 1787 1788 Input Parameters: 1789 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1790 - icntl - index of MUMPS parameter array RINFO() 1791 1792 Output Parameter: 1793 . val - value of MUMPS RINFO(icntl) 1794 1795 Level: beginner 1796 1797 References: MUMPS Users' Guide 1798 1799 .seealso: MatGetFactor() 1800 @*/ 1801 PetscErrorCode MatMumpsGetRinfo(Mat F,PetscInt icntl,PetscReal *val) 1802 { 1803 PetscErrorCode ierr; 1804 1805 PetscFunctionBegin; 1806 PetscValidRealPointer(val,3); 1807 ierr = PetscTryMethod(F,"MatMumpsGetRinfo_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1808 PetscFunctionReturn(0); 1809 } 1810 1811 #undef __FUNCT__ 1812 #define __FUNCT__ "MatMumpsGetRinfog" 1813 /*@ 1814 MatMumpsGetRinfog - Get MUMPS parameter RINFOG() 1815 1816 Logically Collective on Mat 1817 1818 Input Parameters: 1819 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1820 - icntl - index of MUMPS parameter array RINFOG() 1821 1822 Output Parameter: 1823 . val - value of MUMPS RINFOG(icntl) 1824 1825 Level: beginner 1826 1827 References: MUMPS Users' Guide 1828 1829 .seealso: MatGetFactor() 1830 @*/ 1831 PetscErrorCode MatMumpsGetRinfog(Mat F,PetscInt icntl,PetscReal *val) 1832 { 1833 PetscErrorCode ierr; 1834 1835 PetscFunctionBegin; 1836 PetscValidRealPointer(val,3); 1837 ierr = PetscTryMethod(F,"MatMumpsGetRinfog_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 1838 PetscFunctionReturn(0); 1839 } 1840 1841 /*MC 1842 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1843 distributed and sequential matrices via the external package MUMPS. 1844 1845 Works with MATAIJ and MATSBAIJ matrices 1846 1847 Options Database Keys: 1848 + -mat_mumps_icntl_4 <0,...,4> - print level 1849 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1850 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guidec) 1851 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1852 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1853 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1854 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1855 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1856 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1857 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1858 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1859 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1860 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1861 1862 Level: beginner 1863 1864 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1865 1866 M*/ 1867 1868 #undef __FUNCT__ 1869 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1870 static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1871 { 1872 PetscFunctionBegin; 1873 *type = MATSOLVERMUMPS; 1874 PetscFunctionReturn(0); 1875 } 1876 1877 /* MatGetFactor for Seq and MPI AIJ matrices */ 1878 #undef __FUNCT__ 1879 #define __FUNCT__ "MatGetFactor_aij_mumps" 1880 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 1881 { 1882 Mat B; 1883 PetscErrorCode ierr; 1884 Mat_MUMPS *mumps; 1885 PetscBool isSeqAIJ; 1886 1887 PetscFunctionBegin; 1888 /* Create the factorization matrix */ 1889 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1890 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1891 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1892 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1893 if (isSeqAIJ) { 1894 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 1895 } else { 1896 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 1897 } 1898 1899 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1900 1901 B->ops->view = MatView_MUMPS; 1902 B->ops->getinfo = MatGetInfo_MUMPS; 1903 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1904 1905 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1906 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1907 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1908 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1909 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1910 1911 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1912 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1913 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1914 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1915 1916 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 1917 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 1918 if (ftype == MAT_FACTOR_LU) { 1919 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1920 B->factortype = MAT_FACTOR_LU; 1921 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1922 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1923 mumps->sym = 0; 1924 } else { 1925 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1926 B->factortype = MAT_FACTOR_CHOLESKY; 1927 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1928 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1929 if (A->spd_set && A->spd) mumps->sym = 1; 1930 else mumps->sym = 2; 1931 } 1932 1933 mumps->isAIJ = PETSC_TRUE; 1934 mumps->Destroy = B->ops->destroy; 1935 B->ops->destroy = MatDestroy_MUMPS; 1936 B->spptr = (void*)mumps; 1937 1938 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 1939 1940 *F = B; 1941 PetscFunctionReturn(0); 1942 } 1943 1944 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 1945 #undef __FUNCT__ 1946 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 1947 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1948 { 1949 Mat B; 1950 PetscErrorCode ierr; 1951 Mat_MUMPS *mumps; 1952 PetscBool isSeqSBAIJ; 1953 1954 PetscFunctionBegin; 1955 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1956 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"); 1957 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1958 /* Create the factorization matrix */ 1959 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1960 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1961 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1962 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 1963 if (isSeqSBAIJ) { 1964 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 1965 1966 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1967 } else { 1968 ierr = MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);CHKERRQ(ierr); 1969 1970 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1971 } 1972 1973 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1974 B->ops->view = MatView_MUMPS; 1975 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 1976 1977 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1978 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 1979 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 1980 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 1981 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 1982 1983 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 1984 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 1985 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 1986 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 1987 1988 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 1989 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 1990 1991 B->factortype = MAT_FACTOR_CHOLESKY; 1992 if (A->spd_set && A->spd) mumps->sym = 1; 1993 else mumps->sym = 2; 1994 1995 mumps->isAIJ = PETSC_FALSE; 1996 mumps->Destroy = B->ops->destroy; 1997 B->ops->destroy = MatDestroy_MUMPS; 1998 B->spptr = (void*)mumps; 1999 2000 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 2001 2002 *F = B; 2003 PetscFunctionReturn(0); 2004 } 2005 2006 #undef __FUNCT__ 2007 #define __FUNCT__ "MatGetFactor_baij_mumps" 2008 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 2009 { 2010 Mat B; 2011 PetscErrorCode ierr; 2012 Mat_MUMPS *mumps; 2013 PetscBool isSeqBAIJ; 2014 2015 PetscFunctionBegin; 2016 /* Create the factorization matrix */ 2017 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 2018 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 2019 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 2020 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2021 if (isSeqBAIJ) { 2022 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);CHKERRQ(ierr); 2023 } else { 2024 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 2025 } 2026 2027 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 2028 if (ftype == MAT_FACTOR_LU) { 2029 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 2030 B->factortype = MAT_FACTOR_LU; 2031 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 2032 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 2033 mumps->sym = 0; 2034 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 2035 2036 B->ops->view = MatView_MUMPS; 2037 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 2038 2039 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 2040 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 2041 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 2042 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 2043 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 2044 2045 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 2046 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 2047 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 2048 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 2049 2050 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 2051 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 2052 2053 mumps->isAIJ = PETSC_TRUE; 2054 mumps->Destroy = B->ops->destroy; 2055 B->ops->destroy = MatDestroy_MUMPS; 2056 B->spptr = (void*)mumps; 2057 2058 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 2059 2060 *F = B; 2061 PetscFunctionReturn(0); 2062 } 2063 2064 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 2065 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat,MatFactorType,Mat*); 2066 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 2067 2068 #undef __FUNCT__ 2069 #define __FUNCT__ "MatSolverPackageRegister_MUMPS" 2070 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MUMPS(void) 2071 { 2072 PetscErrorCode ierr; 2073 2074 PetscFunctionBegin; 2075 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2076 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2077 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2078 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2079 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPISBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2080 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2081 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2082 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2083 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2084 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2085 PetscFunctionReturn(0); 2086 } 2087 2088