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 int iparm_copy[IPARM_SIZE], mtype_copy, i; 87 mtype_copy = *mtype; 88 pardisoinit(pt, &mtype_copy, iparm_copy); 89 for(i = 0; i < IPARM_SIZE; i++) 90 iparm[i] = iparm_copy[i]; 91 } 92 93 94 /* 95 * Copy the elements of matrix A. 96 * Input: 97 * - Mat A: MATSEQAIJ matrix 98 * - int shift: matrix index. 99 * - 0 for c representation 100 * - 1 for fortran representation 101 * - MatReuse reuse: 102 * - MAT_INITIAL_MATRIX: Create a new aij representation 103 * - MAT_REUSE_MATRIX: Reuse all aij representation and just change values 104 * Output: 105 * - int *nnz: Number of nonzero-elements. 106 * - int **r pointer to i index 107 * - int **c pointer to j elements 108 * - MATRIXTYPE **v: Non-zero elements 109 */ 110 #undef __FUNCT__ 111 #define __FUNCT__ "MatCopy_MKL_PARDISO" 112 PetscErrorCode MatCopy_MKL_PARDISO(Mat A, MatReuse reuse, INT_TYPE *nnz, INT_TYPE **r, INT_TYPE **c, void **v){ 113 114 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 115 116 PetscFunctionBegin; 117 *v=aa->a; 118 if (reuse == MAT_INITIAL_MATRIX) { 119 *r = (INT_TYPE*)aa->i; 120 *c = (INT_TYPE*)aa->j; 121 *nnz = aa->nz; 122 } 123 PetscFunctionReturn(0); 124 } 125 126 127 /* 128 * Free memory for Mat_MKL_PARDISO structure and pointers to objects. 129 */ 130 #undef __FUNCT__ 131 #define __FUNCT__ "MatDestroy_MKL_PARDISO" 132 PetscErrorCode MatDestroy_MKL_PARDISO(Mat A){ 133 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr; 134 PetscErrorCode ierr; 135 136 PetscFunctionBegin; 137 /* Terminate instance, deallocate memories */ 138 if (mat_mkl_pardiso->CleanUp) { 139 mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY; 140 141 142 MKL_PARDISO (mat_mkl_pardiso->pt, 143 &mat_mkl_pardiso->maxfct, 144 &mat_mkl_pardiso->mnum, 145 &mat_mkl_pardiso->mtype, 146 &mat_mkl_pardiso->phase, 147 &mat_mkl_pardiso->n, 148 NULL, 149 NULL, 150 NULL, 151 mat_mkl_pardiso->perm, 152 &mat_mkl_pardiso->nrhs, 153 mat_mkl_pardiso->iparm, 154 &mat_mkl_pardiso->msglvl, 155 NULL, 156 NULL, 157 &mat_mkl_pardiso->err); 158 } 159 ierr = PetscFree(mat_mkl_pardiso->perm);CHKERRQ(ierr); 160 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 161 162 /* clear composed functions */ 163 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);CHKERRQ(ierr); 164 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);CHKERRQ(ierr); 165 PetscFunctionReturn(0); 166 } 167 168 169 /* 170 * Computes Ax = b 171 */ 172 #undef __FUNCT__ 173 #define __FUNCT__ "MatSolve_MKL_PARDISO" 174 PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x){ 175 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr; 176 PetscErrorCode ierr; 177 PetscScalar *barray, *xarray; 178 179 PetscFunctionBegin; 180 181 182 mat_mkl_pardiso->nrhs = 1; 183 ierr = VecGetArray(x,&xarray);CHKERRQ(ierr); 184 ierr = VecGetArray(b,&barray);CHKERRQ(ierr); 185 186 /* solve phase */ 187 /*-------------*/ 188 mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 189 MKL_PARDISO (mat_mkl_pardiso->pt, 190 &mat_mkl_pardiso->maxfct, 191 &mat_mkl_pardiso->mnum, 192 &mat_mkl_pardiso->mtype, 193 &mat_mkl_pardiso->phase, 194 &mat_mkl_pardiso->n, 195 mat_mkl_pardiso->a, 196 mat_mkl_pardiso->ia, 197 mat_mkl_pardiso->ja, 198 mat_mkl_pardiso->perm, 199 &mat_mkl_pardiso->nrhs, 200 mat_mkl_pardiso->iparm, 201 &mat_mkl_pardiso->msglvl, 202 (void*)barray, 203 (void*)xarray, 204 &mat_mkl_pardiso->err); 205 206 207 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); 208 209 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 210 PetscFunctionReturn(0); 211 } 212 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "MatSolveTranspose_MKL_PARDISO" 216 PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x){ 217 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr; 218 PetscErrorCode ierr; 219 220 PetscFunctionBegin; 221 #if defined(PETSC_USE_COMPLEX) 222 mat_mkl_pardiso->iparm[12 - 1] = 1; 223 #else 224 mat_mkl_pardiso->iparm[12 - 1] = 2; 225 #endif 226 ierr = MatSolve_MKL_PARDISO(A,b,x);CHKERRQ(ierr); 227 mat_mkl_pardiso->iparm[12 - 1] = 0; 228 PetscFunctionReturn(0); 229 } 230 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatMatSolve_MKL_PARDISO" 234 PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X){ 235 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->spptr; 236 PetscErrorCode ierr; 237 PetscScalar *barray, *xarray; 238 PetscBool flg; 239 240 PetscFunctionBegin; 241 242 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 243 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix"); 244 ierr = PetscObjectTypeCompare((PetscObject)X,MATSEQDENSE,&flg);CHKERRQ(ierr); 245 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix"); 246 247 ierr = MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);CHKERRQ(ierr); 248 249 if(mat_mkl_pardiso->nrhs > 0){ 250 ierr = MatDenseGetArray(B,&barray); 251 ierr = MatDenseGetArray(X,&xarray); 252 253 /* solve phase */ 254 /*-------------*/ 255 mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 256 MKL_PARDISO (mat_mkl_pardiso->pt, 257 &mat_mkl_pardiso->maxfct, 258 &mat_mkl_pardiso->mnum, 259 &mat_mkl_pardiso->mtype, 260 &mat_mkl_pardiso->phase, 261 &mat_mkl_pardiso->n, 262 mat_mkl_pardiso->a, 263 mat_mkl_pardiso->ia, 264 mat_mkl_pardiso->ja, 265 mat_mkl_pardiso->perm, 266 &mat_mkl_pardiso->nrhs, 267 mat_mkl_pardiso->iparm, 268 &mat_mkl_pardiso->msglvl, 269 (void*)barray, 270 (void*)xarray, 271 &mat_mkl_pardiso->err); 272 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); 273 } 274 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 275 PetscFunctionReturn(0); 276 277 } 278 279 /* 280 * LU Decomposition 281 */ 282 #undef __FUNCT__ 283 #define __FUNCT__ "MatFactorNumeric_MKL_PARDISO" 284 PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info){ 285 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->spptr; 286 PetscErrorCode ierr; 287 288 /* numerical factorization phase */ 289 /*-------------------------------*/ 290 291 PetscFunctionBegin; 292 293 mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN; 294 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); 295 296 /* numerical factorization phase */ 297 /*-------------------------------*/ 298 mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION; 299 MKL_PARDISO (mat_mkl_pardiso->pt, 300 &mat_mkl_pardiso->maxfct, 301 &mat_mkl_pardiso->mnum, 302 &mat_mkl_pardiso->mtype, 303 &mat_mkl_pardiso->phase, 304 &mat_mkl_pardiso->n, 305 mat_mkl_pardiso->a, 306 mat_mkl_pardiso->ia, 307 mat_mkl_pardiso->ja, 308 mat_mkl_pardiso->perm, 309 &mat_mkl_pardiso->nrhs, 310 mat_mkl_pardiso->iparm, 311 &mat_mkl_pardiso->msglvl, 312 NULL, 313 NULL, 314 &mat_mkl_pardiso->err); 315 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); 316 317 mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN; 318 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 319 PetscFunctionReturn(0); 320 } 321 322 /* Sets mkl_pardiso options from the options database */ 323 #undef __FUNCT__ 324 #define __FUNCT__ "PetscSetMKL_PARDISOFromOptions" 325 PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A){ 326 Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr; 327 PetscErrorCode ierr; 328 PetscInt icntl; 329 PetscBool flg; 330 int pt[IPARM_SIZE], threads = 1; 331 332 PetscFunctionBegin; 333 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");CHKERRQ(ierr); 334 ierr = PetscOptionsInt("-mat_mkl_pardiso_65","Number of thrads to use","None",threads,&threads,&flg);CHKERRQ(ierr); 335 if (flg) mkl_set_num_threads(threads); 336 337 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); 338 if (flg) mat_mkl_pardiso->maxfct = icntl; 339 340 ierr = PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);CHKERRQ(ierr); 341 if (flg) mat_mkl_pardiso->mnum = icntl; 342 343 ierr = PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);CHKERRQ(ierr); 344 if (flg) mat_mkl_pardiso->msglvl = icntl; 345 346 ierr = PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);CHKERRQ(ierr); 347 if(flg){ 348 mat_mkl_pardiso->mtype = icntl; 349 MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm); 350 #if defined(PETSC_USE_REAL_SINGLE) 351 mat_mkl_pardiso->iparm[27] = 1; 352 #else 353 mat_mkl_pardiso->iparm[27] = 0; 354 #endif 355 mat_mkl_pardiso->iparm[34] = 1; 356 } 357 ierr = PetscOptionsInt("-mat_mkl_pardiso_1","Use default values","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);CHKERRQ(ierr); 358 359 if(flg && icntl != 0){ 360 ierr = PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);CHKERRQ(ierr); 361 if (flg) mat_mkl_pardiso->iparm[1] = icntl; 362 363 ierr = PetscOptionsInt("-mat_mkl_pardiso_4","Preconditioned CGS/CG","None",mat_mkl_pardiso->iparm[3],&icntl,&flg);CHKERRQ(ierr); 364 if (flg) mat_mkl_pardiso->iparm[3] = icntl; 365 366 ierr = PetscOptionsInt("-mat_mkl_pardiso_5","User permutation","None",mat_mkl_pardiso->iparm[4],&icntl,&flg);CHKERRQ(ierr); 367 if (flg) mat_mkl_pardiso->iparm[4] = icntl; 368 369 ierr = PetscOptionsInt("-mat_mkl_pardiso_6","Write solution on x","None",mat_mkl_pardiso->iparm[5],&icntl,&flg);CHKERRQ(ierr); 370 if (flg) mat_mkl_pardiso->iparm[5] = icntl; 371 372 ierr = PetscOptionsInt("-mat_mkl_pardiso_8","Iterative refinement step","None",mat_mkl_pardiso->iparm[7],&icntl,&flg);CHKERRQ(ierr); 373 if (flg) mat_mkl_pardiso->iparm[7] = icntl; 374 375 ierr = PetscOptionsInt("-mat_mkl_pardiso_10","Pivoting perturbation","None",mat_mkl_pardiso->iparm[9],&icntl,&flg);CHKERRQ(ierr); 376 if (flg) mat_mkl_pardiso->iparm[9] = icntl; 377 378 ierr = PetscOptionsInt("-mat_mkl_pardiso_11","Scaling vectors","None",mat_mkl_pardiso->iparm[10],&icntl,&flg);CHKERRQ(ierr); 379 if (flg) mat_mkl_pardiso->iparm[10] = icntl; 380 381 ierr = PetscOptionsInt("-mat_mkl_pardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_pardiso->iparm[11],&icntl,&flg);CHKERRQ(ierr); 382 if (flg) mat_mkl_pardiso->iparm[11] = icntl; 383 384 ierr = PetscOptionsInt("-mat_mkl_pardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_pardiso->iparm[12],&icntl,&flg);CHKERRQ(ierr); 385 if (flg) mat_mkl_pardiso->iparm[12] = icntl; 386 387 ierr = PetscOptionsInt("-mat_mkl_pardiso_18","Numbers of non-zero elements","None",mat_mkl_pardiso->iparm[17],&icntl,&flg);CHKERRQ(ierr); 388 if (flg) mat_mkl_pardiso->iparm[17] = icntl; 389 390 ierr = PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);CHKERRQ(ierr); 391 if (flg) mat_mkl_pardiso->iparm[18] = icntl; 392 393 ierr = PetscOptionsInt("-mat_mkl_pardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_pardiso->iparm[20],&icntl,&flg);CHKERRQ(ierr); 394 if (flg) mat_mkl_pardiso->iparm[20] = icntl; 395 396 ierr = PetscOptionsInt("-mat_mkl_pardiso_24","Parallel factorization control","None",mat_mkl_pardiso->iparm[23],&icntl,&flg);CHKERRQ(ierr); 397 if (flg) mat_mkl_pardiso->iparm[23] = icntl; 398 399 ierr = PetscOptionsInt("-mat_mkl_pardiso_25","Parallel forward/backward solve control","None",mat_mkl_pardiso->iparm[24],&icntl,&flg);CHKERRQ(ierr); 400 if (flg) mat_mkl_pardiso->iparm[24] = icntl; 401 402 ierr = PetscOptionsInt("-mat_mkl_pardiso_27","Matrix checker","None",mat_mkl_pardiso->iparm[26],&icntl,&flg);CHKERRQ(ierr); 403 if (flg) mat_mkl_pardiso->iparm[26] = icntl; 404 405 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); 406 if (flg) mat_mkl_pardiso->iparm[30] = icntl; 407 408 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); 409 if (flg) mat_mkl_pardiso->iparm[33] = icntl; 410 411 ierr = PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);CHKERRQ(ierr); 412 if (flg) mat_mkl_pardiso->iparm[59] = icntl; 413 } 414 PetscOptionsEnd(); 415 PetscFunctionReturn(0); 416 } 417 418 419 #undef __FUNCT__ 420 #define __FUNCT__ "PetscInitializeMKL_PARDISO" 421 PetscErrorCode PetscInitializeMKL_PARDISO(Mat A, Mat_MKL_PARDISO *mat_mkl_pardiso){ 422 PetscErrorCode ierr; 423 PetscInt i; 424 425 PetscFunctionBegin; 426 mat_mkl_pardiso->CleanUp = PETSC_FALSE; 427 mat_mkl_pardiso->maxfct = 1; 428 mat_mkl_pardiso->mnum = 1; 429 mat_mkl_pardiso->n = A->rmap->N; 430 mat_mkl_pardiso->msglvl = 0; 431 mat_mkl_pardiso->nrhs = 1; 432 mat_mkl_pardiso->err = 0; 433 mat_mkl_pardiso->phase = -1; 434 #if defined(PETSC_USE_COMPLEX) 435 mat_mkl_pardiso->mtype = 13; 436 #else 437 mat_mkl_pardiso->mtype = 11; 438 #endif 439 440 MKL_PARDISO_INIT(mat_mkl_pardiso->pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm); 441 442 #if defined(PETSC_USE_REAL_SINGLE) 443 mat_mkl_pardiso->iparm[27] = 1; 444 #else 445 mat_mkl_pardiso->iparm[27] = 0; 446 #endif 447 448 mat_mkl_pardiso->iparm[34] = 1; 449 450 ierr = PetscMalloc1(A->rmap->N, &mat_mkl_pardiso->perm);CHKERRQ(ierr); 451 for(i = 0; i < A->rmap->N; i++) 452 mat_mkl_pardiso->perm[i] = 0; 453 PetscFunctionReturn(0); 454 } 455 456 457 /* 458 * Symbolic decomposition. Mkl_Pardiso analysis phase. 459 */ 460 #undef __FUNCT__ 461 #define __FUNCT__ "MatLUFactorSymbolic_AIJMKL_PARDISO" 462 PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info){ 463 464 Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->spptr; 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 469 mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN; 470 471 /* Set MKL_PARDISO options from the options database */ 472 ierr = PetscSetMKL_PARDISOFromOptions(F,A);CHKERRQ(ierr); 473 474 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); 475 mat_mkl_pardiso->n = A->rmap->N; 476 477 /* analysis phase */ 478 /*----------------*/ 479 480 mat_mkl_pardiso->phase = JOB_ANALYSIS; 481 482 MKL_PARDISO (mat_mkl_pardiso->pt, 483 &mat_mkl_pardiso->maxfct, 484 &mat_mkl_pardiso->mnum, 485 &mat_mkl_pardiso->mtype, 486 &mat_mkl_pardiso->phase, 487 &mat_mkl_pardiso->n, 488 mat_mkl_pardiso->a, 489 mat_mkl_pardiso->ia, 490 mat_mkl_pardiso->ja, 491 mat_mkl_pardiso->perm, 492 &mat_mkl_pardiso->nrhs, 493 mat_mkl_pardiso->iparm, 494 &mat_mkl_pardiso->msglvl, 495 NULL, 496 NULL, 497 &mat_mkl_pardiso->err); 498 499 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); 500 501 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 502 F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO; 503 F->ops->solve = MatSolve_MKL_PARDISO; 504 F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO; 505 F->ops->matsolve = MatMatSolve_MKL_PARDISO; 506 PetscFunctionReturn(0); 507 } 508 509 510 #undef __FUNCT__ 511 #define __FUNCT__ "MatView_MKL_PARDISO" 512 PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer){ 513 PetscErrorCode ierr; 514 PetscBool iascii; 515 PetscViewerFormat format; 516 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->spptr; 517 PetscInt i; 518 519 PetscFunctionBegin; 520 /* check if matrix is mkl_pardiso type */ 521 if (A->ops->solve != MatSolve_MKL_PARDISO) PetscFunctionReturn(0); 522 523 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 524 if (iascii) { 525 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 526 if (format == PETSC_VIEWER_ASCII_INFO) { 527 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");CHKERRQ(ierr); 528 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase: %d \n",mat_mkl_pardiso->phase);CHKERRQ(ierr); 529 for(i = 1; i <= 64; i++){ 530 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]: %d \n",i, mat_mkl_pardiso->iparm[i - 1]);CHKERRQ(ierr); 531 } 532 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct);CHKERRQ(ierr); 533 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum: %d \n", mat_mkl_pardiso->mnum);CHKERRQ(ierr); 534 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype: %d \n", mat_mkl_pardiso->mtype);CHKERRQ(ierr); 535 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n: %d \n", mat_mkl_pardiso->n);CHKERRQ(ierr); 536 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs);CHKERRQ(ierr); 537 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl);CHKERRQ(ierr); 538 } 539 } 540 PetscFunctionReturn(0); 541 } 542 543 544 #undef __FUNCT__ 545 #define __FUNCT__ "MatGetInfo_MKL_PARDISO" 546 PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info){ 547 Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)A->spptr; 548 549 PetscFunctionBegin; 550 info->block_size = 1.0; 551 info->nz_allocated = mat_mkl_pardiso->nz + 0.0; 552 info->nz_unneeded = 0.0; 553 info->assemblies = 0.0; 554 info->mallocs = 0.0; 555 info->memory = 0.0; 556 info->fill_ratio_given = 0; 557 info->fill_ratio_needed = 0; 558 info->factor_mallocs = 0; 559 PetscFunctionReturn(0); 560 } 561 562 #undef __FUNCT__ 563 #define __FUNCT__ "MatMkl_PardisoSetCntl_MKL_PARDISO" 564 PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival){ 565 Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)F->spptr; 566 PetscFunctionBegin; 567 if(icntl <= 64){ 568 mat_mkl_pardiso->iparm[icntl - 1] = ival; 569 } else { 570 if(icntl == 65) 571 mkl_set_num_threads((int)ival); 572 else if(icntl == 66) 573 mat_mkl_pardiso->maxfct = ival; 574 else if(icntl == 67) 575 mat_mkl_pardiso->mnum = ival; 576 else if(icntl == 68) 577 mat_mkl_pardiso->msglvl = ival; 578 else if(icntl == 69){ 579 int pt[IPARM_SIZE]; 580 mat_mkl_pardiso->mtype = ival; 581 MKL_PARDISO_INIT(&pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm); 582 #if defined(PETSC_USE_REAL_SINGLE) 583 mat_mkl_pardiso->iparm[27] = 1; 584 #else 585 mat_mkl_pardiso->iparm[27] = 0; 586 #endif 587 mat_mkl_pardiso->iparm[34] = 1; 588 } 589 } 590 PetscFunctionReturn(0); 591 } 592 593 #undef __FUNCT__ 594 #define __FUNCT__ "MatMkl_PardisoSetCntl" 595 /*@ 596 MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters 597 598 Logically Collective on Mat 599 600 Input Parameters: 601 + F - the factored matrix obtained by calling MatGetFactor() 602 . icntl - index of Mkl_Pardiso parameter 603 - ival - value of Mkl_Pardiso parameter 604 605 Options Database: 606 . -mat_mkl_pardiso_<icntl> <ival> 607 608 Level: beginner 609 610 References: Mkl_Pardiso Users' Guide 611 612 .seealso: MatGetFactor() 613 @*/ 614 PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival) 615 { 616 PetscErrorCode ierr; 617 618 PetscFunctionBegin; 619 ierr = PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 620 PetscFunctionReturn(0); 621 } 622 623 624 /*MC 625 MATSOLVERMKL_PARDISO - A matrix type providing direct solvers (LU) for 626 sequential matrices via the external package MKL_PARDISO. 627 628 Works with MATSEQAIJ matrices 629 630 Options Database Keys: 631 + -mat_mkl_pardiso_65 - Number of thrads to use 632 . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time 633 . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase 634 . -mat_mkl_pardiso_68 - Message level information 635 . -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 636 . -mat_mkl_pardiso_1 - Use default values 637 . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix 638 . -mat_mkl_pardiso_4 - Preconditioned CGS/CG 639 . -mat_mkl_pardiso_5 - User permutation 640 . -mat_mkl_pardiso_6 - Write solution on x 641 . -mat_mkl_pardiso_8 - Iterative refinement step 642 . -mat_mkl_pardiso_10 - Pivoting perturbation 643 . -mat_mkl_pardiso_11 - Scaling vectors 644 . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A 645 . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching 646 . -mat_mkl_pardiso_18 - Numbers of non-zero elements 647 . -mat_mkl_pardiso_19 - Report number of floating point operations 648 . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices 649 . -mat_mkl_pardiso_24 - Parallel factorization control 650 . -mat_mkl_pardiso_25 - Parallel forward/backward solve control 651 . -mat_mkl_pardiso_27 - Matrix checker 652 . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors 653 . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode 654 - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode 655 656 Level: beginner 657 658 For more information please check mkl_pardiso manual 659 660 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 661 662 M*/ 663 664 665 #undef __FUNCT__ 666 #define __FUNCT__ "MatFactorGetSolverPackage_mkl_pardiso" 667 static PetscErrorCode MatFactorGetSolverPackage_mkl_pardiso(Mat A, const MatSolverPackage *type){ 668 PetscFunctionBegin; 669 *type = MATSOLVERMKL_PARDISO; 670 PetscFunctionReturn(0); 671 } 672 673 674 /* MatGetFactor for Seq AIJ matrices */ 675 #undef __FUNCT__ 676 #define __FUNCT__ "MatGetFactor_aij_mkl_pardiso" 677 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F) 678 { 679 Mat B; 680 PetscErrorCode ierr; 681 Mat_MKL_PARDISO *mat_mkl_pardiso; 682 683 PetscFunctionBegin; 684 /* Create the factorization matrix */ 685 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 686 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 687 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 688 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 689 690 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO; 691 B->ops->destroy = MatDestroy_MKL_PARDISO; 692 B->ops->view = MatView_MKL_PARDISO; 693 B->factortype = ftype; 694 B->ops->getinfo = MatGetInfo_MKL_PARDISO; 695 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 696 697 ierr = PetscNewLog(B,&mat_mkl_pardiso);CHKERRQ(ierr); 698 B->spptr = mat_mkl_pardiso; 699 700 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_mkl_pardiso);CHKERRQ(ierr); 701 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);CHKERRQ(ierr); 702 ierr = PetscInitializeMKL_PARDISO(A, mat_mkl_pardiso);CHKERRQ(ierr); 703 704 *F = B; 705 PetscFunctionReturn(0); 706 } 707 708 #undef __FUNCT__ 709 #define __FUNCT__ "MatSolverPackageRegister_MKL_Pardiso" 710 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_MKL_Pardiso(void) 711 { 712 PetscErrorCode ierr; 713 714 PetscFunctionBegin; 715 ierr = MatSolverPackageRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);CHKERRQ(ierr); 716 PetscFunctionReturn(0); 717 } 718 719