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