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