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