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