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