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