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