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