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