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