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