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