1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the MUMPS sparse solver 5 */ 6 #include "../src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 7 #include "../src/mat/impls/aij/mpi/mpiaij.h" 8 #include "../src/mat/impls/sbaij/seq/sbaij.h" 9 #include "../src/mat/impls/sbaij/mpi/mpisbaij.h" 10 #include "../src/mat/impls/baij/seq/baij.h" 11 #include "../src/mat/impls/baij/mpi/mpibaij.h" 12 13 EXTERN_C_BEGIN 14 #if defined(PETSC_USE_COMPLEX) 15 #include "zmumps_c.h" 16 #else 17 #include "dmumps_c.h" 18 #endif 19 EXTERN_C_END 20 #define JOB_INIT -1 21 #define JOB_END -2 22 /* macros s.t. indices match MUMPS documentation */ 23 #define ICNTL(I) icntl[(I)-1] 24 #define CNTL(I) cntl[(I)-1] 25 #define INFOG(I) infog[(I)-1] 26 #define INFO(I) info[(I)-1] 27 #define RINFOG(I) rinfog[(I)-1] 28 #define RINFO(I) rinfo[(I)-1] 29 30 typedef struct { 31 #if defined(PETSC_USE_COMPLEX) 32 ZMUMPS_STRUC_C id; 33 #else 34 DMUMPS_STRUC_C id; 35 #endif 36 MatStructure matstruc; 37 PetscMPIInt myid,size; 38 PetscInt *irn,*jcn,nz,sym,nSolve; 39 PetscScalar *val; 40 MPI_Comm comm_mumps; 41 VecScatter scat_rhs, scat_sol; 42 PetscTruth isAIJ,CleanUpMUMPS,mumpsview; 43 Vec b_seq,x_seq; 44 PetscErrorCode (*MatDestroy)(Mat); 45 PetscErrorCode (*ConvertToTriples)(Mat, int, PetscTruth, int*, int**, int**, PetscScalar**); 46 } Mat_MUMPS; 47 48 EXTERN PetscErrorCode MatDuplicate_MUMPS(Mat,MatDuplicateOption,Mat*); 49 50 51 /* MatConvertToTriples_A_B */ 52 /*convert Petsc matrix to triples: row[nz], col[nz], val[nz] */ 53 /* 54 input: 55 A - matrix in aij,baij or sbaij (bs=1) format 56 shift - 0: C style output triple; 1: Fortran style output triple. 57 valOnly - FALSE: spaces are allocated and values are set for the triple 58 TRUE: only the values in v array are updated 59 output: 60 nnz - dim of r, c, and v (number of local nonzero entries of A) 61 r, c, v - row and col index, matrix values (matrix triples) 62 */ 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "MatConvertToTriples_seqaij_seqaij" 66 PetscErrorCode MatConvertToTriples_seqaij_seqaij(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 67 { 68 const PetscInt *ai,*aj,*ajj,M=A->rmap->n;; 69 PetscInt nz,rnz,i,j; 70 PetscErrorCode ierr; 71 PetscInt *row,*col; 72 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 73 74 PetscFunctionBegin; 75 *v=aa->a; 76 if (!valOnly){ 77 nz = aa->nz; ai = aa->i; aj = aa->j; 78 *nnz = nz; 79 ierr = PetscMalloc(nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 80 ierr = PetscMalloc(nz*sizeof(PetscInt), &col);CHKERRQ(ierr); 81 nz = 0; 82 for(i=0; i<M; i++) { 83 rnz = ai[i+1] - ai[i]; 84 ajj = aj + ai[i]; 85 for(j=0; j<rnz; j++) { 86 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 87 } 88 } 89 *r = row; *c = col; 90 } 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "MatConvertToTriples_seqbaij_seqaij" 96 PetscErrorCode MatConvertToTriples_seqbaij_seqaij(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 97 { 98 Mat_SeqBAIJ *aa=(Mat_SeqBAIJ*)A->data; 99 const PetscInt *ai,*aj,*ajj,bs=A->rmap->bs,bs2=aa->bs2,M=A->rmap->N/bs; 100 const PetscScalar *av,*v1; 101 PetscScalar *val; 102 PetscInt nz,idx=0,rnz,i,j,k,m,ii; 103 PetscErrorCode ierr; 104 PetscInt *row,*col; 105 106 PetscFunctionBegin; 107 ai = aa->i; aj = aa->j; av = aa->a; 108 if (!valOnly){ 109 nz = bs2*aa->nz; 110 *nnz = nz; 111 ierr = PetscMalloc(nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 112 ierr = PetscMalloc(nz*sizeof(PetscInt), &col);CHKERRQ(ierr); 113 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 114 for(i=0; i<M; i++) { 115 ii = 0; 116 ajj = aj + ai[i]; 117 v1 = av + bs2*ai[i]; 118 rnz = ai[i+1] - ai[i]; 119 for(k=0; k<rnz; k++) { 120 for(j=0; j<bs; j++) { 121 for(m=0; m<bs; m++) { 122 row[idx] = i*bs + m + shift; 123 col[idx] = bs*(ajj[k]) + j + shift; 124 val[idx++] = v1[ii++]; 125 } 126 } 127 } 128 } 129 *r = row; *c = col; *v = val; 130 } else { 131 row = *r; col = *c; val = *v; 132 for(i=0; i<M; i++) { 133 ii = 0; 134 v1 = av + bs2*ai[i]; 135 rnz = ai[i+1] - ai[i]; 136 for(k=0; k<rnz; k++){ 137 for(j=0; j<bs; j++) { 138 for(m=0; m<bs; m++) { 139 val[idx++] = v1[ii++]; 140 } 141 } 142 } 143 } 144 } 145 PetscFunctionReturn(0); 146 } 147 148 #undef __FUNCT__ 149 #define __FUNCT__ "MatConvertToTriples_seqsbaij_seqsbaij" 150 PetscErrorCode MatConvertToTriples_seqsbaij_seqsbaij(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 151 { 152 const PetscInt *ai, *aj,*ajj,M=A->rmap->n; 153 PetscInt nz,rnz,i,j; 154 PetscErrorCode ierr; 155 PetscInt *row,*col; 156 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data; 157 158 PetscFunctionBegin; 159 if (!valOnly){ 160 nz = aa->nz;ai=aa->i; aj=aa->j;*v=aa->a; 161 *nnz = nz; 162 ierr = PetscMalloc(nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 163 ierr = PetscMalloc(nz*sizeof(PetscInt), &col);CHKERRQ(ierr); 164 nz = 0; 165 for(i=0; i<M; i++) { 166 rnz = ai[i+1] - ai[i]; 167 ajj = aj + ai[i]; 168 for(j=0; j<rnz; j++) { 169 row[nz] = i+shift; col[nz++] = ajj[j] + shift; 170 } 171 } 172 *r = row; *c = col; 173 } 174 PetscFunctionReturn(0); 175 } 176 177 #undef __FUNCT__ 178 #define __FUNCT__ "MatConvertToTriples_seqaij_seqsbaij" 179 PetscErrorCode MatConvertToTriples_seqaij_seqsbaij(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 180 { 181 const PetscInt *ai,*aj,*ajj,*adiag,M=A->rmap->n; 182 PetscInt nz,rnz,i,j; 183 const PetscScalar *av,*v1; 184 PetscScalar *val; 185 PetscErrorCode ierr; 186 PetscInt *row,*col; 187 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)A->data; 188 189 PetscFunctionBegin; 190 ai=aa->i; aj=aa->j;av=aa->a; 191 adiag=aa->diag; 192 if (!valOnly){ 193 nz = M + (aa->nz-M)/2; 194 *nnz = nz; 195 ierr = PetscMalloc(nz*sizeof(PetscInt), &row);CHKERRQ(ierr); 196 ierr = PetscMalloc(nz*sizeof(PetscInt), &col);CHKERRQ(ierr); 197 ierr = PetscMalloc(nz*sizeof(PetscScalar), &val);CHKERRQ(ierr); 198 nz = 0; 199 for(i=0; i<M; i++) { 200 rnz = ai[i+1] - adiag[i]; 201 ajj = aj + adiag[i]; 202 for(j=0; j<rnz; j++) { 203 row[nz] = i+shift; col[nz] = ajj[j] + shift; val[nz++] = v1[j]; 204 } 205 } 206 *r = row; *c = col; *v = val; 207 } else { 208 nz = 0; val = *v; 209 for(i=0; i <M; i++) { 210 rnz = ai[i+1] - adiag[i]; 211 ajj = aj + adiag[i]; 212 v1 = av + adiag[i]; 213 for(j=0; j<rnz; j++) { 214 val[nz++] = v1[j]; 215 } 216 } 217 } 218 PetscFunctionReturn(0); 219 } 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "MatConvertToTriples_mpisbaij_mpisbaij" 223 PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 224 { 225 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 226 PetscErrorCode ierr; 227 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 228 PetscInt *row,*col; 229 const PetscScalar *av, *bv,*v1,*v2; 230 PetscScalar *val; 231 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 232 Mat_SeqSBAIJ *aa=(Mat_SeqSBAIJ*)(mat->A)->data; 233 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 234 235 PetscFunctionBegin; 236 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 237 garray = mat->garray; 238 av=aa->a; bv=bb->a; 239 240 if (!valOnly){ 241 nz = aa->nz + bb->nz; 242 *nnz = nz; 243 ierr = PetscMalloc(nz*sizeof(PetscInt) ,&row);CHKERRQ(ierr); 244 ierr = PetscMalloc(nz*sizeof(PetscInt),&col);CHKERRQ(ierr); 245 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 246 *r = row; *c = col; *v = val; 247 } else { 248 row = *r; col = *c; val = *v; 249 } 250 251 jj = 0; irow = rstart; 252 for ( i=0; i<m; i++ ) { 253 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 254 countA = ai[i+1] - ai[i]; 255 countB = bi[i+1] - bi[i]; 256 bjj = bj + bi[i]; 257 v1 = av + ai[i]; 258 v2 = bv + bi[i]; 259 260 /* A-part */ 261 for (j=0; j<countA; j++){ 262 if (!valOnly){ 263 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 264 } 265 val[jj++] = v1[j]; 266 } 267 268 /* B-part */ 269 for(j=0; j < countB; j++){ 270 if (!valOnly){ 271 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 272 } 273 val[jj++] = v2[j]; 274 } 275 irow++; 276 } 277 PetscFunctionReturn(0); 278 } 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpiaij" 282 PetscErrorCode MatConvertToTriples_mpiaij_mpiaij(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 283 { 284 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 285 PetscErrorCode ierr; 286 PetscInt rstart,nz,i,j,jj,irow,countA,countB; 287 PetscInt *row,*col; 288 const PetscScalar *av, *bv,*v1,*v2; 289 PetscScalar *val; 290 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 291 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 292 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 293 294 PetscFunctionBegin; 295 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart= A->rmap->rstart; 296 garray = mat->garray; 297 av=aa->a; bv=bb->a; 298 299 if (!valOnly){ 300 nz = aa->nz + bb->nz; 301 *nnz = nz; 302 ierr = PetscMalloc(nz*sizeof(PetscInt) ,&row);CHKERRQ(ierr); 303 ierr = PetscMalloc(nz*sizeof(PetscInt),&col);CHKERRQ(ierr); 304 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 305 *r = row; *c = col; *v = val; 306 } else { 307 row = *r; col = *c; val = *v; 308 } 309 310 jj = 0; irow = rstart; 311 for ( i=0; i<m; i++ ) { 312 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 313 countA = ai[i+1] - ai[i]; 314 countB = bi[i+1] - bi[i]; 315 bjj = bj + bi[i]; 316 v1 = av + ai[i]; 317 v2 = bv + bi[i]; 318 319 /* A-part */ 320 for (j=0; j<countA; j++){ 321 if (!valOnly){ 322 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 323 } 324 val[jj++] = v1[j]; 325 } 326 327 /* B-part */ 328 for(j=0; j < countB; j++){ 329 if (!valOnly){ 330 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 331 } 332 val[jj++] = v2[j]; 333 } 334 irow++; 335 } 336 PetscFunctionReturn(0); 337 } 338 339 #undef __FUNCT__ 340 #define __FUNCT__ "MatConvertToTriples_mpibaij_mpiaij" 341 PetscErrorCode MatConvertToTriples_mpibaij_mpiaij(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 342 { 343 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data; 344 Mat_SeqBAIJ *aa=(Mat_SeqBAIJ*)(mat->A)->data; 345 Mat_SeqBAIJ *bb=(Mat_SeqBAIJ*)(mat->B)->data; 346 const PetscInt *ai = aa->i, *bi = bb->i, *aj = aa->j, *bj = bb->j,*ajj, *bjj; 347 const PetscInt *garray = mat->garray,mbs=mat->mbs,rstartbs=mat->rstartbs; 348 const PetscInt bs = A->rmap->bs,bs2=mat->bs2; 349 PetscErrorCode ierr; 350 PetscInt nz,i,j,k,n,jj,irow,countA,countB,idx; 351 PetscInt *row,*col; 352 const PetscScalar *av=aa->a, *bv=bb->a,*v1,*v2; 353 PetscScalar *val; 354 355 PetscFunctionBegin; 356 357 if (!valOnly){ 358 nz = bs2*(aa->nz + bb->nz); 359 *nnz = nz; 360 ierr = PetscMalloc(nz*sizeof(PetscInt) ,&row);CHKERRQ(ierr); 361 ierr = PetscMalloc(nz*sizeof(PetscInt),&col);CHKERRQ(ierr); 362 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 363 *r = row; *c = col; *v = val; 364 } else { 365 row = *r; col = *c; val = *v; 366 } 367 368 jj = 0; irow = rstartbs; 369 for ( i=0; i<mbs; i++ ) { 370 countA = ai[i+1] - ai[i]; 371 countB = bi[i+1] - bi[i]; 372 ajj = aj + ai[i]; 373 bjj = bj + bi[i]; 374 v1 = av + bs2*ai[i]; 375 v2 = bv + bs2*bi[i]; 376 377 idx = 0; 378 /* A-part */ 379 for (k=0; k<countA; k++){ 380 for (j=0; j<bs; j++) { 381 for (n=0; n<bs; n++) { 382 if (!valOnly){ 383 row[jj] = bs*irow + n + shift; 384 col[jj] = bs*(rstartbs + ajj[k]) + j + shift; 385 } 386 val[jj++] = v1[idx++]; 387 } 388 } 389 } 390 391 idx = 0; 392 /* B-part */ 393 for(k=0; k<countB; k++){ 394 for (j=0; j<bs; j++) { 395 for (n=0; n<bs; n++) { 396 if (!valOnly){ 397 row[jj] = bs*irow + n + shift; 398 col[jj] = bs*(garray[bjj[k]]) + j + shift; 399 } 400 val[jj++] = bv[idx++]; 401 } 402 } 403 } 404 irow++; 405 } 406 PetscFunctionReturn(0); 407 } 408 409 #undef __FUNCT__ 410 #define __FUNCT__ "MatConvertToTriples_mpiaij_mpisbaij" 411 PetscErrorCode MatConvertToTriples_mpiaij_mpisbaij(Mat A,int shift,PetscTruth valOnly,int *nnz,int **r, int **c, PetscScalar **v) 412 { 413 const PetscInt *ai, *aj,*adiag, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 414 PetscErrorCode ierr; 415 PetscInt rstart,nz,nza,nzb_low,i,j,jj,irow,countA,countB; 416 PetscInt *row,*col; 417 const PetscScalar *av, *bv,*v1,*v2; 418 PetscScalar *val; 419 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 420 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(mat->A)->data; 421 Mat_SeqAIJ *bb=(Mat_SeqAIJ*)(mat->B)->data; 422 423 PetscFunctionBegin; 424 ai=aa->i; aj=aa->j; adiag=aa->diag; 425 bi=bb->i; bj=bb->j; garray = mat->garray; 426 av=aa->a; bv=bb->a; 427 rstart = A->rmap->rstart; 428 429 if (!valOnly){ 430 nza = 0;nzb_low = 0; 431 for(i=0; i<m; i++){ 432 nza = nza + (ai[i+1] - adiag[i]); 433 countB = bi[i+1] - bi[i]; 434 bjj = bj + bi[i]; 435 436 j = 0; 437 while(garray[bjj[j]] < rstart) { 438 if(j == countB) break; 439 j++;nzb_low++; 440 } 441 } 442 /* Total nz = nz for the upper triangular A part + nz for the 2nd B part */ 443 nz = nza + (bb->nz - nzb_low); 444 *nnz = nz; 445 ierr = PetscMalloc(nz*sizeof(PetscInt) ,&row);CHKERRQ(ierr); 446 ierr = PetscMalloc(nz*sizeof(PetscInt),&col);CHKERRQ(ierr); 447 ierr = PetscMalloc(nz*sizeof(PetscScalar),&val);CHKERRQ(ierr); 448 *r = row; *c = col; *v = val; 449 } else { 450 row = *r; col = *c; val = *v; 451 } 452 453 jj = 0; irow = rstart; 454 for ( i=0; i<m; i++ ) { 455 ajj = aj + adiag[i]; /* ptr to the beginning of the diagonal of this row */ 456 v1 = av + adiag[i]; 457 countA = ai[i+1] - adiag[i]; 458 countB = bi[i+1] - bi[i]; 459 bjj = bj + bi[i]; 460 v2 = bv + bi[i]; 461 462 /* A-part */ 463 for (j=0; j<countA; j++){ 464 if (!valOnly){ 465 row[jj] = irow + shift; col[jj] = rstart + ajj[j] + shift; 466 } 467 val[jj++] = v1[j]; 468 } 469 470 /* B-part */ 471 for(j=0; j < countB; j++){ 472 if (garray[bjj[j]] > rstart) { 473 if (!valOnly){ 474 row[jj] = irow + shift; col[jj] = garray[bjj[j]] + shift; 475 } 476 val[jj++] = v2[j]; 477 } 478 } 479 irow++; 480 } 481 PetscFunctionReturn(0); 482 } 483 484 #undef __FUNCT__ 485 #define __FUNCT__ "MatDestroy_MUMPS" 486 PetscErrorCode MatDestroy_MUMPS(Mat A) 487 { 488 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 489 PetscErrorCode ierr; 490 PetscMPIInt size=lu->size; 491 PetscTruth isSeqBAIJ,isMPIBAIJ; 492 493 PetscFunctionBegin; 494 ierr = PetscTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 495 ierr = PetscTypeCompare((PetscObject)A,MATMPIBAIJ,&isMPIBAIJ);CHKERRQ(ierr); 496 if (lu->CleanUpMUMPS) { 497 /* Terminate instance, deallocate memories */ 498 if (size > 1){ 499 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 500 ierr = VecScatterDestroy(lu->scat_rhs);CHKERRQ(ierr); 501 ierr = VecDestroy(lu->b_seq);CHKERRQ(ierr); 502 if (lu->nSolve && lu->scat_sol){ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr);} 503 if (lu->nSolve && lu->x_seq){ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr);} 504 ierr = PetscFree(lu->val);CHKERRQ(ierr); 505 } 506 if( size == 1 && A->factortype == MAT_FACTOR_CHOLESKY && lu->isAIJ) { 507 ierr = PetscFree(lu->val);CHKERRQ(ierr); 508 } 509 if(isSeqBAIJ || isMPIBAIJ) { 510 ierr = PetscFree(lu->val);CHKERRQ(ierr); 511 } 512 lu->id.job=JOB_END; 513 #if defined(PETSC_USE_COMPLEX) 514 zmumps_c(&lu->id); 515 #else 516 dmumps_c(&lu->id); 517 #endif 518 ierr = PetscFree(lu->irn);CHKERRQ(ierr); 519 ierr = PetscFree(lu->jcn);CHKERRQ(ierr); 520 ierr = MPI_Comm_free(&(lu->comm_mumps));CHKERRQ(ierr); 521 } 522 /* clear composed functions */ 523 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatFactorGetSolverPackage_C","",PETSC_NULL);CHKERRQ(ierr); 524 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatMumpsSetIcntl_C","",PETSC_NULL);CHKERRQ(ierr); 525 ierr = (lu->MatDestroy)(A);CHKERRQ(ierr); 526 PetscFunctionReturn(0); 527 } 528 529 #undef __FUNCT__ 530 #define __FUNCT__ "MatSolve_MUMPS" 531 PetscErrorCode MatSolve_MUMPS(Mat A,Vec b,Vec x) 532 { 533 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 534 PetscScalar *array; 535 Vec b_seq; 536 IS is_iden,is_petsc; 537 PetscErrorCode ierr; 538 PetscInt i; 539 540 PetscFunctionBegin; 541 lu->id.nrhs = 1; 542 b_seq = lu->b_seq; 543 if (lu->size > 1){ 544 /* MUMPS only supports centralized rhs. Scatter b into a seqential rhs vector */ 545 ierr = VecScatterBegin(lu->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 546 ierr = VecScatterEnd(lu->scat_rhs,b,b_seq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 547 if (!lu->myid) {ierr = VecGetArray(b_seq,&array);CHKERRQ(ierr);} 548 } else { /* size == 1 */ 549 ierr = VecCopy(b,x);CHKERRQ(ierr); 550 ierr = VecGetArray(x,&array);CHKERRQ(ierr); 551 } 552 if (!lu->myid) { /* define rhs on the host */ 553 lu->id.nrhs = 1; 554 #if defined(PETSC_USE_COMPLEX) 555 lu->id.rhs = (mumps_double_complex*)array; 556 #else 557 lu->id.rhs = array; 558 #endif 559 } 560 561 /* solve phase */ 562 /*-------------*/ 563 lu->id.job = 3; 564 #if defined(PETSC_USE_COMPLEX) 565 zmumps_c(&lu->id); 566 #else 567 dmumps_c(&lu->id); 568 #endif 569 if (lu->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in solve phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 570 571 if (lu->size > 1) { /* convert mumps distributed solution to petsc mpi x */ 572 if (!lu->nSolve){ /* create scatter scat_sol */ 573 ierr = ISCreateStride(PETSC_COMM_SELF,lu->id.lsol_loc,0,1,&is_iden);CHKERRQ(ierr); /* from */ 574 for (i=0; i<lu->id.lsol_loc; i++){ 575 lu->id.isol_loc[i] -= 1; /* change Fortran style to C style */ 576 } 577 ierr = ISCreateGeneral(PETSC_COMM_SELF,lu->id.lsol_loc,lu->id.isol_loc,&is_petsc);CHKERRQ(ierr); /* to */ 578 ierr = VecScatterCreate(lu->x_seq,is_iden,x,is_petsc,&lu->scat_sol);CHKERRQ(ierr); 579 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 580 ierr = ISDestroy(is_petsc);CHKERRQ(ierr); 581 } 582 ierr = VecScatterBegin(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 583 ierr = VecScatterEnd(lu->scat_sol,lu->x_seq,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 584 } 585 lu->nSolve++; 586 PetscFunctionReturn(0); 587 } 588 589 #if !defined(PETSC_USE_COMPLEX) 590 /* 591 input: 592 F: numeric factor 593 output: 594 nneg: total number of negative pivots 595 nzero: 0 596 npos: (global dimension of F) - nneg 597 */ 598 599 #undef __FUNCT__ 600 #define __FUNCT__ "MatGetInertia_SBAIJMUMPS" 601 PetscErrorCode MatGetInertia_SBAIJMUMPS(Mat F,int *nneg,int *nzero,int *npos) 602 { 603 Mat_MUMPS *lu =(Mat_MUMPS*)F->spptr; 604 PetscErrorCode ierr; 605 PetscMPIInt size; 606 607 PetscFunctionBegin; 608 ierr = MPI_Comm_size(((PetscObject)F)->comm,&size);CHKERRQ(ierr); 609 /* 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 */ 610 if (size > 1 && lu->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",lu->id.INFOG(13)); 611 if (nneg){ 612 if (!lu->myid){ 613 *nneg = lu->id.INFOG(12); 614 } 615 ierr = MPI_Bcast(nneg,1,MPI_INT,0,lu->comm_mumps);CHKERRQ(ierr); 616 } 617 if (nzero) *nzero = 0; 618 if (npos) *npos = F->rmap->N - (*nneg); 619 PetscFunctionReturn(0); 620 } 621 #endif /* !defined(PETSC_USE_COMPLEX) */ 622 623 #undef __FUNCT__ 624 #define __FUNCT__ "MatFactorNumeric_MUMPS" 625 PetscErrorCode MatFactorNumeric_MUMPS(Mat F,Mat A,const MatFactorInfo *info) 626 { 627 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 628 PetscErrorCode ierr; 629 PetscTruth valOnly; 630 Mat F_diag; 631 PetscTruth isMPIAIJ; 632 633 PetscFunctionBegin; 634 valOnly = PETSC_TRUE; 635 ierr = (*lu->ConvertToTriples)(A, 1, valOnly, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 636 637 /* numerical factorization phase */ 638 /*-------------------------------*/ 639 lu->id.job = 2; 640 if(!lu->id.ICNTL(18)) { 641 if (!lu->myid) { 642 #if defined(PETSC_USE_COMPLEX) 643 lu->id.a = (mumps_double_complex*)lu->val; 644 #else 645 lu->id.a = lu->val; 646 #endif 647 } 648 } else { 649 #if defined(PETSC_USE_COMPLEX) 650 lu->id.a_loc = (mumps_double_complex*)lu->val; 651 #else 652 lu->id.a_loc = lu->val; 653 #endif 654 } 655 #if defined(PETSC_USE_COMPLEX) 656 zmumps_c(&lu->id); 657 #else 658 dmumps_c(&lu->id); 659 #endif 660 if (lu->id.INFOG(1) < 0) { 661 if (lu->id.INFO(1) == -13) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: Cannot allocate required memory %d megabytes\n",lu->id.INFO(2)); 662 else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in numerical factorization phase: INFO(1)=%d, INFO(2)=%d\n",lu->id.INFO(1),lu->id.INFO(2)); 663 } 664 if (!lu->myid && lu->id.ICNTL(16) > 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB," lu->id.ICNTL(16):=%d\n",lu->id.INFOG(16)); 665 666 if (lu->size > 1){ 667 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 668 if(isMPIAIJ) { 669 F_diag = ((Mat_MPIAIJ *)(F)->data)->A; 670 } else { 671 F_diag = ((Mat_MPISBAIJ *)(F)->data)->A; 672 } 673 F_diag->assembled = PETSC_TRUE; 674 if (lu->nSolve){ 675 ierr = VecScatterDestroy(lu->scat_sol);CHKERRQ(ierr); 676 ierr = PetscFree2(lu->id.sol_loc,lu->id.isol_loc);CHKERRQ(ierr); 677 ierr = VecDestroy(lu->x_seq);CHKERRQ(ierr); 678 } 679 } 680 (F)->assembled = PETSC_TRUE; 681 lu->matstruc = SAME_NONZERO_PATTERN; 682 lu->CleanUpMUMPS = PETSC_TRUE; 683 lu->nSolve = 0; 684 685 if (lu->size > 1){ 686 /* distributed solution */ 687 lu->id.ICNTL(21) = 1; 688 if (!lu->nSolve){ 689 /* Create x_seq=sol_loc for repeated use */ 690 PetscInt lsol_loc; 691 PetscScalar *sol_loc; 692 lsol_loc = lu->id.INFO(23); /* length of sol_loc */ 693 ierr = PetscMalloc2(lsol_loc,PetscScalar,&sol_loc,lsol_loc,PetscInt,&lu->id.isol_loc);CHKERRQ(ierr); 694 lu->id.lsol_loc = lsol_loc; 695 #if defined(PETSC_USE_COMPLEX) 696 lu->id.sol_loc = (mumps_double_complex*)sol_loc; 697 #else 698 lu->id.sol_loc = sol_loc; 699 #endif 700 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,lsol_loc,sol_loc,&lu->x_seq);CHKERRQ(ierr); 701 } 702 } 703 PetscFunctionReturn(0); 704 } 705 706 #undef __FUNCT__ 707 #define __FUNCT__ "PetscSetMUMPSOptions" 708 PetscErrorCode PetscSetMUMPSOptions(Mat F, Mat A) 709 { 710 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 711 PetscErrorCode ierr; 712 PetscInt icntl; 713 PetscTruth flg; 714 715 PetscFunctionBegin; 716 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 717 ierr = PetscOptionsTruth("-mat_mumps_view","View MUMPS parameters","None",lu->mumpsview,&lu->mumpsview,PETSC_NULL);CHKERRQ(ierr); 718 if (lu->size == 1){ 719 lu->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 720 } else { 721 lu->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 722 } 723 724 icntl=-1; 725 lu->id.ICNTL(4) = 0; /* level of printing; overwrite mumps default ICNTL(4)=2 */ 726 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",lu->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 727 if ((flg && icntl > 0) || PetscLogPrintInfo) { 728 lu->id.ICNTL(4)=icntl; /* and use mumps default icntl(i), i=1,2,3 */ 729 } else { /* no output */ 730 lu->id.ICNTL(1) = 0; /* error message, default= 6 */ 731 lu->id.ICNTL(2) = 0; /* output stream for diagnostic printing, statistics, and warning. default=0 */ 732 lu->id.ICNTL(3) = 0; /* output stream for global information, default=6 */ 733 } 734 ierr = PetscOptionsInt("-mat_mumps_icntl_6","ICNTL(6): column permutation and/or scaling to get a zero-free diagonal (0 to 7)","None",lu->id.ICNTL(6),&lu->id.ICNTL(6),PETSC_NULL);CHKERRQ(ierr); 735 icntl=-1; 736 ierr = PetscOptionsInt("-mat_mumps_icntl_7","ICNTL(7): matrix ordering (0 to 7)","None",lu->id.ICNTL(7),&icntl,&flg);CHKERRQ(ierr); 737 if (flg) { 738 if (icntl== 1){ 739 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"); 740 } else { 741 lu->id.ICNTL(7) = icntl; 742 } 743 } 744 ierr = PetscOptionsInt("-mat_mumps_icntl_8","ICNTL(8): scaling strategy (-2 to 7 or 77)","None",lu->id.ICNTL(8),&lu->id.ICNTL(8),PETSC_NULL);CHKERRQ(ierr); 745 ierr = PetscOptionsInt("-mat_mumps_icntl_9","ICNTL(9): A or A^T x=b to be solved. 1: A; otherwise: A^T","None",lu->id.ICNTL(9),&lu->id.ICNTL(9),PETSC_NULL);CHKERRQ(ierr); 746 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",lu->id.ICNTL(10),&lu->id.ICNTL(10),PETSC_NULL);CHKERRQ(ierr); 747 ierr = PetscOptionsInt("-mat_mumps_icntl_11","ICNTL(11): statistics related to the linear system solved (via -ksp_view)","None",lu->id.ICNTL(11),&lu->id.ICNTL(11),PETSC_NULL);CHKERRQ(ierr); 748 ierr = PetscOptionsInt("-mat_mumps_icntl_12","ICNTL(12): efficiency control: defines the ordering strategy with scaling constraints (0 to 3","None",lu->id.ICNTL(12),&lu->id.ICNTL(12),PETSC_NULL);CHKERRQ(ierr); 749 ierr = PetscOptionsInt("-mat_mumps_icntl_13","ICNTL(13): efficiency control: with or without ScaLAPACK","None",lu->id.ICNTL(13),&lu->id.ICNTL(13),PETSC_NULL);CHKERRQ(ierr); 750 ierr = PetscOptionsInt("-mat_mumps_icntl_14","ICNTL(14): percentage of estimated workspace increase","None",lu->id.ICNTL(14),&lu->id.ICNTL(14),PETSC_NULL);CHKERRQ(ierr); 751 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): Schur complement","None",lu->id.ICNTL(19),&lu->id.ICNTL(19),PETSC_NULL);CHKERRQ(ierr); 752 753 ierr = PetscOptionsInt("-mat_mumps_icntl_22","ICNTL(22): in-core/out-of-core facility (0 or 1)","None",lu->id.ICNTL(22),&lu->id.ICNTL(22),PETSC_NULL);CHKERRQ(ierr); 754 ierr = PetscOptionsInt("-mat_mumps_icntl_23","ICNTL(23): max size of the working memory (MB) that can allocate per processor","None",lu->id.ICNTL(23),&lu->id.ICNTL(23),PETSC_NULL);CHKERRQ(ierr); 755 ierr = PetscOptionsInt("-mat_mumps_icntl_24","ICNTL(24): detection of null pivot rows (0 or 1)","None",lu->id.ICNTL(24),&lu->id.ICNTL(24),PETSC_NULL);CHKERRQ(ierr); 756 ierr = PetscOptionsInt("-mat_mumps_icntl_25","ICNTL(25): computation of a null space basis","None",lu->id.ICNTL(25),&lu->id.ICNTL(25),PETSC_NULL);CHKERRQ(ierr); 757 ierr = PetscOptionsInt("-mat_mumps_icntl_26","ICNTL(26): Schur options for right-hand side or solution vector","None",lu->id.ICNTL(26),&lu->id.ICNTL(26),PETSC_NULL);CHKERRQ(ierr); 758 ierr = PetscOptionsInt("-mat_mumps_icntl_27","ICNTL(27): experimental parameter","None",lu->id.ICNTL(27),&lu->id.ICNTL(27),PETSC_NULL);CHKERRQ(ierr); 759 760 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",lu->id.CNTL(1),&lu->id.CNTL(1),PETSC_NULL);CHKERRQ(ierr); 761 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",lu->id.CNTL(2),&lu->id.CNTL(2),PETSC_NULL);CHKERRQ(ierr); 762 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",lu->id.CNTL(3),&lu->id.CNTL(3),PETSC_NULL);CHKERRQ(ierr); 763 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",lu->id.CNTL(4),&lu->id.CNTL(4),PETSC_NULL);CHKERRQ(ierr); 764 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",lu->id.CNTL(5),&lu->id.CNTL(5),PETSC_NULL);CHKERRQ(ierr); 765 PetscOptionsEnd(); 766 PetscFunctionReturn(0); 767 } 768 769 #undef __FUNCT__ 770 #define __FUNCT__ "PetscInitializeMUMPS" 771 PetscErrorCode PetscInitializeMUMPS(Mat F) 772 { 773 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 774 PetscErrorCode ierr; 775 PetscInt icntl; 776 PetscTruth flg; 777 778 PetscFunctionBegin; 779 lu->id.job = JOB_INIT; 780 lu->id.par=1; /* host participates factorizaton and solve */ 781 lu->id.sym=lu->sym; 782 if (lu->sym == 2){ 783 ierr = PetscOptionsInt("-mat_mumps_sym","SYM: (1,2)","None",lu->id.sym,&icntl,&flg);CHKERRQ(ierr); 784 if (flg && icntl == 1) lu->id.sym=icntl; /* matrix is spd */ 785 } 786 #if defined(PETSC_USE_COMPLEX) 787 zmumps_c(&lu->id); 788 #else 789 dmumps_c(&lu->id); 790 #endif 791 PetscFunctionReturn(0); 792 } 793 794 /* Note the Petsc r and c permutations are ignored */ 795 #undef __FUNCT__ 796 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 797 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 798 { 799 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 800 PetscErrorCode ierr; 801 PetscTruth valOnly; 802 Vec b; 803 IS is_iden; 804 const PetscInt M = A->rmap->N; 805 806 PetscFunctionBegin; 807 lu->sym = 0; 808 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 809 810 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 811 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 812 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 813 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 814 815 /* Initialize a MUMPS instance */ 816 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 817 /* Set MUMPS options */ 818 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 819 820 valOnly = PETSC_FALSE; 821 ierr = (*lu->ConvertToTriples)(A, 1, valOnly, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 822 823 /* analysis phase */ 824 /*----------------*/ 825 lu->id.job = 1; 826 lu->id.n = M; 827 switch (lu->id.ICNTL(18)){ 828 case 0: /* centralized assembled matrix input */ 829 if (!lu->myid) { 830 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 831 if (lu->id.ICNTL(6)>1){ 832 #if defined(PETSC_USE_COMPLEX) 833 lu->id.a = (mumps_double_complex*)lu->val; 834 #else 835 lu->id.a = lu->val; 836 #endif 837 } 838 } 839 break; 840 case 3: /* distributed assembled matrix input (size>1) */ 841 lu->id.nz_loc = lu->nz; 842 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 843 if (lu->id.ICNTL(6)>1) { 844 #if defined(PETSC_USE_COMPLEX) 845 lu->id.a_loc = (mumps_double_complex*)lu->val; 846 #else 847 lu->id.a_loc = lu->val; 848 #endif 849 } 850 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 851 if (!lu->myid){ 852 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 853 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 854 } else { 855 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 856 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 857 } 858 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 859 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 860 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 861 862 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 863 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 864 ierr = VecDestroy(b);CHKERRQ(ierr); 865 break; 866 } 867 #if defined(PETSC_USE_COMPLEX) 868 zmumps_c(&lu->id); 869 #else 870 dmumps_c(&lu->id); 871 #endif 872 if (lu->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 873 874 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 875 F->ops->solve = MatSolve_MUMPS; 876 PetscFunctionReturn(0); 877 } 878 879 /* Note the Petsc r and c permutations are ignored */ 880 #undef __FUNCT__ 881 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 882 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 883 { 884 885 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 886 PetscErrorCode ierr; 887 PetscTruth valOnly; 888 Vec b; 889 IS is_iden; 890 const PetscInt M = A->rmap->N; 891 892 PetscFunctionBegin; 893 lu->sym = 0; 894 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 895 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 896 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 897 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 898 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 899 900 /* Initialize a MUMPS instance */ 901 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 902 /* Set MUMPS options */ 903 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 904 905 valOnly = PETSC_FALSE; 906 ierr = (*lu->ConvertToTriples)(A, 1, valOnly, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 907 908 /* analysis phase */ 909 /*----------------*/ 910 lu->id.job = 1; 911 lu->id.n = M; 912 switch (lu->id.ICNTL(18)){ 913 case 0: /* centralized assembled matrix input */ 914 if (!lu->myid) { 915 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 916 if (lu->id.ICNTL(6)>1){ 917 #if defined(PETSC_USE_COMPLEX) 918 lu->id.a = (mumps_double_complex*)lu->val; 919 #else 920 lu->id.a = lu->val; 921 #endif 922 } 923 } 924 break; 925 case 3: /* distributed assembled matrix input (size>1) */ 926 lu->id.nz_loc = lu->nz; 927 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 928 if (lu->id.ICNTL(6)>1) { 929 #if defined(PETSC_USE_COMPLEX) 930 lu->id.a_loc = (mumps_double_complex*)lu->val; 931 #else 932 lu->id.a_loc = lu->val; 933 #endif 934 } 935 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 936 if (!lu->myid){ 937 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 938 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 939 } else { 940 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 941 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 942 } 943 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 944 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 945 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 946 947 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 948 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 949 ierr = VecDestroy(b);CHKERRQ(ierr); 950 break; 951 } 952 #if defined(PETSC_USE_COMPLEX) 953 zmumps_c(&lu->id); 954 #else 955 dmumps_c(&lu->id); 956 #endif 957 if (lu->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 958 959 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 960 F->ops->solve = MatSolve_MUMPS; 961 PetscFunctionReturn(0); 962 } 963 964 /* Note the Petsc r permutation is ignored */ 965 #undef __FUNCT__ 966 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 967 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 968 { 969 Mat_MUMPS *lu = (Mat_MUMPS*)F->spptr; 970 PetscErrorCode ierr; 971 PetscTruth valOnly; 972 Vec b; 973 IS is_iden; 974 const PetscInt M = A->rmap->N; 975 976 PetscFunctionBegin; 977 lu->sym = 2; 978 lu->matstruc = DIFFERENT_NONZERO_PATTERN; 979 ierr = MPI_Comm_rank(((PetscObject)A)->comm, &lu->myid); 980 ierr = MPI_Comm_size(((PetscObject)A)->comm,&lu->size);CHKERRQ(ierr); 981 ierr = MPI_Comm_dup(((PetscObject)A)->comm,&(lu->comm_mumps));CHKERRQ(ierr); 982 lu->id.comm_fortran = MPI_Comm_c2f(lu->comm_mumps); 983 984 /* Initialize a MUMPS instance */ 985 ierr = PetscInitializeMUMPS(F);CHKERRQ(ierr); 986 /* Set MUMPS options */ 987 ierr = PetscSetMUMPSOptions(F,A);CHKERRQ(ierr); 988 989 valOnly = PETSC_FALSE; 990 ierr = (*lu->ConvertToTriples)(A, 1 , valOnly, &lu->nz, &lu->irn, &lu->jcn, &lu->val);CHKERRQ(ierr); 991 992 /* analysis phase */ 993 /*----------------*/ 994 lu->id.job = 1; 995 lu->id.n = M; 996 switch (lu->id.ICNTL(18)){ 997 case 0: /* centralized assembled matrix input */ 998 if (!lu->myid) { 999 lu->id.nz =lu->nz; lu->id.irn=lu->irn; lu->id.jcn=lu->jcn; 1000 if (lu->id.ICNTL(6)>1){ 1001 #if defined(PETSC_USE_COMPLEX) 1002 lu->id.a = (mumps_double_complex*)lu->val; 1003 #else 1004 lu->id.a = lu->val; 1005 #endif 1006 } 1007 } 1008 break; 1009 case 3: /* distributed assembled matrix input (size>1) */ 1010 lu->id.nz_loc = lu->nz; 1011 lu->id.irn_loc=lu->irn; lu->id.jcn_loc=lu->jcn; 1012 if (lu->id.ICNTL(6)>1) { 1013 #if defined(PETSC_USE_COMPLEX) 1014 lu->id.a_loc = (mumps_double_complex*)lu->val; 1015 #else 1016 lu->id.a_loc = lu->val; 1017 #endif 1018 } 1019 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1020 if (!lu->myid){ 1021 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&lu->b_seq);CHKERRQ(ierr); 1022 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1023 } else { 1024 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&lu->b_seq);CHKERRQ(ierr); 1025 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1026 } 1027 ierr = VecCreate(((PetscObject)A)->comm,&b);CHKERRQ(ierr); 1028 ierr = VecSetSizes(b,A->rmap->n,PETSC_DECIDE);CHKERRQ(ierr); 1029 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 1030 1031 ierr = VecScatterCreate(b,is_iden,lu->b_seq,is_iden,&lu->scat_rhs);CHKERRQ(ierr); 1032 ierr = ISDestroy(is_iden);CHKERRQ(ierr); 1033 ierr = VecDestroy(b);CHKERRQ(ierr); 1034 break; 1035 } 1036 #if defined(PETSC_USE_COMPLEX) 1037 zmumps_c(&lu->id); 1038 #else 1039 dmumps_c(&lu->id); 1040 #endif 1041 if (lu->id.INFOG(1) < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MUMPS in analysis phase: INFOG(1)=%d\n",lu->id.INFOG(1)); 1042 1043 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1044 F->ops->solve = MatSolve_MUMPS; 1045 #if !defined(PETSC_USE_COMPLEX) 1046 (F)->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1047 #endif 1048 PetscFunctionReturn(0); 1049 } 1050 1051 #undef __FUNCT__ 1052 #define __FUNCT__ "MatFactorInfo_MUMPS" 1053 PetscErrorCode MatFactorInfo_MUMPS(Mat A,PetscViewer viewer) 1054 { 1055 Mat_MUMPS *lu=(Mat_MUMPS*)A->spptr; 1056 PetscErrorCode ierr; 1057 1058 PetscFunctionBegin; 1059 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",lu->id.ICNTL(1));CHKERRQ(ierr); 1060 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg):%d \n",lu->id.ICNTL(2));CHKERRQ(ierr); 1061 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",lu->id.ICNTL(3));CHKERRQ(ierr); 1062 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",lu->id.ICNTL(4));CHKERRQ(ierr); 1063 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",lu->id.ICNTL(5));CHKERRQ(ierr); 1064 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",lu->id.ICNTL(6));CHKERRQ(ierr); 1065 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (matrix ordering): %d \n",lu->id.ICNTL(7));CHKERRQ(ierr); 1066 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",lu->id.ICNTL(8));CHKERRQ(ierr); 1067 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(9) (A/A^T x=b is solved): %d \n",lu->id.ICNTL(9));CHKERRQ(ierr); 1068 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",lu->id.ICNTL(10));CHKERRQ(ierr); 1069 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",lu->id.ICNTL(11));CHKERRQ(ierr); 1070 if (!lu->myid && lu->id.ICNTL(11)>0) { 1071 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(4) (inf norm of input mat): %g\n",lu->id.RINFOG(4));CHKERRQ(ierr); 1072 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(5) (inf norm of solution): %g\n",lu->id.RINFOG(5));CHKERRQ(ierr); 1073 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(6) (inf norm of residual): %g\n",lu->id.RINFOG(6));CHKERRQ(ierr); 1074 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",lu->id.RINFOG(7),lu->id.RINFOG(8));CHKERRQ(ierr); 1075 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(9) (error estimate): %g \n",lu->id.RINFOG(9));CHKERRQ(ierr); 1076 ierr = PetscPrintf(PETSC_COMM_SELF," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",lu->id.RINFOG(10),lu->id.RINFOG(11));CHKERRQ(ierr); 1077 1078 } 1079 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",lu->id.ICNTL(12));CHKERRQ(ierr); 1080 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",lu->id.ICNTL(13));CHKERRQ(ierr); 1081 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",lu->id.ICNTL(14));CHKERRQ(ierr); 1082 /* ICNTL(15-17) not used */ 1083 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",lu->id.ICNTL(18));CHKERRQ(ierr); 1084 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",lu->id.ICNTL(19));CHKERRQ(ierr); 1085 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",lu->id.ICNTL(20));CHKERRQ(ierr); 1086 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",lu->id.ICNTL(21));CHKERRQ(ierr); 1087 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",lu->id.ICNTL(22));CHKERRQ(ierr); 1088 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",lu->id.ICNTL(23));CHKERRQ(ierr); 1089 1090 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",lu->id.ICNTL(24));CHKERRQ(ierr); 1091 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",lu->id.ICNTL(25));CHKERRQ(ierr); 1092 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",lu->id.ICNTL(26));CHKERRQ(ierr); 1093 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",lu->id.ICNTL(27));CHKERRQ(ierr); 1094 1095 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",lu->id.CNTL(1));CHKERRQ(ierr); 1096 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",lu->id.CNTL(2));CHKERRQ(ierr); 1097 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",lu->id.CNTL(3));CHKERRQ(ierr); 1098 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",lu->id.CNTL(4));CHKERRQ(ierr); 1099 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",lu->id.CNTL(5));CHKERRQ(ierr); 1100 1101 /* infomation local to each processor */ 1102 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr);} 1103 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(1));CHKERRQ(ierr); 1104 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 1105 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr);} 1106 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(2));CHKERRQ(ierr); 1107 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 1108 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr);} 1109 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %g \n",lu->myid,lu->id.RINFO(3));CHKERRQ(ierr); 1110 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 1111 1112 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr);} 1113 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(15));CHKERRQ(ierr); 1114 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 1115 1116 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr);} 1117 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(16));CHKERRQ(ierr); 1118 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 1119 1120 if (!lu->myid) {ierr = PetscPrintf(PETSC_COMM_SELF, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr);} 1121 ierr = PetscSynchronizedPrintf(((PetscObject)A)->comm," [%d] %d \n",lu->myid,lu->id.INFO(23));CHKERRQ(ierr); 1122 ierr = PetscSynchronizedFlush(((PetscObject)A)->comm); 1123 1124 if (!lu->myid){ /* information from the host */ 1125 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",lu->id.RINFOG(1));CHKERRQ(ierr); 1126 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",lu->id.RINFOG(2));CHKERRQ(ierr); 1127 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",lu->id.RINFOG(3));CHKERRQ(ierr); 1128 1129 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(3));CHKERRQ(ierr); 1130 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",lu->id.INFOG(4));CHKERRQ(ierr); 1131 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",lu->id.INFOG(5));CHKERRQ(ierr); 1132 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",lu->id.INFOG(6));CHKERRQ(ierr); 1133 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively uese after analysis): %d \n",lu->id.INFOG(7));CHKERRQ(ierr); 1134 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",lu->id.INFOG(8));CHKERRQ(ierr); 1135 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",lu->id.INFOG(9));CHKERRQ(ierr); 1136 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",lu->id.INFOG(10));CHKERRQ(ierr); 1137 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",lu->id.INFOG(11));CHKERRQ(ierr); 1138 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",lu->id.INFOG(12));CHKERRQ(ierr); 1139 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",lu->id.INFOG(13));CHKERRQ(ierr); 1140 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",lu->id.INFOG(14));CHKERRQ(ierr); 1141 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",lu->id.INFOG(15));CHKERRQ(ierr); 1142 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",lu->id.INFOG(16));CHKERRQ(ierr); 1143 ierr = PetscViewerASCIIPrintf(viewer," INFOG(17) (estimated size of all MUMPS internal data for factorization after analysis: sum over all processors): %d \n",lu->id.INFOG(17));CHKERRQ(ierr); 1144 ierr = PetscViewerASCIIPrintf(viewer," INFOG(18) (size of all MUMPS internal data allocated during factorization: value on the most memory consuming processor): %d \n",lu->id.INFOG(18));CHKERRQ(ierr); 1145 ierr = PetscViewerASCIIPrintf(viewer," INFOG(19) (size of all MUMPS internal data allocated during factorization: sum over all processors): %d \n",lu->id.INFOG(19));CHKERRQ(ierr); 1146 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",lu->id.INFOG(20));CHKERRQ(ierr); 1147 ierr = PetscViewerASCIIPrintf(viewer," INFOG(21) (size in MB of memory effectively used during factorization - value on the most memory consuming processor): %d \n",lu->id.INFOG(21));CHKERRQ(ierr); 1148 ierr = PetscViewerASCIIPrintf(viewer," INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): %d \n",lu->id.INFOG(22));CHKERRQ(ierr); 1149 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",lu->id.INFOG(23));CHKERRQ(ierr); 1150 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",lu->id.INFOG(24));CHKERRQ(ierr); 1151 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",lu->id.INFOG(25));CHKERRQ(ierr); 1152 } 1153 PetscFunctionReturn(0); 1154 } 1155 1156 #undef __FUNCT__ 1157 #define __FUNCT__ "MatView_MUMPS" 1158 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1159 { 1160 PetscErrorCode ierr; 1161 PetscTruth iascii; 1162 PetscViewerFormat format; 1163 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1164 1165 PetscFunctionBegin; 1166 /* check if matrix is mumps type */ 1167 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1168 1169 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 1170 if (iascii) { 1171 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1172 if (format == PETSC_VIEWER_ASCII_INFO || mumps->mumpsview){ 1173 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1174 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1175 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1176 if (mumps->mumpsview){ /* View all MUMPS parameters */ 1177 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 1178 } 1179 } 1180 } 1181 PetscFunctionReturn(0); 1182 } 1183 1184 #undef __FUNCT__ 1185 #define __FUNCT__ "MatGetInfo_MUMPS" 1186 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1187 { 1188 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1189 1190 PetscFunctionBegin; 1191 info->block_size = 1.0; 1192 info->nz_allocated = mumps->id.INFOG(20); 1193 info->nz_used = mumps->id.INFOG(20); 1194 info->nz_unneeded = 0.0; 1195 info->assemblies = 0.0; 1196 info->mallocs = 0.0; 1197 info->memory = 0.0; 1198 info->fill_ratio_given = 0; 1199 info->fill_ratio_needed = 0; 1200 info->factor_mallocs = 0; 1201 PetscFunctionReturn(0); 1202 } 1203 1204 /*MC 1205 MAT_SOLVER_MUMPS - A matrix type providing direct solvers (LU and Cholesky) for 1206 distributed and sequential matrices via the external package MUMPS. 1207 1208 Works with MATAIJ and MATSBAIJ matrices 1209 1210 Options Database Keys: 1211 + -mat_mumps_sym <0,1,2> - 0 the matrix is unsymmetric, 1 symmetric positive definite, 2 symmetric 1212 . -mat_mumps_icntl_4 <0,...,4> - print level 1213 . -mat_mumps_icntl_6 <0,...,7> - matrix prescaling options (see MUMPS User's Guide) 1214 . -mat_mumps_icntl_7 <0,...,7> - matrix orderings (see MUMPS User's Guide) 1215 . -mat_mumps_icntl_9 <1,2> - A or A^T x=b to be solved: 1 denotes A, 2 denotes A^T 1216 . -mat_mumps_icntl_10 <n> - maximum number of iterative refinements 1217 . -mat_mumps_icntl_11 <n> - error analysis, a positive value returns statistics during -ksp_view 1218 . -mat_mumps_icntl_12 <n> - efficiency control (see MUMPS User's Guide) 1219 . -mat_mumps_icntl_13 <n> - efficiency control (see MUMPS User's Guide) 1220 . -mat_mumps_icntl_14 <n> - efficiency control (see MUMPS User's Guide) 1221 . -mat_mumps_icntl_15 <n> - efficiency control (see MUMPS User's Guide) 1222 . -mat_mumps_cntl_1 <delta> - relative pivoting threshold 1223 . -mat_mumps_cntl_2 <tol> - stopping criterion for refinement 1224 - -mat_mumps_cntl_3 <adelta> - absolute pivoting threshold 1225 1226 Level: beginner 1227 1228 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 1229 1230 M*/ 1231 1232 EXTERN_C_BEGIN 1233 #undef __FUNCT__ 1234 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 1235 PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 1236 { 1237 PetscFunctionBegin; 1238 *type = MAT_SOLVER_MUMPS; 1239 PetscFunctionReturn(0); 1240 } 1241 EXTERN_C_END 1242 1243 EXTERN_C_BEGIN 1244 /* 1245 The seq and mpi versions of this function are the same 1246 */ 1247 #undef __FUNCT__ 1248 #define __FUNCT__ "MatGetFactor_seqaij_mumps" 1249 PetscErrorCode MatGetFactor_seqaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1250 { 1251 Mat B; 1252 PetscErrorCode ierr; 1253 Mat_MUMPS *mumps; 1254 1255 PetscFunctionBegin; 1256 /* Create the factorization matrix */ 1257 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1258 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1259 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1260 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1261 1262 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1263 B->ops->view = MatView_MUMPS; 1264 B->ops->getinfo = MatGetInfo_MUMPS; 1265 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1266 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1267 if (ftype == MAT_FACTOR_LU) { 1268 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1269 B->factortype = MAT_FACTOR_LU; 1270 mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 1271 } else { 1272 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1273 B->factortype = MAT_FACTOR_CHOLESKY; 1274 mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 1275 } 1276 1277 mumps->CleanUpMUMPS = PETSC_FALSE; 1278 mumps->isAIJ = PETSC_TRUE; 1279 mumps->scat_rhs = PETSC_NULL; 1280 mumps->scat_sol = PETSC_NULL; 1281 mumps->nSolve = 0; 1282 mumps->MatDestroy = B->ops->destroy; 1283 B->ops->destroy = MatDestroy_MUMPS; 1284 B->spptr = (void*)mumps; 1285 1286 *F = B; 1287 PetscFunctionReturn(0); 1288 } 1289 EXTERN_C_END 1290 1291 EXTERN_C_BEGIN 1292 #undef __FUNCT__ 1293 #define __FUNCT__ "MatGetFactor_seqsbaij_mumps" 1294 PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1295 { 1296 Mat B; 1297 PetscErrorCode ierr; 1298 Mat_MUMPS *mumps; 1299 1300 PetscFunctionBegin; 1301 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1302 /* Create the factorization matrix */ 1303 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1304 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1305 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1306 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1307 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1308 1309 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1310 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1311 B->ops->view = MatView_MUMPS; 1312 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1313 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1314 B->factortype = MAT_FACTOR_CHOLESKY; 1315 1316 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 1317 mumps->CleanUpMUMPS = PETSC_FALSE; 1318 mumps->isAIJ = PETSC_FALSE; 1319 mumps->scat_rhs = PETSC_NULL; 1320 mumps->scat_sol = PETSC_NULL; 1321 mumps->nSolve = 0; 1322 mumps->MatDestroy = B->ops->destroy; 1323 B->ops->destroy = MatDestroy_MUMPS; 1324 B->spptr = (void*)mumps; 1325 1326 *F = B; 1327 PetscFunctionReturn(0); 1328 } 1329 EXTERN_C_END 1330 1331 EXTERN_C_BEGIN 1332 #undef __FUNCT__ 1333 #define __FUNCT__ "MatGetFactor_mpisbaij_mumps" 1334 PetscErrorCode MatGetFactor_mpisbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1335 { 1336 Mat B; 1337 PetscErrorCode ierr; 1338 Mat_MUMPS *mumps; 1339 1340 PetscFunctionBegin; 1341 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 1342 /* Create the factorization matrix */ 1343 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1344 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1345 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1346 ierr = MatSeqSBAIJSetPreallocation(B,1,0,PETSC_NULL);CHKERRQ(ierr); 1347 ierr = MatMPISBAIJSetPreallocation(B,1,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1348 1349 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1350 B->ops->view = MatView_MUMPS; 1351 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1352 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1353 B->factortype = MAT_FACTOR_CHOLESKY; 1354 1355 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1356 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 1357 mumps->CleanUpMUMPS = PETSC_FALSE; 1358 mumps->isAIJ = PETSC_FALSE; 1359 mumps->scat_rhs = PETSC_NULL; 1360 mumps->scat_sol = PETSC_NULL; 1361 mumps->nSolve = 0; 1362 mumps->MatDestroy = B->ops->destroy; 1363 B->ops->destroy = MatDestroy_MUMPS; 1364 B->spptr = (void*)mumps; 1365 1366 *F = B; 1367 PetscFunctionReturn(0); 1368 } 1369 EXTERN_C_END 1370 1371 EXTERN_C_BEGIN 1372 #undef __FUNCT__ 1373 #define __FUNCT__ "MatGetFactor_mpiaij_mumps" 1374 PetscErrorCode MatGetFactor_mpiaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1375 { 1376 Mat B; 1377 PetscErrorCode ierr; 1378 Mat_MUMPS *mumps; 1379 1380 PetscFunctionBegin; 1381 /* Create the factorization matrix */ 1382 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1383 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1384 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1385 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 1386 ierr = MatMPIAIJSetPreallocation(B,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1387 1388 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1389 1390 if (ftype == MAT_FACTOR_LU) { 1391 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 1392 B->factortype = MAT_FACTOR_LU; 1393 mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 1394 } else { 1395 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 1396 B->factortype = MAT_FACTOR_CHOLESKY; 1397 mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 1398 } 1399 1400 B->ops->view = MatView_MUMPS; 1401 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1402 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1403 1404 mumps->CleanUpMUMPS = PETSC_FALSE; 1405 mumps->isAIJ = PETSC_TRUE; 1406 mumps->scat_rhs = PETSC_NULL; 1407 mumps->scat_sol = PETSC_NULL; 1408 mumps->nSolve = 0; 1409 mumps->MatDestroy = B->ops->destroy; 1410 B->ops->destroy = MatDestroy_MUMPS; 1411 B->spptr = (void*)mumps; 1412 1413 *F = B; 1414 PetscFunctionReturn(0); 1415 } 1416 EXTERN_C_END 1417 1418 EXTERN_C_BEGIN 1419 #undef __FUNCT__ 1420 #define __FUNCT__ "MatGetFactor_seqbaij_mumps" 1421 PetscErrorCode MatGetFactor_seqbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1422 { 1423 Mat B; 1424 PetscErrorCode ierr; 1425 Mat_MUMPS *mumps; 1426 1427 PetscFunctionBegin; 1428 /* Create the factorization matrix */ 1429 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1430 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1431 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1432 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,PETSC_NULL);CHKERRQ(ierr); 1433 1434 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1435 if (ftype == MAT_FACTOR_LU) { 1436 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1437 B->factortype = MAT_FACTOR_LU; 1438 mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 1439 } 1440 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cholesky factorization for BAIJ matrices not supported, use AIJ matrix instead\n"); 1441 1442 B->ops->view = MatView_MUMPS; 1443 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1444 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1445 1446 mumps->CleanUpMUMPS = PETSC_FALSE; 1447 mumps->isAIJ = PETSC_TRUE; 1448 mumps->scat_rhs = PETSC_NULL; 1449 mumps->scat_sol = PETSC_NULL; 1450 mumps->nSolve = 0; 1451 mumps->MatDestroy = B->ops->destroy; 1452 B->ops->destroy = MatDestroy_MUMPS; 1453 B->spptr = (void*)mumps; 1454 1455 *F = B; 1456 PetscFunctionReturn(0); 1457 } 1458 EXTERN_C_END 1459 1460 EXTERN_C_BEGIN 1461 #undef __FUNCT__ 1462 #define __FUNCT__ "MatGetFactor_mpibaij_mumps" 1463 PetscErrorCode MatGetFactor_mpibaij_mumps(Mat A,MatFactorType ftype,Mat *F) 1464 { 1465 Mat B; 1466 PetscErrorCode ierr; 1467 Mat_MUMPS *mumps; 1468 1469 PetscFunctionBegin; 1470 /* Create the factorization matrix */ 1471 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 1472 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1473 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1474 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,PETSC_NULL,0,PETSC_NULL);CHKERRQ(ierr); 1475 1476 ierr = PetscNewLog(B,Mat_MUMPS,&mumps);CHKERRQ(ierr); 1477 if (ftype == MAT_FACTOR_LU) { 1478 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 1479 B->factortype = MAT_FACTOR_LU; 1480 mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 1481 } 1482 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cholesky factorization for BAIJ matrices not supported, use AIJ matrix instead\n"); 1483 B->ops->view = MatView_MUMPS; 1484 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mumps",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 1485 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMumpsSetIcntl_C","MatMumpsSetIcntl",MatMumpsSetIcntl);CHKERRQ(ierr); 1486 1487 mumps->CleanUpMUMPS = PETSC_FALSE; 1488 mumps->isAIJ = PETSC_TRUE; 1489 mumps->scat_rhs = PETSC_NULL; 1490 mumps->scat_sol = PETSC_NULL; 1491 mumps->nSolve = 0; 1492 mumps->MatDestroy = B->ops->destroy; 1493 B->ops->destroy = MatDestroy_MUMPS; 1494 B->spptr = (void*)mumps; 1495 1496 *F = B; 1497 PetscFunctionReturn(0); 1498 } 1499 EXTERN_C_END 1500 1501 /* -------------------------------------------------------------------------------------------*/ 1502 /*@ 1503 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 1504 1505 Collective on Mat 1506 1507 Input Parameters: 1508 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1509 . idx - index of MUMPS parameter array ICNTL() 1510 - icntl - value of MUMPS ICNTL(imumps) 1511 1512 Options Database: 1513 . -mat_mumps_icntl_<idx> <icntl> 1514 1515 Level: beginner 1516 1517 References: MUMPS Users' Guide 1518 1519 .seealso: MatGetFactor() 1520 @*/ 1521 #undef __FUNCT__ 1522 #define __FUNCT__ "MatMumpsSetIcntl" 1523 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt idx,PetscInt icntl) 1524 { 1525 Mat_MUMPS *lu =(Mat_MUMPS*)(F)->spptr; 1526 1527 PetscFunctionBegin; 1528 lu->id.ICNTL(idx) = icntl; 1529 PetscFunctionReturn(0); 1530 } 1531 1532