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