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