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