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