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