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