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 if (mumps->size > 1) { 1175 PetscInt lsol_loc; 1176 PetscScalar *sol_loc; 1177 1178 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isMPIAIJ);CHKERRQ(ierr); 1179 if (isMPIAIJ) F_diag = ((Mat_MPIAIJ*)(F)->data)->A; 1180 else F_diag = ((Mat_MPISBAIJ*)(F)->data)->A; 1181 F_diag->assembled = PETSC_TRUE; 1182 1183 /* distributed solution; Create x_seq=sol_loc for repeated use */ 1184 if (mumps->x_seq) { 1185 ierr = VecScatterDestroy(&mumps->scat_sol);CHKERRQ(ierr); 1186 ierr = PetscFree2(mumps->id.sol_loc,mumps->id.isol_loc);CHKERRQ(ierr); 1187 ierr = VecDestroy(&mumps->x_seq);CHKERRQ(ierr); 1188 } 1189 lsol_loc = mumps->id.INFO(23); /* length of sol_loc */ 1190 ierr = PetscMalloc2(lsol_loc,&sol_loc,lsol_loc,&mumps->id.isol_loc);CHKERRQ(ierr); 1191 mumps->id.lsol_loc = lsol_loc; 1192 mumps->id.sol_loc = (MumpsScalar*)sol_loc; 1193 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,lsol_loc,sol_loc,&mumps->x_seq);CHKERRQ(ierr); 1194 } 1195 PetscFunctionReturn(0); 1196 } 1197 1198 /* Sets MUMPS options from the options database */ 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "PetscSetMUMPSFromOptions" 1201 PetscErrorCode PetscSetMUMPSFromOptions(Mat F, Mat A) 1202 { 1203 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1204 PetscErrorCode ierr; 1205 PetscInt icntl,info[40],i,ninfo=40; 1206 PetscBool flg; 1207 1208 PetscFunctionBegin; 1209 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MUMPS Options","Mat");CHKERRQ(ierr); 1210 ierr = PetscOptionsInt("-mat_mumps_icntl_1","ICNTL(1): output stream for error messages","None",mumps->id.ICNTL(1),&icntl,&flg);CHKERRQ(ierr); 1211 if (flg) mumps->id.ICNTL(1) = icntl; 1212 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); 1213 if (flg) mumps->id.ICNTL(2) = icntl; 1214 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); 1215 if (flg) mumps->id.ICNTL(3) = icntl; 1216 1217 ierr = PetscOptionsInt("-mat_mumps_icntl_4","ICNTL(4): level of printing (0 to 4)","None",mumps->id.ICNTL(4),&icntl,&flg);CHKERRQ(ierr); 1218 if (flg) mumps->id.ICNTL(4) = icntl; 1219 if (mumps->id.ICNTL(4) || PetscLogPrintInfo) mumps->id.ICNTL(3) = 6; /* resume MUMPS default id.ICNTL(3) = 6 */ 1220 1221 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); 1222 if (flg) mumps->id.ICNTL(6) = icntl; 1223 1224 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); 1225 if (flg) { 1226 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"); 1227 else mumps->id.ICNTL(7) = icntl; 1228 } 1229 1230 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); 1231 /* 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() */ 1232 ierr = PetscOptionsInt("-mat_mumps_icntl_10","ICNTL(10): max num of refinements","None",mumps->id.ICNTL(10),&mumps->id.ICNTL(10),NULL);CHKERRQ(ierr); 1233 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); 1234 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); 1235 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); 1236 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); 1237 ierr = PetscOptionsInt("-mat_mumps_icntl_19","ICNTL(19): computes the Schur complement","None",mumps->id.ICNTL(19),&mumps->id.ICNTL(19),NULL);CHKERRQ(ierr); 1238 if (mumps->id.ICNTL(19) <= 0 || mumps->id.ICNTL(19) > 3) { /* reset any schur data (if any) */ 1239 ierr = MatMumpsResetSchur_Private(mumps);CHKERRQ(ierr); 1240 } 1241 /* 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 */ 1242 /* 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 */ 1243 1244 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); 1245 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); 1246 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); 1247 if (mumps->id.ICNTL(24)) { 1248 mumps->id.ICNTL(13) = 1; /* turn-off ScaLAPACK to help with the correct detection of null pivots */ 1249 } 1250 1251 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); 1252 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); 1253 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); 1254 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); 1255 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); 1256 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); 1257 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); 1258 /* 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 */ 1259 ierr = PetscOptionsInt("-mat_mumps_icntl_33","ICNTL(33): compute determinant","None",mumps->id.ICNTL(33),&mumps->id.ICNTL(33),NULL);CHKERRQ(ierr); 1260 1261 ierr = PetscOptionsReal("-mat_mumps_cntl_1","CNTL(1): relative pivoting threshold","None",mumps->id.CNTL(1),&mumps->id.CNTL(1),NULL);CHKERRQ(ierr); 1262 ierr = PetscOptionsReal("-mat_mumps_cntl_2","CNTL(2): stopping criterion of refinement","None",mumps->id.CNTL(2),&mumps->id.CNTL(2),NULL);CHKERRQ(ierr); 1263 ierr = PetscOptionsReal("-mat_mumps_cntl_3","CNTL(3): absolute pivoting threshold","None",mumps->id.CNTL(3),&mumps->id.CNTL(3),NULL);CHKERRQ(ierr); 1264 ierr = PetscOptionsReal("-mat_mumps_cntl_4","CNTL(4): value for static pivoting","None",mumps->id.CNTL(4),&mumps->id.CNTL(4),NULL);CHKERRQ(ierr); 1265 ierr = PetscOptionsReal("-mat_mumps_cntl_5","CNTL(5): fixation for null pivots","None",mumps->id.CNTL(5),&mumps->id.CNTL(5),NULL);CHKERRQ(ierr); 1266 1267 ierr = PetscOptionsString("-mat_mumps_ooc_tmpdir", "out of core directory", "None", mumps->id.ooc_tmpdir, mumps->id.ooc_tmpdir, 256, NULL); 1268 1269 ierr = PetscOptionsIntArray("-mat_mumps_view_info","request INFO local to each processor","",info,&ninfo,NULL);CHKERRQ(ierr); 1270 if (ninfo) { 1271 if (ninfo > 40) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"number of INFO %d must <= 40\n",ninfo); 1272 ierr = PetscMalloc1(ninfo,&mumps->info);CHKERRQ(ierr); 1273 mumps->ninfo = ninfo; 1274 for (i=0; i<ninfo; i++) { 1275 if (info[i] < 0 || info[i]>40) { 1276 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index of INFO %d must between 1 and 40\n",ninfo); 1277 } else { 1278 mumps->info[i] = info[i]; 1279 } 1280 } 1281 } 1282 1283 PetscOptionsEnd(); 1284 PetscFunctionReturn(0); 1285 } 1286 1287 #undef __FUNCT__ 1288 #define __FUNCT__ "PetscInitializeMUMPS" 1289 PetscErrorCode PetscInitializeMUMPS(Mat A,Mat_MUMPS *mumps) 1290 { 1291 PetscErrorCode ierr; 1292 1293 PetscFunctionBegin; 1294 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)A), &mumps->myid); 1295 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&mumps->size);CHKERRQ(ierr); 1296 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mumps->comm_mumps));CHKERRQ(ierr); 1297 1298 mumps->id.comm_fortran = MPI_Comm_c2f(mumps->comm_mumps); 1299 1300 mumps->id.job = JOB_INIT; 1301 mumps->id.par = 1; /* host participates factorizaton and solve */ 1302 mumps->id.sym = mumps->sym; 1303 PetscMUMPS_c(&mumps->id); 1304 1305 mumps->CleanUpMUMPS = PETSC_FALSE; 1306 mumps->scat_rhs = NULL; 1307 mumps->scat_sol = NULL; 1308 1309 /* set PETSc-MUMPS default options - override MUMPS default */ 1310 mumps->id.ICNTL(3) = 0; 1311 mumps->id.ICNTL(4) = 0; 1312 if (mumps->size == 1) { 1313 mumps->id.ICNTL(18) = 0; /* centralized assembled matrix input */ 1314 } else { 1315 mumps->id.ICNTL(18) = 3; /* distributed assembled matrix input */ 1316 mumps->id.ICNTL(20) = 0; /* rhs is in dense format */ 1317 mumps->id.ICNTL(21) = 1; /* distributed solution */ 1318 } 1319 1320 /* schur */ 1321 mumps->id.size_schur = 0; 1322 mumps->id.listvar_schur = NULL; 1323 mumps->id.schur = NULL; 1324 mumps->schur_second_solve = PETSC_FALSE; 1325 mumps->sizeredrhs = 0; 1326 mumps->schur_pivots = NULL; 1327 mumps->schur_work = NULL; 1328 mumps->schur_sol = NULL; 1329 mumps->schur_sizesol = 0; 1330 mumps->schur_restored = PETSC_TRUE; 1331 mumps->schur_factored = PETSC_FALSE; 1332 mumps->schur_inverted = PETSC_FALSE; 1333 PetscFunctionReturn(0); 1334 } 1335 1336 /* Note Petsc r(=c) permutation is used when mumps->id.ICNTL(7)==1 with centralized assembled matrix input; otherwise r and c are ignored */ 1337 #undef __FUNCT__ 1338 #define __FUNCT__ "MatLUFactorSymbolic_AIJMUMPS" 1339 PetscErrorCode MatLUFactorSymbolic_AIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1340 { 1341 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1342 PetscErrorCode ierr; 1343 Vec b; 1344 IS is_iden; 1345 const PetscInt M = A->rmap->N; 1346 1347 PetscFunctionBegin; 1348 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1349 1350 /* Set MUMPS options from the options database */ 1351 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1352 1353 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1354 1355 /* analysis phase */ 1356 /*----------------*/ 1357 mumps->id.job = JOB_FACTSYMBOLIC; 1358 mumps->id.n = M; 1359 switch (mumps->id.ICNTL(18)) { 1360 case 0: /* centralized assembled matrix input */ 1361 if (!mumps->myid) { 1362 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1363 if (mumps->id.ICNTL(6)>1) { 1364 mumps->id.a = (MumpsScalar*)mumps->val; 1365 } 1366 if (mumps->id.ICNTL(7) == 1) { /* use user-provide matrix ordering - assuming r = c ordering */ 1367 /* 1368 PetscBool flag; 1369 ierr = ISEqual(r,c,&flag);CHKERRQ(ierr); 1370 if (!flag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"row_perm != col_perm"); 1371 ierr = ISView(r,PETSC_VIEWER_STDOUT_SELF); 1372 */ 1373 if (!mumps->myid) { 1374 const PetscInt *idx; 1375 PetscInt i,*perm_in; 1376 1377 ierr = PetscMalloc1(M,&perm_in);CHKERRQ(ierr); 1378 ierr = ISGetIndices(r,&idx);CHKERRQ(ierr); 1379 1380 mumps->id.perm_in = perm_in; 1381 for (i=0; i<M; i++) perm_in[i] = idx[i]+1; /* perm_in[]: start from 1, not 0! */ 1382 ierr = ISRestoreIndices(r,&idx);CHKERRQ(ierr); 1383 } 1384 } 1385 } 1386 break; 1387 case 3: /* distributed assembled matrix input (size>1) */ 1388 mumps->id.nz_loc = mumps->nz; 1389 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1390 if (mumps->id.ICNTL(6)>1) { 1391 mumps->id.a_loc = (MumpsScalar*)mumps->val; 1392 } 1393 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1394 if (!mumps->myid) { 1395 ierr = VecCreateSeq(PETSC_COMM_SELF,A->rmap->N,&mumps->b_seq);CHKERRQ(ierr); 1396 ierr = ISCreateStride(PETSC_COMM_SELF,A->rmap->N,0,1,&is_iden);CHKERRQ(ierr); 1397 } else { 1398 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1399 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1400 } 1401 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1402 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1403 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1404 ierr = VecDestroy(&b);CHKERRQ(ierr); 1405 break; 1406 } 1407 PetscMUMPS_c(&mumps->id); 1408 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)); 1409 1410 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1411 F->ops->solve = MatSolve_MUMPS; 1412 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1413 F->ops->matsolve = MatMatSolve_MUMPS; 1414 PetscFunctionReturn(0); 1415 } 1416 1417 /* Note the Petsc r and c permutations are ignored */ 1418 #undef __FUNCT__ 1419 #define __FUNCT__ "MatLUFactorSymbolic_BAIJMUMPS" 1420 PetscErrorCode MatLUFactorSymbolic_BAIJMUMPS(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1421 { 1422 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1423 PetscErrorCode ierr; 1424 Vec b; 1425 IS is_iden; 1426 const PetscInt M = A->rmap->N; 1427 1428 PetscFunctionBegin; 1429 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1430 1431 /* Set MUMPS options from the options database */ 1432 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1433 1434 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1435 1436 /* analysis phase */ 1437 /*----------------*/ 1438 mumps->id.job = JOB_FACTSYMBOLIC; 1439 mumps->id.n = M; 1440 switch (mumps->id.ICNTL(18)) { 1441 case 0: /* centralized assembled matrix input */ 1442 if (!mumps->myid) { 1443 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1444 if (mumps->id.ICNTL(6)>1) { 1445 mumps->id.a = (MumpsScalar*)mumps->val; 1446 } 1447 } 1448 break; 1449 case 3: /* distributed assembled matrix input (size>1) */ 1450 mumps->id.nz_loc = mumps->nz; 1451 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1452 if (mumps->id.ICNTL(6)>1) { 1453 mumps->id.a_loc = (MumpsScalar*)mumps->val; 1454 } 1455 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1456 if (!mumps->myid) { 1457 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1458 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1459 } else { 1460 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1461 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1462 } 1463 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1464 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1465 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1466 ierr = VecDestroy(&b);CHKERRQ(ierr); 1467 break; 1468 } 1469 PetscMUMPS_c(&mumps->id); 1470 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)); 1471 1472 F->ops->lufactornumeric = MatFactorNumeric_MUMPS; 1473 F->ops->solve = MatSolve_MUMPS; 1474 F->ops->solvetranspose = MatSolveTranspose_MUMPS; 1475 PetscFunctionReturn(0); 1476 } 1477 1478 /* Note the Petsc r permutation and factor info are ignored */ 1479 #undef __FUNCT__ 1480 #define __FUNCT__ "MatCholeskyFactorSymbolic_MUMPS" 1481 PetscErrorCode MatCholeskyFactorSymbolic_MUMPS(Mat F,Mat A,IS r,const MatFactorInfo *info) 1482 { 1483 Mat_MUMPS *mumps = (Mat_MUMPS*)F->spptr; 1484 PetscErrorCode ierr; 1485 Vec b; 1486 IS is_iden; 1487 const PetscInt M = A->rmap->N; 1488 1489 PetscFunctionBegin; 1490 mumps->matstruc = DIFFERENT_NONZERO_PATTERN; 1491 1492 /* Set MUMPS options from the options database */ 1493 ierr = PetscSetMUMPSFromOptions(F,A);CHKERRQ(ierr); 1494 1495 ierr = (*mumps->ConvertToTriples)(A, 1, MAT_INITIAL_MATRIX, &mumps->nz, &mumps->irn, &mumps->jcn, &mumps->val);CHKERRQ(ierr); 1496 1497 /* analysis phase */ 1498 /*----------------*/ 1499 mumps->id.job = JOB_FACTSYMBOLIC; 1500 mumps->id.n = M; 1501 switch (mumps->id.ICNTL(18)) { 1502 case 0: /* centralized assembled matrix input */ 1503 if (!mumps->myid) { 1504 mumps->id.nz =mumps->nz; mumps->id.irn=mumps->irn; mumps->id.jcn=mumps->jcn; 1505 if (mumps->id.ICNTL(6)>1) { 1506 mumps->id.a = (MumpsScalar*)mumps->val; 1507 } 1508 } 1509 break; 1510 case 3: /* distributed assembled matrix input (size>1) */ 1511 mumps->id.nz_loc = mumps->nz; 1512 mumps->id.irn_loc=mumps->irn; mumps->id.jcn_loc=mumps->jcn; 1513 if (mumps->id.ICNTL(6)>1) { 1514 mumps->id.a_loc = (MumpsScalar*)mumps->val; 1515 } 1516 /* MUMPS only supports centralized rhs. Create scatter scat_rhs for repeated use in MatSolve() */ 1517 if (!mumps->myid) { 1518 ierr = VecCreateSeq(PETSC_COMM_SELF,A->cmap->N,&mumps->b_seq);CHKERRQ(ierr); 1519 ierr = ISCreateStride(PETSC_COMM_SELF,A->cmap->N,0,1,&is_iden);CHKERRQ(ierr); 1520 } else { 1521 ierr = VecCreateSeq(PETSC_COMM_SELF,0,&mumps->b_seq);CHKERRQ(ierr); 1522 ierr = ISCreateStride(PETSC_COMM_SELF,0,0,1,&is_iden);CHKERRQ(ierr); 1523 } 1524 ierr = MatCreateVecs(A,NULL,&b);CHKERRQ(ierr); 1525 ierr = VecScatterCreate(b,is_iden,mumps->b_seq,is_iden,&mumps->scat_rhs);CHKERRQ(ierr); 1526 ierr = ISDestroy(&is_iden);CHKERRQ(ierr); 1527 ierr = VecDestroy(&b);CHKERRQ(ierr); 1528 break; 1529 } 1530 PetscMUMPS_c(&mumps->id); 1531 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)); 1532 1533 F->ops->choleskyfactornumeric = MatFactorNumeric_MUMPS; 1534 F->ops->solve = MatSolve_MUMPS; 1535 F->ops->solvetranspose = MatSolve_MUMPS; 1536 F->ops->matsolve = MatMatSolve_MUMPS; 1537 #if defined(PETSC_USE_COMPLEX) 1538 F->ops->getinertia = NULL; 1539 #else 1540 F->ops->getinertia = MatGetInertia_SBAIJMUMPS; 1541 #endif 1542 PetscFunctionReturn(0); 1543 } 1544 1545 #undef __FUNCT__ 1546 #define __FUNCT__ "MatView_MUMPS" 1547 PetscErrorCode MatView_MUMPS(Mat A,PetscViewer viewer) 1548 { 1549 PetscErrorCode ierr; 1550 PetscBool iascii; 1551 PetscViewerFormat format; 1552 Mat_MUMPS *mumps=(Mat_MUMPS*)A->spptr; 1553 1554 PetscFunctionBegin; 1555 /* check if matrix is mumps type */ 1556 if (A->ops->solve != MatSolve_MUMPS) PetscFunctionReturn(0); 1557 1558 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1559 if (iascii) { 1560 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1561 if (format == PETSC_VIEWER_ASCII_INFO) { 1562 ierr = PetscViewerASCIIPrintf(viewer,"MUMPS run parameters:\n");CHKERRQ(ierr); 1563 ierr = PetscViewerASCIIPrintf(viewer," SYM (matrix type): %d \n",mumps->id.sym);CHKERRQ(ierr); 1564 ierr = PetscViewerASCIIPrintf(viewer," PAR (host participation): %d \n",mumps->id.par);CHKERRQ(ierr); 1565 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(1) (output for error): %d \n",mumps->id.ICNTL(1));CHKERRQ(ierr); 1566 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(2) (output of diagnostic msg): %d \n",mumps->id.ICNTL(2));CHKERRQ(ierr); 1567 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(3) (output for global info): %d \n",mumps->id.ICNTL(3));CHKERRQ(ierr); 1568 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(4) (level of printing): %d \n",mumps->id.ICNTL(4));CHKERRQ(ierr); 1569 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(5) (input mat struct): %d \n",mumps->id.ICNTL(5));CHKERRQ(ierr); 1570 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(6) (matrix prescaling): %d \n",mumps->id.ICNTL(6));CHKERRQ(ierr); 1571 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(7) (sequentia matrix ordering):%d \n",mumps->id.ICNTL(7));CHKERRQ(ierr); 1572 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(8) (scalling strategy): %d \n",mumps->id.ICNTL(8));CHKERRQ(ierr); 1573 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(10) (max num of refinements): %d \n",mumps->id.ICNTL(10));CHKERRQ(ierr); 1574 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(11) (error analysis): %d \n",mumps->id.ICNTL(11));CHKERRQ(ierr); 1575 if (mumps->id.ICNTL(11)>0) { 1576 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(4) (inf norm of input mat): %g\n",mumps->id.RINFOG(4));CHKERRQ(ierr); 1577 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(5) (inf norm of solution): %g\n",mumps->id.RINFOG(5));CHKERRQ(ierr); 1578 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(6) (inf norm of residual): %g\n",mumps->id.RINFOG(6));CHKERRQ(ierr); 1579 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(7),RINFOG(8) (backward error est): %g, %g\n",mumps->id.RINFOG(7),mumps->id.RINFOG(8));CHKERRQ(ierr); 1580 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(9) (error estimate): %g \n",mumps->id.RINFOG(9));CHKERRQ(ierr); 1581 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(10),RINFOG(11)(condition numbers): %g, %g\n",mumps->id.RINFOG(10),mumps->id.RINFOG(11));CHKERRQ(ierr); 1582 } 1583 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(12) (efficiency control): %d \n",mumps->id.ICNTL(12));CHKERRQ(ierr); 1584 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(13) (efficiency control): %d \n",mumps->id.ICNTL(13));CHKERRQ(ierr); 1585 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(14) (percentage of estimated workspace increase): %d \n",mumps->id.ICNTL(14));CHKERRQ(ierr); 1586 /* ICNTL(15-17) not used */ 1587 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(18) (input mat struct): %d \n",mumps->id.ICNTL(18));CHKERRQ(ierr); 1588 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(19) (Shur complement info): %d \n",mumps->id.ICNTL(19));CHKERRQ(ierr); 1589 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(20) (rhs sparse pattern): %d \n",mumps->id.ICNTL(20));CHKERRQ(ierr); 1590 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(21) (solution struct): %d \n",mumps->id.ICNTL(21));CHKERRQ(ierr); 1591 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(22) (in-core/out-of-core facility): %d \n",mumps->id.ICNTL(22));CHKERRQ(ierr); 1592 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(23) (max size of memory can be allocated locally):%d \n",mumps->id.ICNTL(23));CHKERRQ(ierr); 1593 1594 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(24) (detection of null pivot rows): %d \n",mumps->id.ICNTL(24));CHKERRQ(ierr); 1595 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(25) (computation of a null space basis): %d \n",mumps->id.ICNTL(25));CHKERRQ(ierr); 1596 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(26) (Schur options for rhs or solution): %d \n",mumps->id.ICNTL(26));CHKERRQ(ierr); 1597 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(27) (experimental parameter): %d \n",mumps->id.ICNTL(27));CHKERRQ(ierr); 1598 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(28) (use parallel or sequential ordering): %d \n",mumps->id.ICNTL(28));CHKERRQ(ierr); 1599 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(29) (parallel ordering): %d \n",mumps->id.ICNTL(29));CHKERRQ(ierr); 1600 1601 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(30) (user-specified set of entries in inv(A)): %d \n",mumps->id.ICNTL(30));CHKERRQ(ierr); 1602 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(31) (factors is discarded in the solve phase): %d \n",mumps->id.ICNTL(31));CHKERRQ(ierr); 1603 ierr = PetscViewerASCIIPrintf(viewer," ICNTL(33) (compute determinant): %d \n",mumps->id.ICNTL(33));CHKERRQ(ierr); 1604 1605 ierr = PetscViewerASCIIPrintf(viewer," CNTL(1) (relative pivoting threshold): %g \n",mumps->id.CNTL(1));CHKERRQ(ierr); 1606 ierr = PetscViewerASCIIPrintf(viewer," CNTL(2) (stopping criterion of refinement): %g \n",mumps->id.CNTL(2));CHKERRQ(ierr); 1607 ierr = PetscViewerASCIIPrintf(viewer," CNTL(3) (absolute pivoting threshold): %g \n",mumps->id.CNTL(3));CHKERRQ(ierr); 1608 ierr = PetscViewerASCIIPrintf(viewer," CNTL(4) (value of static pivoting): %g \n",mumps->id.CNTL(4));CHKERRQ(ierr); 1609 ierr = PetscViewerASCIIPrintf(viewer," CNTL(5) (fixation for null pivots): %g \n",mumps->id.CNTL(5));CHKERRQ(ierr); 1610 1611 /* infomation local to each processor */ 1612 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(1) (local estimated flops for the elimination after analysis): \n");CHKERRQ(ierr); 1613 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 1614 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(1));CHKERRQ(ierr); 1615 ierr = PetscViewerFlush(viewer); 1616 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(2) (local estimated flops for the assembly after factorization): \n");CHKERRQ(ierr); 1617 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(2));CHKERRQ(ierr); 1618 ierr = PetscViewerFlush(viewer); 1619 ierr = PetscViewerASCIIPrintf(viewer, " RINFO(3) (local estimated flops for the elimination after factorization): \n");CHKERRQ(ierr); 1620 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %g \n",mumps->myid,mumps->id.RINFO(3));CHKERRQ(ierr); 1621 ierr = PetscViewerFlush(viewer); 1622 1623 ierr = PetscViewerASCIIPrintf(viewer, " INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization): \n");CHKERRQ(ierr); 1624 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(15));CHKERRQ(ierr); 1625 ierr = PetscViewerFlush(viewer); 1626 1627 ierr = PetscViewerASCIIPrintf(viewer, " INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization): \n");CHKERRQ(ierr); 1628 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(16));CHKERRQ(ierr); 1629 ierr = PetscViewerFlush(viewer); 1630 1631 ierr = PetscViewerASCIIPrintf(viewer, " INFO(23) (num of pivots eliminated on this processor after factorization): \n");CHKERRQ(ierr); 1632 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(23));CHKERRQ(ierr); 1633 ierr = PetscViewerFlush(viewer); 1634 1635 if (mumps->ninfo && mumps->ninfo <= 40){ 1636 PetscInt i; 1637 for (i=0; i<mumps->ninfo; i++){ 1638 ierr = PetscViewerASCIIPrintf(viewer, " INFO(%d): \n",mumps->info[i]);CHKERRQ(ierr); 1639 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] %d \n",mumps->myid,mumps->id.INFO(mumps->info[i]));CHKERRQ(ierr); 1640 ierr = PetscViewerFlush(viewer); 1641 } 1642 } 1643 1644 1645 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 1646 1647 if (!mumps->myid) { /* information from the host */ 1648 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(1) (global estimated flops for the elimination after analysis): %g \n",mumps->id.RINFOG(1));CHKERRQ(ierr); 1649 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(2) (global estimated flops for the assembly after factorization): %g \n",mumps->id.RINFOG(2));CHKERRQ(ierr); 1650 ierr = PetscViewerASCIIPrintf(viewer," RINFOG(3) (global estimated flops for the elimination after factorization): %g \n",mumps->id.RINFOG(3));CHKERRQ(ierr); 1651 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); 1652 1653 ierr = PetscViewerASCIIPrintf(viewer," INFOG(3) (estimated real workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(3));CHKERRQ(ierr); 1654 ierr = PetscViewerASCIIPrintf(viewer," INFOG(4) (estimated integer workspace for factors on all processors after analysis): %d \n",mumps->id.INFOG(4));CHKERRQ(ierr); 1655 ierr = PetscViewerASCIIPrintf(viewer," INFOG(5) (estimated maximum front size in the complete tree): %d \n",mumps->id.INFOG(5));CHKERRQ(ierr); 1656 ierr = PetscViewerASCIIPrintf(viewer," INFOG(6) (number of nodes in the complete tree): %d \n",mumps->id.INFOG(6));CHKERRQ(ierr); 1657 ierr = PetscViewerASCIIPrintf(viewer," INFOG(7) (ordering option effectively use after analysis): %d \n",mumps->id.INFOG(7));CHKERRQ(ierr); 1658 ierr = PetscViewerASCIIPrintf(viewer," INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): %d \n",mumps->id.INFOG(8));CHKERRQ(ierr); 1659 ierr = PetscViewerASCIIPrintf(viewer," INFOG(9) (total real/complex workspace to store the matrix factors after factorization): %d \n",mumps->id.INFOG(9));CHKERRQ(ierr); 1660 ierr = PetscViewerASCIIPrintf(viewer," INFOG(10) (total integer space store the matrix factors after factorization): %d \n",mumps->id.INFOG(10));CHKERRQ(ierr); 1661 ierr = PetscViewerASCIIPrintf(viewer," INFOG(11) (order of largest frontal matrix after factorization): %d \n",mumps->id.INFOG(11));CHKERRQ(ierr); 1662 ierr = PetscViewerASCIIPrintf(viewer," INFOG(12) (number of off-diagonal pivots): %d \n",mumps->id.INFOG(12));CHKERRQ(ierr); 1663 ierr = PetscViewerASCIIPrintf(viewer," INFOG(13) (number of delayed pivots after factorization): %d \n",mumps->id.INFOG(13));CHKERRQ(ierr); 1664 ierr = PetscViewerASCIIPrintf(viewer," INFOG(14) (number of memory compress after factorization): %d \n",mumps->id.INFOG(14));CHKERRQ(ierr); 1665 ierr = PetscViewerASCIIPrintf(viewer," INFOG(15) (number of steps of iterative refinement after solution): %d \n",mumps->id.INFOG(15));CHKERRQ(ierr); 1666 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); 1667 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); 1668 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); 1669 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); 1670 ierr = PetscViewerASCIIPrintf(viewer," INFOG(20) (estimated number of entries in the factors): %d \n",mumps->id.INFOG(20));CHKERRQ(ierr); 1671 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); 1672 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); 1673 ierr = PetscViewerASCIIPrintf(viewer," INFOG(23) (after analysis: value of ICNTL(6) effectively used): %d \n",mumps->id.INFOG(23));CHKERRQ(ierr); 1674 ierr = PetscViewerASCIIPrintf(viewer," INFOG(24) (after analysis: value of ICNTL(12) effectively used): %d \n",mumps->id.INFOG(24));CHKERRQ(ierr); 1675 ierr = PetscViewerASCIIPrintf(viewer," INFOG(25) (after factorization: number of pivots modified by static pivoting): %d \n",mumps->id.INFOG(25));CHKERRQ(ierr); 1676 ierr = PetscViewerASCIIPrintf(viewer," INFOG(28) (after factorization: number of null pivots encountered): %d\n",mumps->id.INFOG(28));CHKERRQ(ierr); 1677 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); 1678 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); 1679 ierr = PetscViewerASCIIPrintf(viewer," INFOG(32) (after analysis: type of analysis done): %d\n",mumps->id.INFOG(32));CHKERRQ(ierr); 1680 ierr = PetscViewerASCIIPrintf(viewer," INFOG(33) (value used for ICNTL(8)): %d\n",mumps->id.INFOG(33));CHKERRQ(ierr); 1681 ierr = PetscViewerASCIIPrintf(viewer," INFOG(34) (exponent of the determinant if determinant is requested): %d\n",mumps->id.INFOG(34));CHKERRQ(ierr); 1682 } 1683 } 1684 } 1685 PetscFunctionReturn(0); 1686 } 1687 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "MatGetInfo_MUMPS" 1690 PetscErrorCode MatGetInfo_MUMPS(Mat A,MatInfoType flag,MatInfo *info) 1691 { 1692 Mat_MUMPS *mumps =(Mat_MUMPS*)A->spptr; 1693 1694 PetscFunctionBegin; 1695 info->block_size = 1.0; 1696 info->nz_allocated = mumps->id.INFOG(20); 1697 info->nz_used = mumps->id.INFOG(20); 1698 info->nz_unneeded = 0.0; 1699 info->assemblies = 0.0; 1700 info->mallocs = 0.0; 1701 info->memory = 0.0; 1702 info->fill_ratio_given = 0; 1703 info->fill_ratio_needed = 0; 1704 info->factor_mallocs = 0; 1705 PetscFunctionReturn(0); 1706 } 1707 1708 /* -------------------------------------------------------------------------------------------*/ 1709 #undef __FUNCT__ 1710 #define __FUNCT__ "MatMumpsSetSchurIndices_MUMPS" 1711 PetscErrorCode MatMumpsSetSchurIndices_MUMPS(Mat F,PetscInt size,PetscInt idxs[]) 1712 { 1713 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1714 PetscErrorCode ierr; 1715 1716 PetscFunctionBegin; 1717 if (mumps->size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MUMPS parallel Schur complements not yet supported from PETSc\n"); 1718 if (mumps->id.size_schur != size) { 1719 ierr = PetscFree2(mumps->id.listvar_schur,mumps->id.schur);CHKERRQ(ierr); 1720 mumps->id.size_schur = size; 1721 mumps->id.schur_lld = size; 1722 ierr = PetscMalloc2(size,&mumps->id.listvar_schur,size*size,&mumps->id.schur);CHKERRQ(ierr); 1723 } 1724 ierr = PetscMemcpy(mumps->id.listvar_schur,idxs,size*sizeof(PetscInt));CHKERRQ(ierr); 1725 if (F->factortype == MAT_FACTOR_LU) { 1726 mumps->id.ICNTL(19) = 3; /* MUMPS returns full matrix */ 1727 } else { 1728 mumps->id.ICNTL(19) = 2; /* MUMPS returns lower triangular part */ 1729 } 1730 /* set a special value of ICNTL (not handled my MUMPS) to be used in the solve phase by PETSc */ 1731 mumps->id.ICNTL(26) = -1; 1732 PetscFunctionReturn(0); 1733 } 1734 1735 #undef __FUNCT__ 1736 #define __FUNCT__ "MatMumpsSetSchurIndices" 1737 /*@ 1738 MatMumpsSetSchurIndices - Set indices defining the Schur complement that MUMPS will compute during the factorization steps 1739 1740 Logically Collective on Mat 1741 1742 Input Parameters: 1743 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1744 . size - size of the Schur complement indices 1745 - idxs[] - array of Schur complement indices 1746 1747 Notes: 1748 The user has to free the array idxs[] since the indices are copied by the routine. 1749 MUMPS Schur complement mode is currently implemented for sequential matrices. 1750 1751 Level: advanced 1752 1753 References: MUMPS Users' Guide 1754 1755 .seealso: MatGetFactor(), MatMumpsCreateSchurComplement(), MatMumpsGetSchurComplement() 1756 @*/ 1757 PetscErrorCode MatMumpsSetSchurIndices(Mat F,PetscInt size,PetscInt idxs[]) 1758 { 1759 PetscErrorCode ierr; 1760 1761 PetscFunctionBegin; 1762 PetscValidHeaderSpecific(F,MAT_CLASSID,1); 1763 if (size) PetscValidIntPointer(idxs,3); 1764 ierr = PetscTryMethod(F,"MatMumpsSetSchurIndices_C",(Mat,PetscInt,PetscInt[]),(F,size,idxs));CHKERRQ(ierr); 1765 PetscFunctionReturn(0); 1766 } 1767 1768 /* -------------------------------------------------------------------------------------------*/ 1769 #undef __FUNCT__ 1770 #define __FUNCT__ "MatMumpsCreateSchurComplement_MUMPS" 1771 PetscErrorCode MatMumpsCreateSchurComplement_MUMPS(Mat F,Mat* S) 1772 { 1773 Mat St; 1774 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1775 PetscScalar *array; 1776 #if defined(PETSC_USE_COMPLEX) 1777 PetscScalar im = PetscSqrtScalar(-1.0); 1778 #endif 1779 PetscErrorCode ierr; 1780 1781 PetscFunctionBegin; 1782 if (!mumps->CleanUpMUMPS) { /* CleanUpMUMPS is set to true after numerical factorization */ 1783 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Numerical factorization phase not yet performed! You should call MatFactorSymbolic/Numeric before"); 1784 } else if (!mumps->id.ICNTL(19)) { 1785 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur complement mode not selected! You should call MatMumpsSetSchurIndices to enable it"); 1786 } else if (!mumps->id.size_schur) { 1787 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur indices not set! You should call MatMumpsSetSchurIndices before"); 1788 } else if (!mumps->schur_restored) { 1789 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur matrix has not been restored using MatMumpsRestoreSchurComplement"); 1790 } 1791 ierr = MatCreate(PetscObjectComm((PetscObject)F),&St);CHKERRQ(ierr); 1792 ierr = MatSetSizes(St,PETSC_DECIDE,PETSC_DECIDE,mumps->id.size_schur,mumps->id.size_schur);CHKERRQ(ierr); 1793 ierr = MatSetType(St,MATDENSE);CHKERRQ(ierr); 1794 ierr = MatSetUp(St);CHKERRQ(ierr); 1795 ierr = MatDenseGetArray(St,&array);CHKERRQ(ierr); 1796 if (!mumps->sym) { /* MUMPS always return a full matrix */ 1797 if (mumps->id.ICNTL(19) == 1) { /* stored by rows */ 1798 PetscInt i,j,N=mumps->id.size_schur; 1799 for (i=0;i<N;i++) { 1800 for (j=0;j<N;j++) { 1801 #if !defined(PETSC_USE_COMPLEX) 1802 PetscScalar val = mumps->id.schur[i*N+j]; 1803 #else 1804 PetscScalar val = mumps->id.schur[i*N+j].r + im*mumps->id.schur[i*N+j].i; 1805 #endif 1806 array[j*N+i] = val; 1807 } 1808 } 1809 } else { /* stored by columns */ 1810 ierr = PetscMemcpy(array,mumps->id.schur,mumps->id.size_schur*mumps->id.size_schur*sizeof(PetscScalar));CHKERRQ(ierr); 1811 } 1812 } else { /* either full or lower-triangular (not packed) */ 1813 if (mumps->id.ICNTL(19) == 2) { /* lower triangular stored by columns */ 1814 PetscInt i,j,N=mumps->id.size_schur; 1815 for (i=0;i<N;i++) { 1816 for (j=i;j<N;j++) { 1817 #if !defined(PETSC_USE_COMPLEX) 1818 PetscScalar val = mumps->id.schur[i*N+j]; 1819 #else 1820 PetscScalar val = mumps->id.schur[i*N+j].r + im*mumps->id.schur[i*N+j].i; 1821 #endif 1822 array[i*N+j] = val; 1823 } 1824 for (j=i;j<N;j++) { 1825 #if !defined(PETSC_USE_COMPLEX) 1826 PetscScalar val = mumps->id.schur[i*N+j]; 1827 #else 1828 PetscScalar val = mumps->id.schur[i*N+j].r + im*mumps->id.schur[i*N+j].i; 1829 #endif 1830 array[j*N+i] = val; 1831 } 1832 } 1833 } else if (mumps->id.ICNTL(19) == 3) { /* full matrix */ 1834 ierr = PetscMemcpy(array,mumps->id.schur,mumps->id.size_schur*mumps->id.size_schur*sizeof(PetscScalar));CHKERRQ(ierr); 1835 } else { /* ICNTL(19) == 1 lower triangular stored by rows */ 1836 PetscInt i,j,N=mumps->id.size_schur; 1837 for (i=0;i<N;i++) { 1838 for (j=0;j<i+1;j++) { 1839 #if !defined(PETSC_USE_COMPLEX) 1840 PetscScalar val = mumps->id.schur[i*N+j]; 1841 #else 1842 PetscScalar val = mumps->id.schur[i*N+j].r + im*mumps->id.schur[i*N+j].i; 1843 #endif 1844 array[i*N+j] = val; 1845 } 1846 for (j=0;j<i+1;j++) { 1847 #if !defined(PETSC_USE_COMPLEX) 1848 PetscScalar val = mumps->id.schur[i*N+j]; 1849 #else 1850 PetscScalar val = mumps->id.schur[i*N+j].r + im*mumps->id.schur[i*N+j].i; 1851 #endif 1852 array[j*N+i] = val; 1853 } 1854 } 1855 } 1856 } 1857 ierr = MatDenseRestoreArray(St,&array);CHKERRQ(ierr); 1858 *S = St; 1859 PetscFunctionReturn(0); 1860 } 1861 1862 #undef __FUNCT__ 1863 #define __FUNCT__ "MatMumpsCreateSchurComplement" 1864 /*@ 1865 MatMumpsCreateSchurComplement - Create a Schur complement matrix object using Schur data computed by MUMPS during the factorization step 1866 1867 Logically Collective on Mat 1868 1869 Input Parameters: 1870 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1871 . *S - location where to return the Schur complement (MATDENSE) 1872 1873 Notes: 1874 MUMPS Schur complement mode is currently implemented for sequential matrices. 1875 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. 1876 1877 Level: advanced 1878 1879 References: MUMPS Users' Guide 1880 1881 .seealso: MatGetFactor(), MatMumpsSetSchurIndices(), MatMumpsGetSchurComplement() 1882 @*/ 1883 PetscErrorCode MatMumpsCreateSchurComplement(Mat F,Mat* S) 1884 { 1885 PetscErrorCode ierr; 1886 1887 PetscFunctionBegin; 1888 PetscValidHeaderSpecific(F,MAT_CLASSID,1); 1889 ierr = PetscTryMethod(F,"MatMumpsCreateSchurComplement_C",(Mat,Mat*),(F,S));CHKERRQ(ierr); 1890 PetscFunctionReturn(0); 1891 } 1892 1893 /* -------------------------------------------------------------------------------------------*/ 1894 #undef __FUNCT__ 1895 #define __FUNCT__ "MatMumpsGetSchurComplement_MUMPS" 1896 PetscErrorCode MatMumpsGetSchurComplement_MUMPS(Mat F,Mat* S) 1897 { 1898 Mat St; 1899 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1900 PetscErrorCode ierr; 1901 1902 PetscFunctionBegin; 1903 if (!mumps->CleanUpMUMPS) { /* CleanUpMUMPS is set to true after numerical factorization */ 1904 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Numerical factorization phase not yet performed! You should call MatFactorSymbolic/Numeric before"); 1905 } else if (!mumps->id.ICNTL(19)) { 1906 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur complement mode not selected! You should call MatMumpsSetSchurIndices to enable it"); 1907 } else if (!mumps->id.size_schur) { 1908 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur indices not set! You should call MatMumpsSetSchurIndices before"); 1909 } else if (!mumps->schur_restored) { 1910 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur matrix has not been restored using MatMumpsRestoreSchurComplement"); 1911 } 1912 /* It should be the responsibility of the user to handle different ICNTL(19) cases if they want to work with the raw data */ 1913 /* should I also add errors when the Schur complement has been already factored? */ 1914 ierr = MatCreate(PetscObjectComm((PetscObject)F),&St);CHKERRQ(ierr); 1915 ierr = MatCreateSeqDense(PetscObjectComm((PetscObject)F),mumps->id.size_schur,mumps->id.size_schur,(PetscScalar*)mumps->id.schur,&St);CHKERRQ(ierr); 1916 *S = St; 1917 mumps->schur_restored = PETSC_FALSE; 1918 PetscFunctionReturn(0); 1919 } 1920 1921 #undef __FUNCT__ 1922 #define __FUNCT__ "MatMumpsGetSchurComplement" 1923 /*@ 1924 MatMumpsGetSchurComplement - Get a Schur complement matrix object using the current status of the raw Schur data computed by MUMPS during the factorization step 1925 1926 Logically Collective on Mat 1927 1928 Input Parameters: 1929 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1930 . *S - location where to return the Schur complement (MATDENSE) 1931 1932 Notes: 1933 MUMPS Schur complement mode is currently implemented for sequential matrices. 1934 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. 1935 1936 Level: advanced 1937 1938 References: MUMPS Users' Guide 1939 1940 .seealso: MatGetFactor(), MatMumpsSetSchurIndices(), MatMumpsRestoreSchurComplement(), MatMumpsCreateSchurComplement() 1941 @*/ 1942 PetscErrorCode MatMumpsGetSchurComplement(Mat F,Mat* S) 1943 { 1944 PetscErrorCode ierr; 1945 1946 PetscFunctionBegin; 1947 PetscValidHeaderSpecific(F,MAT_CLASSID,1); 1948 ierr = PetscUseMethod(F,"MatMumpsGetSchurComplement_C",(Mat,Mat*),(F,S));CHKERRQ(ierr); 1949 PetscFunctionReturn(0); 1950 } 1951 1952 /* -------------------------------------------------------------------------------------------*/ 1953 #undef __FUNCT__ 1954 #define __FUNCT__ "MatMumpsRestoreSchurComplement_MUMPS" 1955 PetscErrorCode MatMumpsRestoreSchurComplement_MUMPS(Mat F,Mat* S) 1956 { 1957 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 1958 PetscErrorCode ierr; 1959 1960 PetscFunctionBegin; 1961 if (!mumps->CleanUpMUMPS) { /* CleanUpMUMPS is set to true after numerical factorization */ 1962 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Numerical factorization phase not yet performed! You should call MatFactorSymbolic/Numeric before"); 1963 } else if (!mumps->id.ICNTL(19)) { 1964 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur complement mode not selected! You should call MatMumpsSetSchurIndices to enable it"); 1965 } else if (!mumps->id.size_schur) { 1966 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur indices not set! You should call MatMumpsSetSchurIndices before"); 1967 } else if (!mumps->schur_restored) { 1968 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur matrix has not been restored using MatMumpsRestoreSchurComplement"); 1969 } 1970 ierr = MatDestroy(S);CHKERRQ(ierr); 1971 *S = NULL; 1972 mumps->schur_restored = PETSC_TRUE; 1973 PetscFunctionReturn(0); 1974 } 1975 1976 #undef __FUNCT__ 1977 #define __FUNCT__ "MatMumpsRestoreSchurComplement" 1978 /*@ 1979 MatMumpsRestoreSchurComplement - Restore the Schur complement matrix object obtained from a call to MatGetSchurComplement 1980 1981 Logically Collective on Mat 1982 1983 Input Parameters: 1984 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 1985 . *S - location where the Schur complement is stored 1986 1987 Notes: 1988 MUMPS Schur complement mode is currently implemented for sequential matrices. 1989 1990 Level: advanced 1991 1992 References: MUMPS Users' Guide 1993 1994 .seealso: MatGetFactor(), MatMumpsSetSchurIndices(), MatMumpsGetSchurComplement(), MatMumpsCreateSchurComplement() 1995 @*/ 1996 PetscErrorCode MatMumpsRestoreSchurComplement(Mat F,Mat* S) 1997 { 1998 PetscErrorCode ierr; 1999 2000 PetscFunctionBegin; 2001 PetscValidHeaderSpecific(F,MAT_CLASSID,1); 2002 ierr = PetscUseMethod(F,"MatMumpsRestoreSchurComplement_C",(Mat,Mat*),(F,S));CHKERRQ(ierr); 2003 PetscFunctionReturn(0); 2004 } 2005 2006 /* -------------------------------------------------------------------------------------------*/ 2007 #undef __FUNCT__ 2008 #define __FUNCT__ "MatMumpsInvertSchurComplement_MUMPS" 2009 PetscErrorCode MatMumpsInvertSchurComplement_MUMPS(Mat F) 2010 { 2011 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2012 PetscErrorCode ierr; 2013 2014 PetscFunctionBegin; 2015 if (!mumps->id.ICNTL(19)) { /* do nothing */ 2016 PetscFunctionReturn(0); 2017 } 2018 if (!mumps->CleanUpMUMPS) { /* CleanUpMUMPS is set to true after numerical factorization */ 2019 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Numerical factorization phase not yet performed! You should call MatFactorSymbolic/Numeric before"); 2020 } else if (!mumps->id.size_schur) { 2021 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur indices not set! You should call MatMumpsSetSchurIndices before"); 2022 } else if (!mumps->schur_restored) { 2023 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur matrix has not been restored using MatMumpsRestoreSchurComplement"); 2024 } 2025 ierr = MatMumpsInvertSchur_Private(mumps);CHKERRQ(ierr); 2026 PetscFunctionReturn(0); 2027 } 2028 2029 #undef __FUNCT__ 2030 #define __FUNCT__ "MatMumpsInvertSchurComplement" 2031 /*@ 2032 MatMumpsInvertSchurComplement - Invert the raw Schur data computed by MUMPS during the factorization step 2033 2034 Logically Collective on Mat 2035 2036 Input Parameters: 2037 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2038 2039 Notes: 2040 MUMPS Schur complement mode is currently implemented for sequential matrices. 2041 The routine uses the pointer to the raw data of the Schur Complement stored within MUMPS data strutures. 2042 2043 Level: advanced 2044 2045 References: MUMPS Users' Guide 2046 2047 .seealso: MatGetFactor(), MatMumpsSetSchurIndices() 2048 @*/ 2049 PetscErrorCode MatMumpsInvertSchurComplement(Mat F) 2050 { 2051 PetscErrorCode ierr; 2052 2053 PetscFunctionBegin; 2054 PetscValidHeaderSpecific(F,MAT_CLASSID,1); 2055 ierr = PetscTryMethod(F,"MatMumpsInvertSchurComplement_C",(Mat),(F));CHKERRQ(ierr); 2056 PetscFunctionReturn(0); 2057 } 2058 2059 /* -------------------------------------------------------------------------------------------*/ 2060 #undef __FUNCT__ 2061 #define __FUNCT__ "MatMumpsSolveSchurComplement_MUMPS" 2062 PetscErrorCode MatMumpsSolveSchurComplement_MUMPS(Mat F, Vec rhs, Vec sol) 2063 { 2064 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2065 MumpsScalar *orhs; 2066 PetscScalar *osol,*nrhs,*nsol; 2067 PetscErrorCode ierr; 2068 2069 PetscFunctionBegin; 2070 if (!mumps->CleanUpMUMPS) { /* CleanUpMUMPS is set to true after numerical factorization */ 2071 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Numerical factorization phase not yet performed! You should call MatFactorSymbolic/Numeric before"); 2072 } else if (!mumps->id.ICNTL(19)) { 2073 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur complement mode not selected! You should call MatMumpsSetSchurIndices to enable it"); 2074 } else if (!mumps->id.size_schur) { 2075 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur indices not set! You should call MatMumpsSetSchurIndices before"); 2076 } else if (!mumps->schur_restored) { 2077 SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ORDER,"Schur matrix has not been restored using MatMumpsRestoreSchurComplement"); 2078 } 2079 /* swap pointers */ 2080 orhs = mumps->id.redrhs; 2081 osol = mumps->schur_sol; 2082 ierr = VecGetArray(rhs,&nrhs);CHKERRQ(ierr); 2083 ierr = VecGetArray(sol,&nsol);CHKERRQ(ierr); 2084 mumps->id.redrhs = (MumpsScalar*)nrhs; 2085 mumps->schur_sol = nsol; 2086 /* solve Schur complement */ 2087 mumps->id.nrhs = 1; 2088 ierr = MatMumpsSolveSchur_Private(mumps,PETSC_FALSE);CHKERRQ(ierr); 2089 /* restore pointers */ 2090 ierr = VecRestoreArray(rhs,&nrhs);CHKERRQ(ierr); 2091 ierr = VecRestoreArray(sol,&nsol);CHKERRQ(ierr); 2092 mumps->id.redrhs = orhs; 2093 mumps->schur_sol = osol; 2094 PetscFunctionReturn(0); 2095 } 2096 2097 #undef __FUNCT__ 2098 #define __FUNCT__ "MatMumpsSolveSchurComplement" 2099 /*@ 2100 MatMumpsSolveSchurComplement - Solve the Schur complement system computed by MUMPS during the factorization step 2101 2102 Logically Collective on Mat 2103 2104 Input Parameters: 2105 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2106 . rhs - location where the right hand side of the Schur complement system is stored 2107 - sol - location where the solution of the Schur complement system has to be returned 2108 2109 Notes: 2110 MUMPS Schur complement mode is currently implemented for sequential matrices. 2111 The sizes of the vectors should match the size of the Schur complement 2112 2113 Level: advanced 2114 2115 References: MUMPS Users' Guide 2116 2117 .seealso: MatGetFactor(), MatMumpsSetSchurIndices() 2118 @*/ 2119 PetscErrorCode MatMumpsSolveSchurComplement(Mat F, Vec rhs, Vec sol) 2120 { 2121 PetscErrorCode ierr; 2122 2123 PetscFunctionBegin; 2124 PetscValidHeaderSpecific(F,MAT_CLASSID,1); 2125 PetscValidHeaderSpecific(rhs,VEC_CLASSID,2); 2126 PetscValidHeaderSpecific(sol,VEC_CLASSID,2); 2127 PetscCheckSameComm(F,1,rhs,2); 2128 PetscCheckSameComm(F,1,sol,3); 2129 ierr = PetscTryMethod(F,"MatMumpsSolveSchurComplement_C",(Mat,Vec,Vec),(F,rhs,sol));CHKERRQ(ierr); 2130 PetscFunctionReturn(0); 2131 } 2132 2133 /* -------------------------------------------------------------------------------------------*/ 2134 #undef __FUNCT__ 2135 #define __FUNCT__ "MatMumpsSetIcntl_MUMPS" 2136 PetscErrorCode MatMumpsSetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt ival) 2137 { 2138 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2139 2140 PetscFunctionBegin; 2141 mumps->id.ICNTL(icntl) = ival; 2142 PetscFunctionReturn(0); 2143 } 2144 2145 #undef __FUNCT__ 2146 #define __FUNCT__ "MatMumpsGetIcntl_MUMPS" 2147 PetscErrorCode MatMumpsGetIcntl_MUMPS(Mat F,PetscInt icntl,PetscInt *ival) 2148 { 2149 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2150 2151 PetscFunctionBegin; 2152 *ival = mumps->id.ICNTL(icntl); 2153 PetscFunctionReturn(0); 2154 } 2155 2156 #undef __FUNCT__ 2157 #define __FUNCT__ "MatMumpsSetIcntl" 2158 /*@ 2159 MatMumpsSetIcntl - Set MUMPS parameter ICNTL() 2160 2161 Logically Collective on Mat 2162 2163 Input Parameters: 2164 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2165 . icntl - index of MUMPS parameter array ICNTL() 2166 - ival - value of MUMPS ICNTL(icntl) 2167 2168 Options Database: 2169 . -mat_mumps_icntl_<icntl> <ival> 2170 2171 Level: beginner 2172 2173 References: MUMPS Users' Guide 2174 2175 .seealso: MatGetFactor() 2176 @*/ 2177 PetscErrorCode MatMumpsSetIcntl(Mat F,PetscInt icntl,PetscInt ival) 2178 { 2179 PetscErrorCode ierr; 2180 2181 PetscFunctionBegin; 2182 PetscValidLogicalCollectiveInt(F,icntl,2); 2183 PetscValidLogicalCollectiveInt(F,ival,3); 2184 ierr = PetscTryMethod(F,"MatMumpsSetIcntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 2185 PetscFunctionReturn(0); 2186 } 2187 2188 #undef __FUNCT__ 2189 #define __FUNCT__ "MatMumpsGetIcntl" 2190 /*@ 2191 MatMumpsGetIcntl - Get MUMPS parameter ICNTL() 2192 2193 Logically Collective on Mat 2194 2195 Input Parameters: 2196 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2197 - icntl - index of MUMPS parameter array ICNTL() 2198 2199 Output Parameter: 2200 . ival - value of MUMPS ICNTL(icntl) 2201 2202 Level: beginner 2203 2204 References: MUMPS Users' Guide 2205 2206 .seealso: MatGetFactor() 2207 @*/ 2208 PetscErrorCode MatMumpsGetIcntl(Mat F,PetscInt icntl,PetscInt *ival) 2209 { 2210 PetscErrorCode ierr; 2211 2212 PetscFunctionBegin; 2213 PetscValidLogicalCollectiveInt(F,icntl,2); 2214 PetscValidIntPointer(ival,3); 2215 ierr = PetscTryMethod(F,"MatMumpsGetIcntl_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 2216 PetscFunctionReturn(0); 2217 } 2218 2219 /* -------------------------------------------------------------------------------------------*/ 2220 #undef __FUNCT__ 2221 #define __FUNCT__ "MatMumpsSetCntl_MUMPS" 2222 PetscErrorCode MatMumpsSetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal val) 2223 { 2224 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2225 2226 PetscFunctionBegin; 2227 mumps->id.CNTL(icntl) = val; 2228 PetscFunctionReturn(0); 2229 } 2230 2231 #undef __FUNCT__ 2232 #define __FUNCT__ "MatMumpsGetCntl_MUMPS" 2233 PetscErrorCode MatMumpsGetCntl_MUMPS(Mat F,PetscInt icntl,PetscReal *val) 2234 { 2235 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2236 2237 PetscFunctionBegin; 2238 *val = mumps->id.CNTL(icntl); 2239 PetscFunctionReturn(0); 2240 } 2241 2242 #undef __FUNCT__ 2243 #define __FUNCT__ "MatMumpsSetCntl" 2244 /*@ 2245 MatMumpsSetCntl - Set MUMPS parameter CNTL() 2246 2247 Logically Collective on Mat 2248 2249 Input Parameters: 2250 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2251 . icntl - index of MUMPS parameter array CNTL() 2252 - val - value of MUMPS CNTL(icntl) 2253 2254 Options Database: 2255 . -mat_mumps_cntl_<icntl> <val> 2256 2257 Level: beginner 2258 2259 References: MUMPS Users' Guide 2260 2261 .seealso: MatGetFactor() 2262 @*/ 2263 PetscErrorCode MatMumpsSetCntl(Mat F,PetscInt icntl,PetscReal val) 2264 { 2265 PetscErrorCode ierr; 2266 2267 PetscFunctionBegin; 2268 PetscValidLogicalCollectiveInt(F,icntl,2); 2269 PetscValidLogicalCollectiveReal(F,val,3); 2270 ierr = PetscTryMethod(F,"MatMumpsSetCntl_C",(Mat,PetscInt,PetscReal),(F,icntl,val));CHKERRQ(ierr); 2271 PetscFunctionReturn(0); 2272 } 2273 2274 #undef __FUNCT__ 2275 #define __FUNCT__ "MatMumpsGetCntl" 2276 /*@ 2277 MatMumpsGetCntl - Get MUMPS parameter CNTL() 2278 2279 Logically Collective on Mat 2280 2281 Input Parameters: 2282 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2283 - icntl - index of MUMPS parameter array CNTL() 2284 2285 Output Parameter: 2286 . val - value of MUMPS CNTL(icntl) 2287 2288 Level: beginner 2289 2290 References: MUMPS Users' Guide 2291 2292 .seealso: MatGetFactor() 2293 @*/ 2294 PetscErrorCode MatMumpsGetCntl(Mat F,PetscInt icntl,PetscReal *val) 2295 { 2296 PetscErrorCode ierr; 2297 2298 PetscFunctionBegin; 2299 PetscValidLogicalCollectiveInt(F,icntl,2); 2300 PetscValidRealPointer(val,3); 2301 ierr = PetscTryMethod(F,"MatMumpsGetCntl_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 2302 PetscFunctionReturn(0); 2303 } 2304 2305 #undef __FUNCT__ 2306 #define __FUNCT__ "MatMumpsGetInfo_MUMPS" 2307 PetscErrorCode MatMumpsGetInfo_MUMPS(Mat F,PetscInt icntl,PetscInt *info) 2308 { 2309 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2310 2311 PetscFunctionBegin; 2312 *info = mumps->id.INFO(icntl); 2313 PetscFunctionReturn(0); 2314 } 2315 2316 #undef __FUNCT__ 2317 #define __FUNCT__ "MatMumpsGetInfog_MUMPS" 2318 PetscErrorCode MatMumpsGetInfog_MUMPS(Mat F,PetscInt icntl,PetscInt *infog) 2319 { 2320 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2321 2322 PetscFunctionBegin; 2323 *infog = mumps->id.INFOG(icntl); 2324 PetscFunctionReturn(0); 2325 } 2326 2327 #undef __FUNCT__ 2328 #define __FUNCT__ "MatMumpsGetRinfo_MUMPS" 2329 PetscErrorCode MatMumpsGetRinfo_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfo) 2330 { 2331 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2332 2333 PetscFunctionBegin; 2334 *rinfo = mumps->id.RINFO(icntl); 2335 PetscFunctionReturn(0); 2336 } 2337 2338 #undef __FUNCT__ 2339 #define __FUNCT__ "MatMumpsGetRinfog_MUMPS" 2340 PetscErrorCode MatMumpsGetRinfog_MUMPS(Mat F,PetscInt icntl,PetscReal *rinfog) 2341 { 2342 Mat_MUMPS *mumps =(Mat_MUMPS*)F->spptr; 2343 2344 PetscFunctionBegin; 2345 *rinfog = mumps->id.RINFOG(icntl); 2346 PetscFunctionReturn(0); 2347 } 2348 2349 #undef __FUNCT__ 2350 #define __FUNCT__ "MatMumpsGetInfo" 2351 /*@ 2352 MatMumpsGetInfo - Get MUMPS parameter INFO() 2353 2354 Logically Collective on Mat 2355 2356 Input Parameters: 2357 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2358 - icntl - index of MUMPS parameter array INFO() 2359 2360 Output Parameter: 2361 . ival - value of MUMPS INFO(icntl) 2362 2363 Level: beginner 2364 2365 References: MUMPS Users' Guide 2366 2367 .seealso: MatGetFactor() 2368 @*/ 2369 PetscErrorCode MatMumpsGetInfo(Mat F,PetscInt icntl,PetscInt *ival) 2370 { 2371 PetscErrorCode ierr; 2372 2373 PetscFunctionBegin; 2374 PetscValidIntPointer(ival,3); 2375 ierr = PetscTryMethod(F,"MatMumpsGetInfo_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 2376 PetscFunctionReturn(0); 2377 } 2378 2379 #undef __FUNCT__ 2380 #define __FUNCT__ "MatMumpsGetInfog" 2381 /*@ 2382 MatMumpsGetInfog - Get MUMPS parameter INFOG() 2383 2384 Logically Collective on Mat 2385 2386 Input Parameters: 2387 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2388 - icntl - index of MUMPS parameter array INFOG() 2389 2390 Output Parameter: 2391 . ival - value of MUMPS INFOG(icntl) 2392 2393 Level: beginner 2394 2395 References: MUMPS Users' Guide 2396 2397 .seealso: MatGetFactor() 2398 @*/ 2399 PetscErrorCode MatMumpsGetInfog(Mat F,PetscInt icntl,PetscInt *ival) 2400 { 2401 PetscErrorCode ierr; 2402 2403 PetscFunctionBegin; 2404 PetscValidIntPointer(ival,3); 2405 ierr = PetscTryMethod(F,"MatMumpsGetInfog_C",(Mat,PetscInt,PetscInt*),(F,icntl,ival));CHKERRQ(ierr); 2406 PetscFunctionReturn(0); 2407 } 2408 2409 #undef __FUNCT__ 2410 #define __FUNCT__ "MatMumpsGetRinfo" 2411 /*@ 2412 MatMumpsGetRinfo - Get MUMPS parameter RINFO() 2413 2414 Logically Collective on Mat 2415 2416 Input Parameters: 2417 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2418 - icntl - index of MUMPS parameter array RINFO() 2419 2420 Output Parameter: 2421 . val - value of MUMPS RINFO(icntl) 2422 2423 Level: beginner 2424 2425 References: MUMPS Users' Guide 2426 2427 .seealso: MatGetFactor() 2428 @*/ 2429 PetscErrorCode MatMumpsGetRinfo(Mat F,PetscInt icntl,PetscReal *val) 2430 { 2431 PetscErrorCode ierr; 2432 2433 PetscFunctionBegin; 2434 PetscValidRealPointer(val,3); 2435 ierr = PetscTryMethod(F,"MatMumpsGetRinfo_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 2436 PetscFunctionReturn(0); 2437 } 2438 2439 #undef __FUNCT__ 2440 #define __FUNCT__ "MatMumpsGetRinfog" 2441 /*@ 2442 MatMumpsGetRinfog - Get MUMPS parameter RINFOG() 2443 2444 Logically Collective on Mat 2445 2446 Input Parameters: 2447 + F - the factored matrix obtained by calling MatGetFactor() from PETSc-MUMPS interface 2448 - icntl - index of MUMPS parameter array RINFOG() 2449 2450 Output Parameter: 2451 . val - value of MUMPS RINFOG(icntl) 2452 2453 Level: beginner 2454 2455 References: MUMPS Users' Guide 2456 2457 .seealso: MatGetFactor() 2458 @*/ 2459 PetscErrorCode MatMumpsGetRinfog(Mat F,PetscInt icntl,PetscReal *val) 2460 { 2461 PetscErrorCode ierr; 2462 2463 PetscFunctionBegin; 2464 PetscValidRealPointer(val,3); 2465 ierr = PetscTryMethod(F,"MatMumpsGetRinfog_C",(Mat,PetscInt,PetscReal*),(F,icntl,val));CHKERRQ(ierr); 2466 PetscFunctionReturn(0); 2467 } 2468 2469 /*MC 2470 MATSOLVERMUMPS - A matrix type providing direct solvers (LU and Cholesky) for 2471 distributed and sequential matrices via the external package MUMPS. 2472 2473 Works with MATAIJ and MATSBAIJ matrices 2474 2475 Options Database Keys: 2476 + -mat_mumps_icntl_1 <6>: ICNTL(1): output stream for error messages (None) 2477 . -mat_mumps_icntl_2 <0>: ICNTL(2): output stream for diagnostic printing, statistics, and warning (None) 2478 . -mat_mumps_icntl_3 <0>: ICNTL(3): output stream for global information, collected on the host (None) 2479 . -mat_mumps_icntl_4 <0>: ICNTL(4): level of printing (0 to 4) (None) 2480 . -mat_mumps_icntl_6 <7>: ICNTL(6): permutes to a zero-free diagonal and/or scale the matrix (0 to 7) (None) 2481 . -mat_mumps_icntl_7 <7>: ICNTL(7): computes a symmetric permutation in sequential analysis (0 to 7). 3=Scotch, 4=PORD, 5=Metis (None) 2482 . -mat_mumps_icntl_8 <77>: ICNTL(8): scaling strategy (-2 to 8 or 77) (None) 2483 . -mat_mumps_icntl_10 <0>: ICNTL(10): max num of refinements (None) 2484 . -mat_mumps_icntl_11 <0>: ICNTL(11): statistics related to an error analysis (via -ksp_view) (None) 2485 . -mat_mumps_icntl_12 <1>: ICNTL(12): an ordering strategy for symmetric matrices (0 to 3) (None) 2486 . -mat_mumps_icntl_13 <0>: ICNTL(13): parallelism of the root node (enable ScaLAPACK) and its splitting (None) 2487 . -mat_mumps_icntl_14 <20>: ICNTL(14): percentage increase in the estimated working space (None) 2488 . -mat_mumps_icntl_19 <0>: ICNTL(19): computes the Schur complement (None) 2489 . -mat_mumps_icntl_22 <0>: ICNTL(22): in-core/out-of-core factorization and solve (0 or 1) (None) 2490 . -mat_mumps_icntl_23 <0>: ICNTL(23): max size of the working memory (MB) that can allocate per processor (None) 2491 . -mat_mumps_icntl_24 <0>: ICNTL(24): detection of null pivot rows (0 or 1) (None) 2492 . -mat_mumps_icntl_25 <0>: ICNTL(25): compute a solution of a deficient matrix and a null space basis (None) 2493 . -mat_mumps_icntl_26 <0>: ICNTL(26): drives the solution phase if a Schur complement matrix (None) 2494 . -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) 2495 . -mat_mumps_icntl_29 <0>: ICNTL(29): parallel ordering 1 = ptscotch, 2 = parmetis (None) 2496 . -mat_mumps_icntl_30 <0>: ICNTL(30): compute user-specified set of entries in inv(A) (None) 2497 . -mat_mumps_icntl_31 <0>: ICNTL(31): indicates which factors may be discarded during factorization (None) 2498 . -mat_mumps_icntl_33 <0>: ICNTL(33): compute determinant (None) 2499 . -mat_mumps_cntl_1 <0.01>: CNTL(1): relative pivoting threshold (None) 2500 . -mat_mumps_cntl_2 <1.49012e-08>: CNTL(2): stopping criterion of refinement (None) 2501 . -mat_mumps_cntl_3 <0>: CNTL(3): absolute pivoting threshold (None) 2502 . -mat_mumps_cntl_4 <-1>: CNTL(4): value for static pivoting (None) 2503 - -mat_mumps_cntl_5 <0>: CNTL(5): fixation for null pivots (None) 2504 2505 Level: beginner 2506 2507 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 2508 2509 M*/ 2510 2511 #undef __FUNCT__ 2512 #define __FUNCT__ "MatFactorGetSolverPackage_mumps" 2513 static PetscErrorCode MatFactorGetSolverPackage_mumps(Mat A,const MatSolverPackage *type) 2514 { 2515 PetscFunctionBegin; 2516 *type = MATSOLVERMUMPS; 2517 PetscFunctionReturn(0); 2518 } 2519 2520 /* MatGetFactor for Seq and MPI AIJ matrices */ 2521 #undef __FUNCT__ 2522 #define __FUNCT__ "MatGetFactor_aij_mumps" 2523 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat A,MatFactorType ftype,Mat *F) 2524 { 2525 Mat B; 2526 PetscErrorCode ierr; 2527 Mat_MUMPS *mumps; 2528 PetscBool isSeqAIJ; 2529 2530 PetscFunctionBegin; 2531 /* Create the factorization matrix */ 2532 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 2533 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 2534 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 2535 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2536 if (isSeqAIJ) { 2537 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 2538 } else { 2539 ierr = MatMPIAIJSetPreallocation(B,0,NULL,0,NULL);CHKERRQ(ierr); 2540 } 2541 2542 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 2543 2544 B->ops->view = MatView_MUMPS; 2545 B->ops->getinfo = MatGetInfo_MUMPS; 2546 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 2547 2548 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 2549 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 2550 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 2551 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 2552 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 2553 2554 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 2555 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 2556 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 2557 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 2558 2559 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 2560 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsInvertSchurComplement_C",MatMumpsInvertSchurComplement_MUMPS);CHKERRQ(ierr); 2561 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsCreateSchurComplement_C",MatMumpsCreateSchurComplement_MUMPS);CHKERRQ(ierr); 2562 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 2563 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsRestoreSchurComplement_C",MatMumpsRestoreSchurComplement_MUMPS);CHKERRQ(ierr); 2564 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSolveSchurComplement_C",MatMumpsSolveSchurComplement_MUMPS);CHKERRQ(ierr); 2565 2566 if (ftype == MAT_FACTOR_LU) { 2567 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMUMPS; 2568 B->factortype = MAT_FACTOR_LU; 2569 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqaij; 2570 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij; 2571 mumps->sym = 0; 2572 } else { 2573 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 2574 B->factortype = MAT_FACTOR_CHOLESKY; 2575 if (isSeqAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqaij_seqsbaij; 2576 else mumps->ConvertToTriples = MatConvertToTriples_mpiaij_mpisbaij; 2577 #if defined(PETSC_USE_COMPLEX) 2578 mumps->sym = 2; 2579 #else 2580 if (A->spd_set && A->spd) mumps->sym = 1; 2581 else mumps->sym = 2; 2582 #endif 2583 } 2584 2585 mumps->isAIJ = PETSC_TRUE; 2586 mumps->Destroy = B->ops->destroy; 2587 B->ops->destroy = MatDestroy_MUMPS; 2588 B->spptr = (void*)mumps; 2589 2590 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 2591 2592 *F = B; 2593 PetscFunctionReturn(0); 2594 } 2595 2596 /* MatGetFactor for Seq and MPI SBAIJ matrices */ 2597 #undef __FUNCT__ 2598 #define __FUNCT__ "MatGetFactor_sbaij_mumps" 2599 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat A,MatFactorType ftype,Mat *F) 2600 { 2601 Mat B; 2602 PetscErrorCode ierr; 2603 Mat_MUMPS *mumps; 2604 PetscBool isSeqSBAIJ; 2605 2606 PetscFunctionBegin; 2607 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot use PETSc SBAIJ matrices with MUMPS LU, use AIJ matrix"); 2608 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"); 2609 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 2610 /* Create the factorization matrix */ 2611 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 2612 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 2613 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2614 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 2615 if (isSeqSBAIJ) { 2616 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 2617 2618 mumps->ConvertToTriples = MatConvertToTriples_seqsbaij_seqsbaij; 2619 } else { 2620 ierr = MatMPISBAIJSetPreallocation(B,1,0,NULL,0,NULL);CHKERRQ(ierr); 2621 2622 mumps->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij; 2623 } 2624 2625 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MUMPS; 2626 B->ops->view = MatView_MUMPS; 2627 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 2628 2629 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 2630 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 2631 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 2632 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 2633 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 2634 2635 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 2636 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 2637 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 2638 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 2639 2640 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 2641 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsInvertSchurComplement_C",MatMumpsInvertSchurComplement_MUMPS);CHKERRQ(ierr); 2642 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsCreateSchurComplement_C",MatMumpsCreateSchurComplement_MUMPS);CHKERRQ(ierr); 2643 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 2644 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsRestoreSchurComplement_C",MatMumpsRestoreSchurComplement_MUMPS);CHKERRQ(ierr); 2645 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSolveSchurComplement_C",MatMumpsSolveSchurComplement_MUMPS);CHKERRQ(ierr); 2646 2647 B->factortype = MAT_FACTOR_CHOLESKY; 2648 #if defined(PETSC_USE_COMPLEX) 2649 mumps->sym = 2; 2650 #else 2651 if (A->spd_set && A->spd) mumps->sym = 1; 2652 else mumps->sym = 2; 2653 #endif 2654 2655 mumps->isAIJ = PETSC_FALSE; 2656 mumps->Destroy = B->ops->destroy; 2657 B->ops->destroy = MatDestroy_MUMPS; 2658 B->spptr = (void*)mumps; 2659 2660 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 2661 2662 *F = B; 2663 PetscFunctionReturn(0); 2664 } 2665 2666 #undef __FUNCT__ 2667 #define __FUNCT__ "MatGetFactor_baij_mumps" 2668 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat A,MatFactorType ftype,Mat *F) 2669 { 2670 Mat B; 2671 PetscErrorCode ierr; 2672 Mat_MUMPS *mumps; 2673 PetscBool isSeqBAIJ; 2674 2675 PetscFunctionBegin; 2676 /* Create the factorization matrix */ 2677 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 2678 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 2679 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 2680 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2681 if (isSeqBAIJ) { 2682 ierr = MatSeqBAIJSetPreallocation(B,A->rmap->bs,0,NULL);CHKERRQ(ierr); 2683 } else { 2684 ierr = MatMPIBAIJSetPreallocation(B,A->rmap->bs,0,NULL,0,NULL);CHKERRQ(ierr); 2685 } 2686 2687 ierr = PetscNewLog(B,&mumps);CHKERRQ(ierr); 2688 if (ftype == MAT_FACTOR_LU) { 2689 B->ops->lufactorsymbolic = MatLUFactorSymbolic_BAIJMUMPS; 2690 B->factortype = MAT_FACTOR_LU; 2691 if (isSeqBAIJ) mumps->ConvertToTriples = MatConvertToTriples_seqbaij_seqaij; 2692 else mumps->ConvertToTriples = MatConvertToTriples_mpibaij_mpiaij; 2693 mumps->sym = 0; 2694 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot use PETSc BAIJ matrices with MUMPS Cholesky, use SBAIJ or AIJ matrix instead\n"); 2695 2696 B->ops->view = MatView_MUMPS; 2697 B->ops->getdiagonal = MatGetDiagonal_MUMPS; 2698 2699 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mumps);CHKERRQ(ierr); 2700 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetIcntl_C",MatMumpsSetIcntl_MUMPS);CHKERRQ(ierr); 2701 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetIcntl_C",MatMumpsGetIcntl_MUMPS);CHKERRQ(ierr); 2702 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetCntl_C",MatMumpsSetCntl_MUMPS);CHKERRQ(ierr); 2703 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetCntl_C",MatMumpsGetCntl_MUMPS);CHKERRQ(ierr); 2704 2705 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfo_C",MatMumpsGetInfo_MUMPS);CHKERRQ(ierr); 2706 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetInfog_C",MatMumpsGetInfog_MUMPS);CHKERRQ(ierr); 2707 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfo_C",MatMumpsGetRinfo_MUMPS);CHKERRQ(ierr); 2708 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetRinfog_C",MatMumpsGetRinfog_MUMPS);CHKERRQ(ierr); 2709 2710 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSetSchurIndices_C",MatMumpsSetSchurIndices_MUMPS);CHKERRQ(ierr); 2711 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsInvertSchurComplement_C",MatMumpsInvertSchurComplement_MUMPS);CHKERRQ(ierr); 2712 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsCreateSchurComplement_C",MatMumpsCreateSchurComplement_MUMPS);CHKERRQ(ierr); 2713 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsGetSchurComplement_C",MatMumpsGetSchurComplement_MUMPS);CHKERRQ(ierr); 2714 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsRestoreSchurComplement_C",MatMumpsRestoreSchurComplement_MUMPS);CHKERRQ(ierr); 2715 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMumpsSolveSchurComplement_C",MatMumpsSolveSchurComplement_MUMPS);CHKERRQ(ierr); 2716 2717 mumps->isAIJ = PETSC_TRUE; 2718 mumps->Destroy = B->ops->destroy; 2719 B->ops->destroy = MatDestroy_MUMPS; 2720 B->spptr = (void*)mumps; 2721 2722 ierr = PetscInitializeMUMPS(A,mumps);CHKERRQ(ierr); 2723 2724 *F = B; 2725 PetscFunctionReturn(0); 2726 } 2727 2728 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 2729 PETSC_EXTERN PetscErrorCode MatGetFactor_baij_mumps(Mat,MatFactorType,Mat*); 2730 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 2731 2732 #undef __FUNCT__ 2733 #define __FUNCT__ "MatSolverPackageRegister_MUMPS" 2734 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MUMPS(void) 2735 { 2736 PetscErrorCode ierr; 2737 2738 PetscFunctionBegin; 2739 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2740 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2741 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ,MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2742 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPIBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2743 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATMPISBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2744 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2745 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mumps);CHKERRQ(ierr); 2746 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ,MAT_FACTOR_LU,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2747 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_baij_mumps);CHKERRQ(ierr); 2748 ierr = MatSolverPackageRegister(MATSOLVERMUMPS,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 2749 PetscFunctionReturn(0); 2750 } 2751 2752