1 #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64) 2 #define MKL_ILP64 3 #endif 4 5 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 6 #include <../src/mat/impls/dense/seq/dense.h> 7 8 #include <stdio.h> 9 #include <stdlib.h> 10 #include <math.h> 11 #include <mkl.h> 12 13 /* 14 * Possible mkl_pardiso phases that controls the execution of the solver. 15 * For more information check mkl_pardiso manual. 16 */ 17 #define JOB_ANALYSIS 11 18 #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12 19 #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13 20 #define JOB_NUMERICAL_FACTORIZATION 22 21 #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23 22 #define JOB_SOLVE_ITERATIVE_REFINEMENT 33 23 #define JOB_SOLVE_FORWARD_SUBSTITUTION 331 24 #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332 25 #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333 26 #define JOB_RELEASE_OF_LU_MEMORY 0 27 #define JOB_RELEASE_OF_ALL_MEMORY -1 28 29 #define IPARM_SIZE 64 30 31 #if defined(PETSC_USE_64BIT_INDICES) 32 #if defined(PETSC_HAVE_LIBMKL_INTEL_ILP64) 33 /* sizeof(MKL_INT) == sizeof(long long int) if ilp64*/ 34 #define INT_TYPE long long int 35 #define MKL_PARDISO pardiso 36 #define MKL_PARDISO_INIT pardisoinit 37 #else 38 #define INT_TYPE long long int 39 #define MKL_PARDISO pardiso_64 40 #define MKL_PARDISO_INIT pardiso_64init 41 #endif 42 #else 43 #define INT_TYPE int 44 #define MKL_PARDISO pardiso 45 #define MKL_PARDISO_INIT pardisoinit 46 #endif 47 48 49 /* 50 * Internal data structure. 51 * For more information check mkl_pardiso manual. 52 */ 53 typedef struct { 54 55 /* Configuration vector*/ 56 INT_TYPE iparm[IPARM_SIZE]; 57 58 /* 59 * Internal mkl_pardiso memory location. 60 * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak. 61 */ 62 void *pt[IPARM_SIZE]; 63 64 /* Basic mkl_pardiso info*/ 65 INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err; 66 67 /* Matrix structure*/ 68 void *a; 69 INT_TYPE *ia, *ja; 70 71 /* Number of non-zero elements*/ 72 INT_TYPE nz; 73 74 /* Row permutaton vector*/ 75 INT_TYPE *perm; 76 77 /* Define if matrix preserves sparse structure.*/ 78 MatStructure matstruc; 79 80 /* True if mkl_pardiso function have been used.*/ 81 PetscBool CleanUp; 82 } Mat_MKL_PARDISO; 83 84 85 void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm []) 86 { 87 int iparm_copy[IPARM_SIZE], mtype_copy, i; 88 89 mtype_copy = *mtype; 90 pardisoinit(pt, &mtype_copy, iparm_copy); 91 for(i = 0; i < IPARM_SIZE; i++){ 92 iparm[i] = iparm_copy[i]; 93 } 94 } 95 96 97 /* 98 * Copy the elements of matrix A. 99 * Input: 100 * - Mat A: MATSEQAIJ matrix 101 * - int shift: matrix index. 102 * - 0 for c representation 103 * - 1 for fortran representation 104 * - MatReuse reuse: 105 * - MAT_INITIAL_MATRIX: Create a new aij representation 106 * - MAT_REUSE_MATRIX: Reuse all aij representation and just change values 107 * Output: 108 * - int *nnz: Number of nonzero-elements. 109 * - int **r pointer to i index 110 * - int **c pointer to j elements 111 * - MATRIXTYPE **v: Non-zero elements 112 */ 113 #undef __FUNCT__ 114 #define __FUNCT__ "MatCopy_MKL_PARDISO" 115 PetscErrorCode MatCopy_MKL_PARDISO(Mat A, MatReuse reuse, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, void **v) 116 { 117 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 118 119 PetscFunctionBegin; 120 *v=aa->a; 121 if (reuse == MAT_INITIAL_MATRIX) { 122 *r = (INT_TYPE*)aa->i; 123 *c = (INT_TYPE*)aa->j; 124 *nnz = aa->nz; 125 } 126 PetscFunctionReturn(0); 127 } 128 129 130 /* 131 * Free memory for Mat_MKL_PARDISO structure and pointers to objects. 132 */ 133 #undef __FUNCT__ 134 #define __FUNCT__ "MatDestroy_MKL_PARDISO" 135 PetscErrorCode MatDestroy_MKL_PARDISO(Mat A) 136 { 137 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr; 138 PetscErrorCode ierr; 139 140 PetscFunctionBegin; 141 /* Terminate instance, deallocate memories */ 142 if (mat_mkl_pardiso->CleanUp) { 143 mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY; 144 145 MKL_PARDISO (mat_mkl_pardiso->pt, 146 &mat_mkl_pardiso->maxfct, 147 &mat_mkl_pardiso->mnum, 148 &mat_mkl_pardiso->mtype, 149 &mat_mkl_pardiso->phase, 150 &mat_mkl_pardiso->n, 151 NULL, 152 NULL, 153 NULL, 154 mat_mkl_pardiso->perm, 155 &mat_mkl_pardiso->nrhs, 156 mat_mkl_pardiso->iparm, 157 &mat_mkl_pardiso->msglvl, 158 NULL, 159 NULL, 160 &mat_mkl_pardiso->err); 161 } 162 ierr = PetscFree(mat_mkl_pardiso->perm);CHKERRQ(ierr); 163 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 164 165 /* clear composed functions */ 166 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);CHKERRQ(ierr); 167 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);CHKERRQ(ierr); 168 PetscFunctionReturn(0); 169 } 170 171 /* 172 * Computes Ax = b 173 */ 174 #undef __FUNCT__ 175 #define __FUNCT__ "MatSolve_MKL_PARDISO" 176 PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x) 177 { 178 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr; 179 PetscErrorCode ierr; 180 PetscScalar *xarray; 181 const PetscScalar *barray; 182 183 PetscFunctionBegin; 184 mat_mkl_pardiso->nrhs = 1; 185 ierr = VecGetArray(x,&xarray);CHKERRQ(ierr); 186 ierr = VecGetArrayRead(b,&barray);CHKERRQ(ierr); 187 188 /* solve phase */ 189 /*-------------*/ 190 mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 191 MKL_PARDISO (mat_mkl_pardiso->pt, 192 &mat_mkl_pardiso->maxfct, 193 &mat_mkl_pardiso->mnum, 194 &mat_mkl_pardiso->mtype, 195 &mat_mkl_pardiso->phase, 196 &mat_mkl_pardiso->n, 197 mat_mkl_pardiso->a, 198 mat_mkl_pardiso->ia, 199 mat_mkl_pardiso->ja, 200 mat_mkl_pardiso->perm, 201 &mat_mkl_pardiso->nrhs, 202 mat_mkl_pardiso->iparm, 203 &mat_mkl_pardiso->msglvl, 204 (void*)barray, 205 (void*)xarray, 206 &mat_mkl_pardiso->err); 207 208 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err); 209 ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr); 210 ierr = VecRestoreArrayRead(b,&barray);CHKERRQ(ierr); 211 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 212 PetscFunctionReturn(0); 213 } 214 215 216 #undef __FUNCT__ 217 #define __FUNCT__ "MatSolveTranspose_MKL_PARDISO" 218 PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x) 219 { 220 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr; 221 PetscErrorCode ierr; 222 223 PetscFunctionBegin; 224 #if defined(PETSC_USE_COMPLEX) 225 mat_mkl_pardiso->iparm[12 - 1] = 1; 226 #else 227 mat_mkl_pardiso->iparm[12 - 1] = 2; 228 #endif 229 ierr = MatSolve_MKL_PARDISO(A,b,x);CHKERRQ(ierr); 230 mat_mkl_pardiso->iparm[12 - 1] = 0; 231 PetscFunctionReturn(0); 232 } 233 234 235 #undef __FUNCT__ 236 #define __FUNCT__ "MatMatSolve_MKL_PARDISO" 237 PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X) 238 { 239 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr; 240 PetscErrorCode ierr; 241 PetscScalar *barray, *xarray; 242 PetscBool flg; 243 244 PetscFunctionBegin; 245 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 246 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix"); 247 ierr = PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&flg);CHKERRQ(ierr); 248 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix"); 249 250 ierr = MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);CHKERRQ(ierr); 251 252 if(mat_mkl_pardiso->nrhs > 0){ 253 ierr = MatDenseGetArray(B,&barray); 254 ierr = MatDenseGetArray(X,&xarray); 255 256 /* solve phase */ 257 /*-------------*/ 258 mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 259 MKL_PARDISO (mat_mkl_pardiso->pt, 260 &mat_mkl_pardiso->maxfct, 261 &mat_mkl_pardiso->mnum, 262 &mat_mkl_pardiso->mtype, 263 &mat_mkl_pardiso->phase, 264 &mat_mkl_pardiso->n, 265 mat_mkl_pardiso->a, 266 mat_mkl_pardiso->ia, 267 mat_mkl_pardiso->ja, 268 mat_mkl_pardiso->perm, 269 &mat_mkl_pardiso->nrhs, 270 mat_mkl_pardiso->iparm, 271 &mat_mkl_pardiso->msglvl, 272 (void*)barray, 273 (void*)xarray, 274 &mat_mkl_pardiso->err); 275 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err); 276 } 277 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 278 PetscFunctionReturn(0); 279 } 280 281 /* 282 * LU Decomposition 283 */ 284 #undef __FUNCT__ 285 #define __FUNCT__ "MatFactorNumeric_MKL_PARDISO" 286 PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info) 287 { 288 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->spptr; 289 PetscErrorCode ierr; 290 291 /* numerical factorization phase */ 292 /*-------------------------------*/ 293 PetscFunctionBegin; 294 mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN; 295 ierr = MatCopy_MKL_PARDISO(A, MAT_REUSE_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);CHKERRQ(ierr); 296 297 /* numerical factorization phase */ 298 /*-------------------------------*/ 299 mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION; 300 MKL_PARDISO (mat_mkl_pardiso->pt, 301 &mat_mkl_pardiso->maxfct, 302 &mat_mkl_pardiso->mnum, 303 &mat_mkl_pardiso->mtype, 304 &mat_mkl_pardiso->phase, 305 &mat_mkl_pardiso->n, 306 mat_mkl_pardiso->a, 307 mat_mkl_pardiso->ia, 308 mat_mkl_pardiso->ja, 309 mat_mkl_pardiso->perm, 310 &mat_mkl_pardiso->nrhs, 311 mat_mkl_pardiso->iparm, 312 &mat_mkl_pardiso->msglvl, 313 NULL, 314 NULL, 315 &mat_mkl_pardiso->err); 316 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual\n",mat_mkl_pardiso->err); 317 318 mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN; 319 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 320 PetscFunctionReturn(0); 321 } 322 323 /* Sets mkl_pardiso options from the options database */ 324 #undef __FUNCT__ 325 #define __FUNCT__ "PetscSetMKL_PARDISOFromOptions" 326 PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A) 327 { 328 Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr; 329 PetscErrorCode ierr; 330 PetscInt icntl; 331 PetscBool flg; 332 int pt[IPARM_SIZE], threads = 1; 333 334 PetscFunctionBegin; 335 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");CHKERRQ(ierr); 336 ierr = PetscOptionsInt("-mat_mkl_pardiso_65","Number of threads to use","None",threads,&threads,&flg);CHKERRQ(ierr); 337 if (flg) mkl_set_num_threads(threads); 338 339 ierr = PetscOptionsInt("-mat_mkl_pardiso_66","Maximum number of factors with identical sparsity structure that must be kept in memory at the same time","None",mat_mkl_pardiso->maxfct,&icntl,&flg);CHKERRQ(ierr); 340 if (flg) mat_mkl_pardiso->maxfct = icntl; 341 342 ierr = PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);CHKERRQ(ierr); 343 if (flg) mat_mkl_pardiso->mnum = icntl; 344 345 ierr = PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);CHKERRQ(ierr); 346 if (flg) mat_mkl_pardiso->msglvl = icntl; 347 348 ierr = PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);CHKERRQ(ierr); 349 if(flg){ 350 mat_mkl_pardiso->mtype = icntl; 351 MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm); 352 #if defined(PETSC_USE_REAL_SINGLE) 353 mat_mkl_pardiso->iparm[27] = 1; 354 #else 355 mat_mkl_pardiso->iparm[27] = 0; 356 #endif 357 mat_mkl_pardiso->iparm[34] = 1; 358 } 359 ierr = PetscOptionsInt("-mat_mkl_pardiso_1","Use default values","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);CHKERRQ(ierr); 360 361 if(flg && icntl != 0){ 362 ierr = PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);CHKERRQ(ierr); 363 if (flg) mat_mkl_pardiso->iparm[1] = icntl; 364 365 ierr = PetscOptionsInt("-mat_mkl_pardiso_4","Preconditioned CGS/CG","None",mat_mkl_pardiso->iparm[3],&icntl,&flg);CHKERRQ(ierr); 366 if (flg) mat_mkl_pardiso->iparm[3] = icntl; 367 368 ierr = PetscOptionsInt("-mat_mkl_pardiso_5","User permutation","None",mat_mkl_pardiso->iparm[4],&icntl,&flg);CHKERRQ(ierr); 369 if (flg) mat_mkl_pardiso->iparm[4] = icntl; 370 371 ierr = PetscOptionsInt("-mat_mkl_pardiso_6","Write solution on x","None",mat_mkl_pardiso->iparm[5],&icntl,&flg);CHKERRQ(ierr); 372 if (flg) mat_mkl_pardiso->iparm[5] = icntl; 373 374 ierr = PetscOptionsInt("-mat_mkl_pardiso_8","Iterative refinement step","None",mat_mkl_pardiso->iparm[7],&icntl,&flg);CHKERRQ(ierr); 375 if (flg) mat_mkl_pardiso->iparm[7] = icntl; 376 377 ierr = PetscOptionsInt("-mat_mkl_pardiso_10","Pivoting perturbation","None",mat_mkl_pardiso->iparm[9],&icntl,&flg);CHKERRQ(ierr); 378 if (flg) mat_mkl_pardiso->iparm[9] = icntl; 379 380 ierr = PetscOptionsInt("-mat_mkl_pardiso_11","Scaling vectors","None",mat_mkl_pardiso->iparm[10],&icntl,&flg);CHKERRQ(ierr); 381 if (flg) mat_mkl_pardiso->iparm[10] = icntl; 382 383 ierr = PetscOptionsInt("-mat_mkl_pardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_pardiso->iparm[11],&icntl,&flg);CHKERRQ(ierr); 384 if (flg) mat_mkl_pardiso->iparm[11] = icntl; 385 386 ierr = PetscOptionsInt("-mat_mkl_pardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_pardiso->iparm[12],&icntl,&flg);CHKERRQ(ierr); 387 if (flg) mat_mkl_pardiso->iparm[12] = icntl; 388 389 ierr = PetscOptionsInt("-mat_mkl_pardiso_18","Numbers of non-zero elements","None",mat_mkl_pardiso->iparm[17],&icntl,&flg);CHKERRQ(ierr); 390 if (flg) mat_mkl_pardiso->iparm[17] = icntl; 391 392 ierr = PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);CHKERRQ(ierr); 393 if (flg) mat_mkl_pardiso->iparm[18] = icntl; 394 395 ierr = PetscOptionsInt("-mat_mkl_pardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_pardiso->iparm[20],&icntl,&flg);CHKERRQ(ierr); 396 if (flg) mat_mkl_pardiso->iparm[20] = icntl; 397 398 ierr = PetscOptionsInt("-mat_mkl_pardiso_24","Parallel factorization control","None",mat_mkl_pardiso->iparm[23],&icntl,&flg);CHKERRQ(ierr); 399 if (flg) mat_mkl_pardiso->iparm[23] = icntl; 400 401 ierr = PetscOptionsInt("-mat_mkl_pardiso_25","Parallel forward/backward solve control","None",mat_mkl_pardiso->iparm[24],&icntl,&flg);CHKERRQ(ierr); 402 if (flg) mat_mkl_pardiso->iparm[24] = icntl; 403 404 ierr = PetscOptionsInt("-mat_mkl_pardiso_27","Matrix checker","None",mat_mkl_pardiso->iparm[26],&icntl,&flg);CHKERRQ(ierr); 405 if (flg) mat_mkl_pardiso->iparm[26] = icntl; 406 407 ierr = PetscOptionsInt("-mat_mkl_pardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_pardiso->iparm[30],&icntl,&flg);CHKERRQ(ierr); 408 if (flg) mat_mkl_pardiso->iparm[30] = icntl; 409 410 ierr = PetscOptionsInt("-mat_mkl_pardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_pardiso->iparm[33],&icntl,&flg);CHKERRQ(ierr); 411 if (flg) mat_mkl_pardiso->iparm[33] = icntl; 412 413 ierr = PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);CHKERRQ(ierr); 414 if (flg) mat_mkl_pardiso->iparm[59] = icntl; 415 } 416 PetscOptionsEnd(); 417 PetscFunctionReturn(0); 418 } 419 420 #undef __FUNCT__ 421 #define __FUNCT__ "MatFactorMKL_PARDISOInitialize_Private" 422 PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso) 423 { 424 PetscErrorCode ierr; 425 PetscInt i; 426 427 PetscFunctionBegin; 428 for ( i = 0; i < IPARM_SIZE; i++ ){ 429 mat_mkl_pardiso->iparm[i] = 0; 430 } 431 432 for ( i = 0; i < IPARM_SIZE; i++ ){ 433 mat_mkl_pardiso->pt[i] = 0; 434 } 435 436 /*Default options for both sym and unsym */ 437 mat_mkl_pardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */ 438 mat_mkl_pardiso->iparm[ 1] = 2; /* Metis reordering */ 439 mat_mkl_pardiso->iparm[ 5] = 0; /* Write solution into x */ 440 mat_mkl_pardiso->iparm[ 7] = 2; /* Max number of iterative refinement steps */ 441 mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */ 442 mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */ 443 #if 0 444 mat_mkl_pardiso->iparm[23] = 1; /* Parallel factorization control*/ 445 #endif 446 mat_mkl_pardiso->iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */ 447 mat_mkl_pardiso->iparm[39] = 0; /* Input: matrix/rhs/solution stored on master */ 448 449 mat_mkl_pardiso->CleanUp = PETSC_FALSE; 450 mat_mkl_pardiso->maxfct = 1; /* Maximum number of numerical factorizations. */ 451 mat_mkl_pardiso->mnum = 1; /* Which factorization to use. */ 452 mat_mkl_pardiso->msglvl = 0; /* 0: do not print 1: Print statistical information in file */ 453 mat_mkl_pardiso->phase = -1; 454 mat_mkl_pardiso->err = 0; 455 456 mat_mkl_pardiso->n = A->rmap->N; 457 mat_mkl_pardiso->nrhs = 1; 458 mat_mkl_pardiso->err = 0; 459 mat_mkl_pardiso->phase = -1; 460 461 if(ftype == MAT_FACTOR_LU){ 462 /*Default type for non-sym*/ 463 #if defined(PETSC_USE_COMPLEX) 464 mat_mkl_pardiso->mtype = 13; 465 #else 466 mat_mkl_pardiso->mtype = 11; 467 #endif 468 469 mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */ 470 mat_mkl_pardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ 471 mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */ 472 473 } else { 474 /*Default type for sym*/ 475 #if defined(PETSC_USE_COMPLEX) 476 mat_mkl_pardiso ->mtype = 3; 477 #else 478 mat_mkl_pardiso ->mtype = -2; 479 #endif 480 mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */ 481 mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */ 482 mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */ 483 /* mat_mkl_pardiso->iparm[20] = 1; */ /* Apply 1x1 and 2x2 Bunch-Kaufman pivoting during the factorization process */ 484 #if defined(PETSC_USE_DEBUG) 485 mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */ 486 #endif 487 } 488 ierr = PetscMalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);CHKERRQ(ierr); 489 for(i = 0; i < A->rmap->N; i++){ 490 mat_mkl_pardiso->perm[i] = 0; 491 } 492 PetscFunctionReturn(0); 493 } 494 495 /* 496 * Symbolic decomposition. Mkl_Pardiso analysis phase. 497 */ 498 #undef __FUNCT__ 499 #define __FUNCT__ "MatFactorSymbolic_AIJMKL_PARDISO_Private" 500 PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info) 501 { 502 Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr; 503 PetscErrorCode ierr; 504 505 PetscFunctionBegin; 506 mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN; 507 508 /* Set MKL_PARDISO options from the options database */ 509 ierr = PetscSetMKL_PARDISOFromOptions(F,A);CHKERRQ(ierr); 510 511 ierr = MatCopy_MKL_PARDISO(A, MAT_INITIAL_MATRIX, &mat_mkl_pardiso->nz, &mat_mkl_pardiso->ia, &mat_mkl_pardiso->ja, &mat_mkl_pardiso->a);CHKERRQ(ierr); 512 mat_mkl_pardiso->n = A->rmap->N; 513 514 /* analysis phase */ 515 /*----------------*/ 516 mat_mkl_pardiso->phase = JOB_ANALYSIS; 517 518 MKL_PARDISO (mat_mkl_pardiso->pt, 519 &mat_mkl_pardiso->maxfct, 520 &mat_mkl_pardiso->mnum, 521 &mat_mkl_pardiso->mtype, 522 &mat_mkl_pardiso->phase, 523 &mat_mkl_pardiso->n, 524 mat_mkl_pardiso->a, 525 mat_mkl_pardiso->ia, 526 mat_mkl_pardiso->ja, 527 mat_mkl_pardiso->perm, 528 &mat_mkl_pardiso->nrhs, 529 mat_mkl_pardiso->iparm, 530 &mat_mkl_pardiso->msglvl, 531 NULL, 532 NULL, 533 &mat_mkl_pardiso->err); 534 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d\n. Please check manual",mat_mkl_pardiso->err); 535 536 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 537 538 if(F->factortype == MAT_FACTOR_LU){ 539 F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO; 540 } else { 541 F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO; 542 } 543 F->ops->solve = MatSolve_MKL_PARDISO; 544 F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO; 545 F->ops->matsolve = MatMatSolve_MKL_PARDISO; 546 PetscFunctionReturn(0); 547 } 548 549 #undef __FUNCT__ 550 #define __FUNCT__ "MatLUFactorSymbolic_AIJMKL_PARDISO" 551 PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 552 { 553 PetscErrorCode ierr; 554 555 PetscFunctionBegin; 556 ierr = MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);CHKERRQ(ierr); 557 PetscFunctionReturn(0); 558 } 559 560 #undef __FUNCT__ 561 #define __FUNCT__ "MatCholeskyFactorSymbolic_AIJMKL_PARDISO" 562 PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,const MatFactorInfo *info) 563 { 564 PetscErrorCode ierr; 565 566 PetscFunctionBegin; 567 ierr = MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);CHKERRQ(ierr); 568 PetscFunctionReturn(0); 569 } 570 571 #undef __FUNCT__ 572 #define __FUNCT__ "MatView_MKL_PARDISO" 573 PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer) 574 { 575 PetscErrorCode ierr; 576 PetscBool iascii; 577 PetscViewerFormat format; 578 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr; 579 PetscInt i; 580 581 PetscFunctionBegin; 582 /* check if matrix is mkl_pardiso type */ 583 if (A->ops->solve != MatSolve_MKL_PARDISO) PetscFunctionReturn(0); 584 585 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 586 if (iascii) { 587 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 588 if (format == PETSC_VIEWER_ASCII_INFO) { 589 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");CHKERRQ(ierr); 590 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase: %d \n",mat_mkl_pardiso->phase);CHKERRQ(ierr); 591 for(i = 1; i <= 64; i++){ 592 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]: %d \n",i, mat_mkl_pardiso->iparm[i - 1]);CHKERRQ(ierr); 593 } 594 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct);CHKERRQ(ierr); 595 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum: %d \n", mat_mkl_pardiso->mnum);CHKERRQ(ierr); 596 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype: %d \n", mat_mkl_pardiso->mtype);CHKERRQ(ierr); 597 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n: %d \n", mat_mkl_pardiso->n);CHKERRQ(ierr); 598 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs);CHKERRQ(ierr); 599 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl);CHKERRQ(ierr); 600 } 601 } 602 PetscFunctionReturn(0); 603 } 604 605 606 #undef __FUNCT__ 607 #define __FUNCT__ "MatGetInfo_MKL_PARDISO" 608 PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info) 609 { 610 Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)A->spptr; 611 612 PetscFunctionBegin; 613 info->block_size = 1.0; 614 info->nz_allocated = mat_mkl_pardiso->nz + 0.0; 615 info->nz_unneeded = 0.0; 616 info->assemblies = 0.0; 617 info->mallocs = 0.0; 618 info->memory = 0.0; 619 info->fill_ratio_given = 0; 620 info->fill_ratio_needed = 0; 621 info->factor_mallocs = 0; 622 PetscFunctionReturn(0); 623 } 624 625 #undef __FUNCT__ 626 #define __FUNCT__ "MatMkl_PardisoSetCntl_MKL_PARDISO" 627 PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival) 628 { 629 Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)F->spptr; 630 631 PetscFunctionBegin; 632 if(icntl <= 64){ 633 mat_mkl_pardiso->iparm[icntl - 1] = ival; 634 } else { 635 if(icntl == 65) 636 mkl_set_num_threads((int)ival); 637 else if(icntl == 66) 638 mat_mkl_pardiso->maxfct = ival; 639 else if(icntl == 67) 640 mat_mkl_pardiso->mnum = ival; 641 else if(icntl == 68) 642 mat_mkl_pardiso->msglvl = ival; 643 else if(icntl == 69){ 644 int pt[IPARM_SIZE]; 645 mat_mkl_pardiso->mtype = ival; 646 MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm); 647 #if defined(PETSC_USE_REAL_SINGLE) 648 mat_mkl_pardiso->iparm[27] = 1; 649 #else 650 mat_mkl_pardiso->iparm[27] = 0; 651 #endif 652 mat_mkl_pardiso->iparm[34] = 1; 653 } 654 } 655 PetscFunctionReturn(0); 656 } 657 658 #undef __FUNCT__ 659 #define __FUNCT__ "MatMkl_PardisoSetCntl" 660 /*@ 661 MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters 662 663 Logically Collective on Mat 664 665 Input Parameters: 666 + F - the factored matrix obtained by calling MatGetFactor() 667 . icntl - index of Mkl_Pardiso parameter 668 - ival - value of Mkl_Pardiso parameter 669 670 Options Database: 671 . -mat_mkl_pardiso_<icntl> <ival> 672 673 Level: beginner 674 675 References: Mkl_Pardiso Users' Guide 676 677 .seealso: MatGetFactor() 678 @*/ 679 PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival) 680 { 681 PetscErrorCode ierr; 682 683 PetscFunctionBegin; 684 ierr = PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 685 PetscFunctionReturn(0); 686 } 687 688 689 /*MC 690 MATSOLVERMKL_PARDISO - A matrix type providing direct solvers (LU) for 691 sequential matrices via the external package MKL_PARDISO. 692 693 Works with MATSEQAIJ matrices 694 695 Options Database Keys: 696 + -mat_mkl_pardiso_65 - Number of thrads to use 697 . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time 698 . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase 699 . -mat_mkl_pardiso_68 - Message level information 700 . -mat_mkl_pardiso_69 - Defines the matrix type. IMPORTANT: When you set this flag, iparm parameters are going to be set to the default ones for the matrix type 701 . -mat_mkl_pardiso_1 - Use default values 702 . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix 703 . -mat_mkl_pardiso_4 - Preconditioned CGS/CG 704 . -mat_mkl_pardiso_5 - User permutation 705 . -mat_mkl_pardiso_6 - Write solution on x 706 . -mat_mkl_pardiso_8 - Iterative refinement step 707 . -mat_mkl_pardiso_10 - Pivoting perturbation 708 . -mat_mkl_pardiso_11 - Scaling vectors 709 . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A 710 . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching 711 . -mat_mkl_pardiso_18 - Numbers of non-zero elements 712 . -mat_mkl_pardiso_19 - Report number of floating point operations 713 . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices 714 . -mat_mkl_pardiso_24 - Parallel factorization control 715 . -mat_mkl_pardiso_25 - Parallel forward/backward solve control 716 . -mat_mkl_pardiso_27 - Matrix checker 717 . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors 718 . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode 719 - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode 720 721 Level: beginner 722 723 For more information please check mkl_pardiso manual 724 725 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 726 727 M*/ 728 #undef __FUNCT__ 729 #define __FUNCT__ "MatFactorGetSolverPackage_mkl_pardiso" 730 static PetscErrorCode MatFactorGetSolverPackage_mkl_pardiso(Mat A, const MatSolverPackage *type) 731 { 732 PetscFunctionBegin; 733 *type = MATSOLVERMKL_PARDISO; 734 PetscFunctionReturn(0); 735 } 736 737 /* MatGetFactor for Seq sbAIJ matrices */ 738 #undef __FUNCT__ 739 #define __FUNCT__ "MatGetFactor_sbaij_mkl_pardiso" 740 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F) 741 { 742 Mat B; 743 PetscErrorCode ierr; 744 Mat_MKL_PARDISO *mat_mkl_pardiso; 745 PetscBool isSeqSBAIJ; 746 PetscInt bs; 747 748 PetscFunctionBegin; 749 /* Create the factorization matrix */ 750 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 751 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 752 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 753 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 754 if (isSeqSBAIJ) { 755 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 756 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Is not allowed other types of matrices apart from MATSEQSBAIJ."); 757 758 ierr = MatGetBlockSize(A,&bs); CHKERRQ(ierr); 759 760 if(bs != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQSBAIJ with block size other than 1 is not supported by Pardiso"); 761 762 if(ftype != MAT_FACTOR_CHOLESKY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQAIJ should be used only with MAT_FACTOR_CHOLESKY."); 763 764 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO; 765 B->factortype = MAT_FACTOR_CHOLESKY; 766 B->ops->destroy = MatDestroy_MKL_PARDISO; 767 B->ops->view = MatView_MKL_PARDISO; 768 B->factortype = ftype; 769 B->ops->getinfo = MatGetInfo_MKL_PARDISO; 770 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 771 772 ierr = PetscNewLog(B,&mat_mkl_pardiso);CHKERRQ(ierr); 773 B->spptr = mat_mkl_pardiso; 774 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);CHKERRQ(ierr); 775 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);CHKERRQ(ierr); 776 ierr = MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);CHKERRQ(ierr); 777 *F = B; 778 PetscFunctionReturn(0); 779 } 780 781 /* MatGetFactor for Seq AIJ matrices */ 782 #undef __FUNCT__ 783 #define __FUNCT__ "MatGetFactor_aij_mkl_pardiso" 784 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F) 785 { 786 Mat B; 787 PetscErrorCode ierr; 788 Mat_MKL_PARDISO *mat_mkl_pardiso; 789 PetscBool isSeqAIJ; 790 791 PetscFunctionBegin; 792 /* Create the factorization matrix */ 793 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 794 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 795 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 796 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 797 if (isSeqAIJ) { 798 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 799 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Is not allowed other types of matrices apart from MATSEQAIJ."); 800 801 if(ftype != MAT_FACTOR_LU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrice MATSEQAIJ should be used only with MAT_FACTOR_LU."); 802 803 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO; 804 B->factortype = MAT_FACTOR_LU; 805 B->ops->destroy = MatDestroy_MKL_PARDISO; 806 B->ops->view = MatView_MKL_PARDISO; 807 B->factortype = ftype; 808 B->ops->getinfo = MatGetInfo_MKL_PARDISO; 809 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 810 811 ierr = PetscNewLog(B,&mat_mkl_pardiso);CHKERRQ(ierr); 812 B->spptr = mat_mkl_pardiso; 813 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);CHKERRQ(ierr); 814 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);CHKERRQ(ierr); 815 ierr = MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);CHKERRQ(ierr); 816 817 *F = B; 818 PetscFunctionReturn(0); 819 } 820 821 #undef __FUNCT__ 822 #define __FUNCT__ "MatSolverPackageRegister_MKL_Pardiso" 823 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MKL_Pardiso(void) 824 { 825 PetscErrorCode ierr; 826 827 PetscFunctionBegin; 828 ierr = MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_aij_mkl_pardiso );CHKERRQ(ierr); 829 ierr = MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_sbaij_mkl_pardiso);CHKERRQ(ierr); 830 PetscFunctionReturn(0); 831 } 832 833