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