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