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