1 2 #include <petscsys.h> 3 #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I "petscmat.h" I*/ 4 #include <../src/mat/impls/sbaij/mpi/mpisbaij.h> 5 6 #if defined(PETSC_HAVE_MKL_INTEL_ILP64) 7 #define MKL_ILP64 8 #endif 9 #include <mkl.h> 10 #include <mkl_cluster_sparse_solver.h> 11 12 /* 13 * Possible mkl_cpardiso phases that controls the execution of the solver. 14 * For more information check mkl_cpardiso manual. 15 */ 16 #define JOB_ANALYSIS 11 17 #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION 12 18 #define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13 19 #define JOB_NUMERICAL_FACTORIZATION 22 20 #define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 23 21 #define JOB_SOLVE_ITERATIVE_REFINEMENT 33 22 #define JOB_SOLVE_FORWARD_SUBSTITUTION 331 23 #define JOB_SOLVE_DIAGONAL_SUBSTITUTION 332 24 #define JOB_SOLVE_BACKWARD_SUBSTITUTION 333 25 #define JOB_RELEASE_OF_LU_MEMORY 0 26 #define JOB_RELEASE_OF_ALL_MEMORY -1 27 28 #define IPARM_SIZE 64 29 #define INT_TYPE MKL_INT 30 31 static const char *Err_MSG_CPardiso(int errNo) { 32 switch (errNo) { 33 case -1: 34 return "input inconsistent"; break; 35 case -2: 36 return "not enough memory"; break; 37 case -3: 38 return "reordering problem"; break; 39 case -4: 40 return "zero pivot, numerical factorization or iterative refinement problem"; break; 41 case -5: 42 return "unclassified (internal) error"; break; 43 case -6: 44 return "preordering failed (matrix types 11, 13 only)"; break; 45 case -7: 46 return "diagonal matrix problem"; break; 47 case -8: 48 return "32-bit integer overflow problem"; break; 49 case -9: 50 return "not enough memory for OOC"; break; 51 case -10: 52 return "problems with opening OOC temporary files"; break; 53 case -11: 54 return "read/write problems with the OOC data file"; break; 55 default : 56 return "unknown error"; 57 } 58 } 59 60 /* 61 * Internal data structure. 62 * For more information check mkl_cpardiso manual. 63 */ 64 65 typedef struct { 66 67 /* Configuration vector */ 68 INT_TYPE iparm[IPARM_SIZE]; 69 70 /* 71 * Internal mkl_cpardiso memory location. 72 * After the first call to mkl_cpardiso do not modify pt, as that could cause a serious memory leak. 73 */ 74 void *pt[IPARM_SIZE]; 75 76 MPI_Fint comm_mkl_cpardiso; 77 78 /* Basic mkl_cpardiso info*/ 79 INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err; 80 81 /* Matrix structure */ 82 PetscScalar *a; 83 84 INT_TYPE *ia, *ja; 85 86 /* Number of non-zero elements */ 87 INT_TYPE nz; 88 89 /* Row permutaton vector*/ 90 INT_TYPE *perm; 91 92 /* Define is matrix preserve sparce structure. */ 93 MatStructure matstruc; 94 95 PetscErrorCode (*ConvertToTriples)(Mat, MatReuse, PetscInt*, PetscInt**, PetscInt**, PetscScalar**); 96 97 /* True if mkl_cpardiso function have been used. */ 98 PetscBool CleanUp; 99 } Mat_MKL_CPARDISO; 100 101 /* 102 * Copy the elements of matrix A. 103 * Input: 104 * - Mat A: MATSEQAIJ matrix 105 * - int shift: matrix index. 106 * - 0 for c representation 107 * - 1 for fortran representation 108 * - MatReuse reuse: 109 * - MAT_INITIAL_MATRIX: Create a new aij representation 110 * - MAT_REUSE_MATRIX: Reuse all aij representation and just change values 111 * Output: 112 * - int *nnz: Number of nonzero-elements. 113 * - int **r pointer to i index 114 * - int **c pointer to j elements 115 * - MATRIXTYPE **v: Non-zero elements 116 */ 117 PetscErrorCode MatCopy_seqaij_seqaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v) 118 { 119 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)A->data; 120 121 PetscFunctionBegin; 122 *v=aa->a; 123 if (reuse == MAT_INITIAL_MATRIX) { 124 *r = (INT_TYPE*)aa->i; 125 *c = (INT_TYPE*)aa->j; 126 *nnz = aa->nz; 127 } 128 PetscFunctionReturn(0); 129 } 130 131 PetscErrorCode MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v) 132 { 133 const PetscInt *ai, *aj, *bi, *bj,*garray,m=A->rmap->n,*ajj,*bjj; 134 PetscErrorCode ierr; 135 PetscInt rstart,nz,i,j,countA,countB; 136 PetscInt *row,*col; 137 const PetscScalar *av, *bv; 138 PetscScalar *val; 139 Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data; 140 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(mat->A)->data; 141 Mat_SeqAIJ *bb = (Mat_SeqAIJ*)(mat->B)->data; 142 PetscInt colA_start,jB,jcol; 143 144 PetscFunctionBegin; 145 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart; 146 av=aa->a; bv=bb->a; 147 148 garray = mat->garray; 149 150 if (reuse == MAT_INITIAL_MATRIX) { 151 nz = aa->nz + bb->nz; 152 *nnz = nz; 153 ierr = PetscMalloc3(m+1,&row,nz,&col,nz,&val);CHKERRQ(ierr); 154 *r = row; *c = col; *v = val; 155 } else { 156 row = *r; col = *c; val = *v; 157 } 158 159 nz = 0; 160 for (i=0; i<m; i++) { 161 row[i] = nz; 162 countA = ai[i+1] - ai[i]; 163 countB = bi[i+1] - bi[i]; 164 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 165 bjj = bj + bi[i]; 166 167 /* B part, smaller col index */ 168 colA_start = rstart + ajj[0]; /* the smallest global col index of A */ 169 jB = 0; 170 for (j=0; j<countB; j++) { 171 jcol = garray[bjj[j]]; 172 if (jcol > colA_start) break; 173 col[nz] = jcol; 174 val[nz++] = *bv++; 175 } 176 jB = j; 177 178 /* A part */ 179 for (j=0; j<countA; j++) { 180 col[nz] = rstart + ajj[j]; 181 val[nz++] = *av++; 182 } 183 184 /* B part, larger col index */ 185 for (j=jB; j<countB; j++) { 186 col[nz] = garray[bjj[j]]; 187 val[nz++] = *bv++; 188 } 189 } 190 row[m] = nz; 191 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v) 196 { 197 const PetscInt *ai, *aj, *bi, *bj,*garray,bs=A->rmap->bs,bs2=bs*bs,m=A->rmap->n/bs,*ajj,*bjj; 198 PetscErrorCode ierr; 199 PetscInt rstart,nz,i,j,countA,countB; 200 PetscInt *row,*col; 201 const PetscScalar *av, *bv; 202 PetscScalar *val; 203 Mat_MPIBAIJ *mat = (Mat_MPIBAIJ*)A->data; 204 Mat_SeqBAIJ *aa = (Mat_SeqBAIJ*)(mat->A)->data; 205 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 206 PetscInt colA_start,jB,jcol; 207 208 PetscFunctionBegin; 209 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart/bs; 210 av=aa->a; bv=bb->a; 211 212 garray = mat->garray; 213 214 if (reuse == MAT_INITIAL_MATRIX) { 215 nz = aa->nz + bb->nz; 216 *nnz = nz; 217 ierr = PetscMalloc3(m+1,&row,nz,&col,nz*bs2,&val);CHKERRQ(ierr); 218 *r = row; *c = col; *v = val; 219 } else { 220 row = *r; col = *c; val = *v; 221 } 222 223 nz = 0; 224 for (i=0; i<m; i++) { 225 row[i] = nz+1; 226 countA = ai[i+1] - ai[i]; 227 countB = bi[i+1] - bi[i]; 228 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 229 bjj = bj + bi[i]; 230 231 /* B part, smaller col index */ 232 colA_start = rstart + (countA > 0 ? ajj[0] : 0); /* the smallest global col index of A */ 233 jB = 0; 234 for (j=0; j<countB; j++) { 235 jcol = garray[bjj[j]]; 236 if (jcol > colA_start) break; 237 col[nz++] = jcol + 1; 238 } 239 jB = j; 240 ierr = PetscArraycpy(val,bv,jB*bs2);CHKERRQ(ierr); 241 val += jB*bs2; 242 bv += jB*bs2; 243 244 /* A part */ 245 for (j=0; j<countA; j++) col[nz++] = rstart + ajj[j] + 1; 246 ierr = PetscArraycpy(val,av,countA*bs2);CHKERRQ(ierr); 247 val += countA*bs2; 248 av += countA*bs2; 249 250 /* B part, larger col index */ 251 for (j=jB; j<countB; j++) col[nz++] = garray[bjj[j]] + 1; 252 ierr = PetscArraycpy(val,bv,(countB-jB)*bs2);CHKERRQ(ierr); 253 val += (countB-jB)*bs2; 254 bv += (countB-jB)*bs2; 255 } 256 row[m] = nz+1; 257 258 PetscFunctionReturn(0); 259 } 260 261 PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v) 262 { 263 const PetscInt *ai, *aj, *bi, *bj,*garray,bs=A->rmap->bs,bs2=bs*bs,m=A->rmap->n/bs,*ajj,*bjj; 264 PetscErrorCode ierr; 265 PetscInt rstart,nz,i,j,countA,countB; 266 PetscInt *row,*col; 267 const PetscScalar *av, *bv; 268 PetscScalar *val; 269 Mat_MPISBAIJ *mat = (Mat_MPISBAIJ*)A->data; 270 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)(mat->A)->data; 271 Mat_SeqBAIJ *bb = (Mat_SeqBAIJ*)(mat->B)->data; 272 273 PetscFunctionBegin; 274 ai=aa->i; aj=aa->j; bi=bb->i; bj=bb->j; rstart=A->rmap->rstart/bs; 275 av=aa->a; bv=bb->a; 276 277 garray = mat->garray; 278 279 if (reuse == MAT_INITIAL_MATRIX) { 280 nz = aa->nz + bb->nz; 281 *nnz = nz; 282 ierr = PetscMalloc3(m+1,&row,nz,&col,nz*bs2,&val);CHKERRQ(ierr); 283 *r = row; *c = col; *v = val; 284 } else { 285 row = *r; col = *c; val = *v; 286 } 287 288 nz = 0; 289 for (i=0; i<m; i++) { 290 row[i] = nz+1; 291 countA = ai[i+1] - ai[i]; 292 countB = bi[i+1] - bi[i]; 293 ajj = aj + ai[i]; /* ptr to the beginning of this row */ 294 bjj = bj + bi[i]; 295 296 /* A part */ 297 for (j=0; j<countA; j++) col[nz++] = rstart + ajj[j] + 1; 298 ierr = PetscArraycpy(val,av,countA*bs2);CHKERRQ(ierr); 299 val += countA*bs2; 300 av += countA*bs2; 301 302 /* B part, larger col index */ 303 for (j=0; j<countB; j++) col[nz++] = garray[bjj[j]] + 1; 304 ierr = PetscArraycpy(val,bv,countB*bs2);CHKERRQ(ierr); 305 val += countB*bs2; 306 bv += countB*bs2; 307 } 308 row[m] = nz+1; 309 310 PetscFunctionReturn(0); 311 } 312 313 /* 314 * Free memory for Mat_MKL_CPARDISO structure and pointers to objects. 315 */ 316 PetscErrorCode MatDestroy_MKL_CPARDISO(Mat A) 317 { 318 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data; 319 MPI_Comm comm; 320 PetscErrorCode ierr; 321 322 PetscFunctionBegin; 323 /* Terminate instance, deallocate memories */ 324 if (mat_mkl_cpardiso->CleanUp) { 325 mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY; 326 327 cluster_sparse_solver ( 328 mat_mkl_cpardiso->pt, 329 &mat_mkl_cpardiso->maxfct, 330 &mat_mkl_cpardiso->mnum, 331 &mat_mkl_cpardiso->mtype, 332 &mat_mkl_cpardiso->phase, 333 &mat_mkl_cpardiso->n, 334 NULL, 335 NULL, 336 NULL, 337 mat_mkl_cpardiso->perm, 338 &mat_mkl_cpardiso->nrhs, 339 mat_mkl_cpardiso->iparm, 340 &mat_mkl_cpardiso->msglvl, 341 NULL, 342 NULL, 343 &mat_mkl_cpardiso->comm_mkl_cpardiso, 344 (PetscInt*)&mat_mkl_cpardiso->err); 345 } 346 347 if (mat_mkl_cpardiso->ConvertToTriples != MatCopy_seqaij_seqaij_MKL_CPARDISO) { 348 ierr = PetscFree3(mat_mkl_cpardiso->ia,mat_mkl_cpardiso->ja,mat_mkl_cpardiso->a);CHKERRQ(ierr); 349 } 350 comm = MPI_Comm_f2c(mat_mkl_cpardiso->comm_mkl_cpardiso); 351 ierr = MPI_Comm_free(&comm);CHKERRMPI(ierr); 352 ierr = PetscFree(A->data);CHKERRQ(ierr); 353 354 /* clear composed functions */ 355 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);CHKERRQ(ierr); 356 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMkl_CPardisoSetCntl_C",NULL);CHKERRQ(ierr); 357 PetscFunctionReturn(0); 358 } 359 360 /* 361 * Computes Ax = b 362 */ 363 PetscErrorCode MatSolve_MKL_CPARDISO(Mat A,Vec b,Vec x) 364 { 365 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data; 366 PetscErrorCode ierr; 367 PetscScalar *xarray; 368 const PetscScalar *barray; 369 370 PetscFunctionBegin; 371 mat_mkl_cpardiso->nrhs = 1; 372 ierr = VecGetArray(x,&xarray);CHKERRQ(ierr); 373 ierr = VecGetArrayRead(b,&barray);CHKERRQ(ierr); 374 375 /* solve phase */ 376 /*-------------*/ 377 mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 378 cluster_sparse_solver ( 379 mat_mkl_cpardiso->pt, 380 &mat_mkl_cpardiso->maxfct, 381 &mat_mkl_cpardiso->mnum, 382 &mat_mkl_cpardiso->mtype, 383 &mat_mkl_cpardiso->phase, 384 &mat_mkl_cpardiso->n, 385 mat_mkl_cpardiso->a, 386 mat_mkl_cpardiso->ia, 387 mat_mkl_cpardiso->ja, 388 mat_mkl_cpardiso->perm, 389 &mat_mkl_cpardiso->nrhs, 390 mat_mkl_cpardiso->iparm, 391 &mat_mkl_cpardiso->msglvl, 392 (void*)barray, 393 (void*)xarray, 394 &mat_mkl_cpardiso->comm_mkl_cpardiso, 395 (PetscInt*)&mat_mkl_cpardiso->err); 396 397 if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err)); 398 399 ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr); 400 ierr = VecRestoreArrayRead(b,&barray);CHKERRQ(ierr); 401 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 402 PetscFunctionReturn(0); 403 } 404 405 PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A,Vec b,Vec x) 406 { 407 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 #if defined(PETSC_USE_COMPLEX) 412 mat_mkl_cpardiso->iparm[12 - 1] = 1; 413 #else 414 mat_mkl_cpardiso->iparm[12 - 1] = 2; 415 #endif 416 ierr = MatSolve_MKL_CPARDISO(A,b,x);CHKERRQ(ierr); 417 mat_mkl_cpardiso->iparm[12 - 1] = 0; 418 PetscFunctionReturn(0); 419 } 420 421 PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A,Mat B,Mat X) 422 { 423 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data; 424 PetscErrorCode ierr; 425 PetscScalar *xarray; 426 const PetscScalar *barray; 427 428 PetscFunctionBegin; 429 ierr = MatGetSize(B,NULL,(PetscInt*)&mat_mkl_cpardiso->nrhs);CHKERRQ(ierr); 430 431 if (mat_mkl_cpardiso->nrhs > 0) { 432 ierr = MatDenseGetArrayRead(B,&barray);CHKERRQ(ierr); 433 ierr = MatDenseGetArray(X,&xarray);CHKERRQ(ierr); 434 435 if (barray == xarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"B and X cannot share the same memory location"); 436 437 /* solve phase */ 438 /*-------------*/ 439 mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 440 cluster_sparse_solver ( 441 mat_mkl_cpardiso->pt, 442 &mat_mkl_cpardiso->maxfct, 443 &mat_mkl_cpardiso->mnum, 444 &mat_mkl_cpardiso->mtype, 445 &mat_mkl_cpardiso->phase, 446 &mat_mkl_cpardiso->n, 447 mat_mkl_cpardiso->a, 448 mat_mkl_cpardiso->ia, 449 mat_mkl_cpardiso->ja, 450 mat_mkl_cpardiso->perm, 451 &mat_mkl_cpardiso->nrhs, 452 mat_mkl_cpardiso->iparm, 453 &mat_mkl_cpardiso->msglvl, 454 (void*)barray, 455 (void*)xarray, 456 &mat_mkl_cpardiso->comm_mkl_cpardiso, 457 (PetscInt*)&mat_mkl_cpardiso->err); 458 if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err)); 459 ierr = MatDenseRestoreArrayRead(B,&barray);CHKERRQ(ierr); 460 ierr = MatDenseRestoreArray(X,&xarray);CHKERRQ(ierr); 461 462 } 463 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 464 PetscFunctionReturn(0); 465 } 466 467 /* 468 * LU Decomposition 469 */ 470 PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F,Mat A,const MatFactorInfo *info) 471 { 472 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(F)->data; 473 PetscErrorCode ierr; 474 475 PetscFunctionBegin; 476 mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN; 477 ierr = (*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_REUSE_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);CHKERRQ(ierr); 478 479 mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION; 480 cluster_sparse_solver ( 481 mat_mkl_cpardiso->pt, 482 &mat_mkl_cpardiso->maxfct, 483 &mat_mkl_cpardiso->mnum, 484 &mat_mkl_cpardiso->mtype, 485 &mat_mkl_cpardiso->phase, 486 &mat_mkl_cpardiso->n, 487 mat_mkl_cpardiso->a, 488 mat_mkl_cpardiso->ia, 489 mat_mkl_cpardiso->ja, 490 mat_mkl_cpardiso->perm, 491 &mat_mkl_cpardiso->nrhs, 492 mat_mkl_cpardiso->iparm, 493 &mat_mkl_cpardiso->msglvl, 494 NULL, 495 NULL, 496 &mat_mkl_cpardiso->comm_mkl_cpardiso, 497 &mat_mkl_cpardiso->err); 498 if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\". Please check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err)); 499 500 mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN; 501 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 502 PetscFunctionReturn(0); 503 } 504 505 /* Sets mkl_cpardiso options from the options database */ 506 PetscErrorCode PetscSetMKL_CPARDISOFromOptions(Mat F, Mat A) 507 { 508 Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data; 509 PetscErrorCode ierr; 510 PetscInt icntl,threads; 511 PetscBool flg; 512 513 PetscFunctionBegin; 514 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_CPARDISO Options","Mat");CHKERRQ(ierr); 515 ierr = PetscOptionsInt("-mat_mkl_cpardiso_65","Suggested number of threads to use within MKL_CPARDISO","None",threads,&threads,&flg);CHKERRQ(ierr); 516 if (flg) mkl_set_num_threads((int)threads); 517 518 ierr = PetscOptionsInt("-mat_mkl_cpardiso_66","Maximum number of factors with identical sparsity structure that must be kept in memory at the same time","None",mat_mkl_cpardiso->maxfct,&icntl,&flg);CHKERRQ(ierr); 519 if (flg) mat_mkl_cpardiso->maxfct = icntl; 520 521 ierr = PetscOptionsInt("-mat_mkl_cpardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_cpardiso->mnum,&icntl,&flg);CHKERRQ(ierr); 522 if (flg) mat_mkl_cpardiso->mnum = icntl; 523 524 ierr = PetscOptionsInt("-mat_mkl_cpardiso_68","Message level information","None",mat_mkl_cpardiso->msglvl,&icntl,&flg);CHKERRQ(ierr); 525 if (flg) mat_mkl_cpardiso->msglvl = icntl; 526 527 ierr = PetscOptionsInt("-mat_mkl_cpardiso_69","Defines the matrix type","None",mat_mkl_cpardiso->mtype,&icntl,&flg);CHKERRQ(ierr); 528 if (flg) mat_mkl_cpardiso->mtype = icntl; 529 ierr = PetscOptionsInt("-mat_mkl_cpardiso_1","Use default values","None",mat_mkl_cpardiso->iparm[0],&icntl,&flg);CHKERRQ(ierr); 530 531 if (flg && icntl != 0) { 532 ierr = PetscOptionsInt("-mat_mkl_cpardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_cpardiso->iparm[1],&icntl,&flg);CHKERRQ(ierr); 533 if (flg) mat_mkl_cpardiso->iparm[1] = icntl; 534 535 ierr = PetscOptionsInt("-mat_mkl_cpardiso_4","Preconditioned CGS/CG","None",mat_mkl_cpardiso->iparm[3],&icntl,&flg);CHKERRQ(ierr); 536 if (flg) mat_mkl_cpardiso->iparm[3] = icntl; 537 538 ierr = PetscOptionsInt("-mat_mkl_cpardiso_5","User permutation","None",mat_mkl_cpardiso->iparm[4],&icntl,&flg);CHKERRQ(ierr); 539 if (flg) mat_mkl_cpardiso->iparm[4] = icntl; 540 541 ierr = PetscOptionsInt("-mat_mkl_cpardiso_6","Write solution on x","None",mat_mkl_cpardiso->iparm[5],&icntl,&flg);CHKERRQ(ierr); 542 if (flg) mat_mkl_cpardiso->iparm[5] = icntl; 543 544 ierr = PetscOptionsInt("-mat_mkl_cpardiso_8","Iterative refinement step","None",mat_mkl_cpardiso->iparm[7],&icntl,&flg);CHKERRQ(ierr); 545 if (flg) mat_mkl_cpardiso->iparm[7] = icntl; 546 547 ierr = PetscOptionsInt("-mat_mkl_cpardiso_10","Pivoting perturbation","None",mat_mkl_cpardiso->iparm[9],&icntl,&flg);CHKERRQ(ierr); 548 if (flg) mat_mkl_cpardiso->iparm[9] = icntl; 549 550 ierr = PetscOptionsInt("-mat_mkl_cpardiso_11","Scaling vectors","None",mat_mkl_cpardiso->iparm[10],&icntl,&flg);CHKERRQ(ierr); 551 if (flg) mat_mkl_cpardiso->iparm[10] = icntl; 552 553 ierr = PetscOptionsInt("-mat_mkl_cpardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_cpardiso->iparm[11],&icntl,&flg);CHKERRQ(ierr); 554 if (flg) mat_mkl_cpardiso->iparm[11] = icntl; 555 556 ierr = PetscOptionsInt("-mat_mkl_cpardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_cpardiso->iparm[12],&icntl, 557 &flg);CHKERRQ(ierr); 558 if (flg) mat_mkl_cpardiso->iparm[12] = icntl; 559 560 ierr = PetscOptionsInt("-mat_mkl_cpardiso_18","Numbers of non-zero elements","None",mat_mkl_cpardiso->iparm[17],&icntl, 561 &flg);CHKERRQ(ierr); 562 if (flg) mat_mkl_cpardiso->iparm[17] = icntl; 563 564 ierr = PetscOptionsInt("-mat_mkl_cpardiso_19","Report number of floating point operations","None",mat_mkl_cpardiso->iparm[18],&icntl,&flg);CHKERRQ(ierr); 565 if (flg) mat_mkl_cpardiso->iparm[18] = icntl; 566 567 ierr = PetscOptionsInt("-mat_mkl_cpardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_cpardiso->iparm[20],&icntl,&flg);CHKERRQ(ierr); 568 if (flg) mat_mkl_cpardiso->iparm[20] = icntl; 569 570 ierr = PetscOptionsInt("-mat_mkl_cpardiso_24","Parallel factorization control","None",mat_mkl_cpardiso->iparm[23],&icntl,&flg);CHKERRQ(ierr); 571 if (flg) mat_mkl_cpardiso->iparm[23] = icntl; 572 573 ierr = PetscOptionsInt("-mat_mkl_cpardiso_25","Parallel forward/backward solve control","None",mat_mkl_cpardiso->iparm[24],&icntl,&flg);CHKERRQ(ierr); 574 if (flg) mat_mkl_cpardiso->iparm[24] = icntl; 575 576 ierr = PetscOptionsInt("-mat_mkl_cpardiso_27","Matrix checker","None",mat_mkl_cpardiso->iparm[26],&icntl,&flg);CHKERRQ(ierr); 577 if (flg) mat_mkl_cpardiso->iparm[26] = icntl; 578 579 ierr = PetscOptionsInt("-mat_mkl_cpardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_cpardiso->iparm[30],&icntl,&flg);CHKERRQ(ierr); 580 if (flg) mat_mkl_cpardiso->iparm[30] = icntl; 581 582 ierr = PetscOptionsInt("-mat_mkl_cpardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_cpardiso->iparm[33],&icntl,&flg);CHKERRQ(ierr); 583 if (flg) mat_mkl_cpardiso->iparm[33] = icntl; 584 585 ierr = PetscOptionsInt("-mat_mkl_cpardiso_60","Intel MKL_CPARDISO mode","None",mat_mkl_cpardiso->iparm[59],&icntl,&flg);CHKERRQ(ierr); 586 if (flg) mat_mkl_cpardiso->iparm[59] = icntl; 587 } 588 589 ierr = PetscOptionsEnd();CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 593 PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso) 594 { 595 PetscErrorCode ierr; 596 PetscInt bs; 597 PetscBool match; 598 PetscMPIInt size; 599 MPI_Comm comm; 600 601 PetscFunctionBegin; 602 603 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&comm);CHKERRMPI(ierr); 604 ierr = MPI_Comm_size(comm, &size);CHKERRMPI(ierr); 605 mat_mkl_cpardiso->comm_mkl_cpardiso = MPI_Comm_c2f(comm); 606 607 mat_mkl_cpardiso->CleanUp = PETSC_FALSE; 608 mat_mkl_cpardiso->maxfct = 1; 609 mat_mkl_cpardiso->mnum = 1; 610 mat_mkl_cpardiso->n = A->rmap->N; 611 if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36]; 612 mat_mkl_cpardiso->msglvl = 0; 613 mat_mkl_cpardiso->nrhs = 1; 614 mat_mkl_cpardiso->err = 0; 615 mat_mkl_cpardiso->phase = -1; 616 #if defined(PETSC_USE_COMPLEX) 617 mat_mkl_cpardiso->mtype = 13; 618 #else 619 mat_mkl_cpardiso->mtype = 11; 620 #endif 621 622 #if defined(PETSC_USE_REAL_SINGLE) 623 mat_mkl_cpardiso->iparm[27] = 1; 624 #else 625 mat_mkl_cpardiso->iparm[27] = 0; 626 #endif 627 628 mat_mkl_cpardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */ 629 mat_mkl_cpardiso->iparm[ 1] = 2; /* Use METIS for fill-in reordering */ 630 mat_mkl_cpardiso->iparm[ 5] = 0; /* Write solution into x */ 631 mat_mkl_cpardiso->iparm[ 7] = 2; /* Max number of iterative refinement steps */ 632 mat_mkl_cpardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */ 633 mat_mkl_cpardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ 634 mat_mkl_cpardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */ 635 mat_mkl_cpardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */ 636 mat_mkl_cpardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */ 637 mat_mkl_cpardiso->iparm[26] = 1; /* Check input data for correctness */ 638 639 mat_mkl_cpardiso->iparm[39] = 0; 640 if (size > 1) { 641 mat_mkl_cpardiso->iparm[39] = 2; 642 mat_mkl_cpardiso->iparm[40] = A->rmap->rstart; 643 mat_mkl_cpardiso->iparm[41] = A->rmap->rend-1; 644 } 645 ierr = PetscObjectTypeCompareAny((PetscObject)A,&match,MATMPIBAIJ,MATMPISBAIJ,"");CHKERRQ(ierr); 646 if (match) { 647 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 648 mat_mkl_cpardiso->iparm[36] = bs; 649 mat_mkl_cpardiso->iparm[40] /= bs; 650 mat_mkl_cpardiso->iparm[41] /= bs; 651 mat_mkl_cpardiso->iparm[40]++; 652 mat_mkl_cpardiso->iparm[41]++; 653 mat_mkl_cpardiso->iparm[34] = 0; /* Fortran style */ 654 } else { 655 mat_mkl_cpardiso->iparm[34] = 1; /* C style */ 656 } 657 658 mat_mkl_cpardiso->perm = 0; 659 PetscFunctionReturn(0); 660 } 661 662 /* 663 * Symbolic decomposition. Mkl_Pardiso analysis phase. 664 */ 665 PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 666 { 667 Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data; 668 PetscErrorCode ierr; 669 670 PetscFunctionBegin; 671 mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN; 672 673 /* Set MKL_CPARDISO options from the options database */ 674 ierr = PetscSetMKL_CPARDISOFromOptions(F,A);CHKERRQ(ierr); 675 ierr = (*mat_mkl_cpardiso->ConvertToTriples)(A,MAT_INITIAL_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);CHKERRQ(ierr); 676 677 mat_mkl_cpardiso->n = A->rmap->N; 678 if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36]; 679 680 /* analysis phase */ 681 /*----------------*/ 682 mat_mkl_cpardiso->phase = JOB_ANALYSIS; 683 684 cluster_sparse_solver ( 685 mat_mkl_cpardiso->pt, 686 &mat_mkl_cpardiso->maxfct, 687 &mat_mkl_cpardiso->mnum, 688 &mat_mkl_cpardiso->mtype, 689 &mat_mkl_cpardiso->phase, 690 &mat_mkl_cpardiso->n, 691 mat_mkl_cpardiso->a, 692 mat_mkl_cpardiso->ia, 693 mat_mkl_cpardiso->ja, 694 mat_mkl_cpardiso->perm, 695 &mat_mkl_cpardiso->nrhs, 696 mat_mkl_cpardiso->iparm, 697 &mat_mkl_cpardiso->msglvl, 698 NULL, 699 NULL, 700 &mat_mkl_cpardiso->comm_mkl_cpardiso, 701 (PetscInt*)&mat_mkl_cpardiso->err); 702 703 if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\".Check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err)); 704 705 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 706 F->ops->lufactornumeric = MatFactorNumeric_MKL_CPARDISO; 707 F->ops->solve = MatSolve_MKL_CPARDISO; 708 F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO; 709 F->ops->matsolve = MatMatSolve_MKL_CPARDISO; 710 PetscFunctionReturn(0); 711 } 712 713 PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS perm,const MatFactorInfo *info) 714 { 715 Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data; 716 PetscErrorCode ierr; 717 718 PetscFunctionBegin; 719 mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN; 720 721 /* Set MKL_CPARDISO options from the options database */ 722 ierr = PetscSetMKL_CPARDISOFromOptions(F,A);CHKERRQ(ierr); 723 ierr = (*mat_mkl_cpardiso->ConvertToTriples)(A,MAT_INITIAL_MATRIX,&mat_mkl_cpardiso->nz,&mat_mkl_cpardiso->ia,&mat_mkl_cpardiso->ja,&mat_mkl_cpardiso->a);CHKERRQ(ierr); 724 725 mat_mkl_cpardiso->n = A->rmap->N; 726 if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36]; 727 #if defined(PETSC_USE_COMPLEX) 728 SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead",((PetscObject)A)->type_name); 729 #endif 730 if (A->spd_set && A->spd) mat_mkl_cpardiso->mtype = 2; 731 else mat_mkl_cpardiso->mtype = -2; 732 733 /* analysis phase */ 734 /*----------------*/ 735 mat_mkl_cpardiso->phase = JOB_ANALYSIS; 736 737 cluster_sparse_solver ( 738 mat_mkl_cpardiso->pt, 739 &mat_mkl_cpardiso->maxfct, 740 &mat_mkl_cpardiso->mnum, 741 &mat_mkl_cpardiso->mtype, 742 &mat_mkl_cpardiso->phase, 743 &mat_mkl_cpardiso->n, 744 mat_mkl_cpardiso->a, 745 mat_mkl_cpardiso->ia, 746 mat_mkl_cpardiso->ja, 747 mat_mkl_cpardiso->perm, 748 &mat_mkl_cpardiso->nrhs, 749 mat_mkl_cpardiso->iparm, 750 &mat_mkl_cpardiso->msglvl, 751 NULL, 752 NULL, 753 &mat_mkl_cpardiso->comm_mkl_cpardiso, 754 (PetscInt*)&mat_mkl_cpardiso->err); 755 756 if (mat_mkl_cpardiso->err < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_CPARDISO: err=%d, msg = \"%s\".Check manual\n",mat_mkl_cpardiso->err,Err_MSG_CPardiso(mat_mkl_cpardiso->err)); 757 758 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 759 F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_CPARDISO; 760 F->ops->solve = MatSolve_MKL_CPARDISO; 761 F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO; 762 F->ops->matsolve = MatMatSolve_MKL_CPARDISO; 763 PetscFunctionReturn(0); 764 } 765 766 PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer) 767 { 768 PetscErrorCode ierr; 769 PetscBool iascii; 770 PetscViewerFormat format; 771 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data; 772 PetscInt i; 773 774 PetscFunctionBegin; 775 /* check if matrix is mkl_cpardiso type */ 776 if (A->ops->solve != MatSolve_MKL_CPARDISO) PetscFunctionReturn(0); 777 778 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 779 if (iascii) { 780 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 781 if (format == PETSC_VIEWER_ASCII_INFO) { 782 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO run parameters:\n");CHKERRQ(ierr); 783 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO phase: %d \n",mat_mkl_cpardiso->phase);CHKERRQ(ierr); 784 for (i = 1; i <= 64; i++) { 785 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO iparm[%d]: %d \n",i, mat_mkl_cpardiso->iparm[i - 1]);CHKERRQ(ierr); 786 } 787 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO maxfct: %d \n", mat_mkl_cpardiso->maxfct);CHKERRQ(ierr); 788 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mnum: %d \n", mat_mkl_cpardiso->mnum);CHKERRQ(ierr); 789 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mtype: %d \n", mat_mkl_cpardiso->mtype);CHKERRQ(ierr); 790 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO n: %d \n", mat_mkl_cpardiso->n);CHKERRQ(ierr); 791 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO nrhs: %d \n", mat_mkl_cpardiso->nrhs);CHKERRQ(ierr); 792 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO msglvl: %d \n", mat_mkl_cpardiso->msglvl);CHKERRQ(ierr); 793 } 794 } 795 PetscFunctionReturn(0); 796 } 797 798 PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info) 799 { 800 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data; 801 802 PetscFunctionBegin; 803 info->block_size = 1.0; 804 info->nz_allocated = mat_mkl_cpardiso->nz + 0.0; 805 info->nz_unneeded = 0.0; 806 info->assemblies = 0.0; 807 info->mallocs = 0.0; 808 info->memory = 0.0; 809 info->fill_ratio_given = 0; 810 info->fill_ratio_needed = 0; 811 info->factor_mallocs = 0; 812 PetscFunctionReturn(0); 813 } 814 815 PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F,PetscInt icntl,PetscInt ival) 816 { 817 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)F->data; 818 819 PetscFunctionBegin; 820 if (icntl <= 64) { 821 mat_mkl_cpardiso->iparm[icntl - 1] = ival; 822 } else { 823 if (icntl == 65) mkl_set_num_threads((int)ival); 824 else if (icntl == 66) mat_mkl_cpardiso->maxfct = ival; 825 else if (icntl == 67) mat_mkl_cpardiso->mnum = ival; 826 else if (icntl == 68) mat_mkl_cpardiso->msglvl = ival; 827 else if (icntl == 69) mat_mkl_cpardiso->mtype = ival; 828 } 829 PetscFunctionReturn(0); 830 } 831 832 /*@ 833 MatMkl_CPardisoSetCntl - Set Mkl_Pardiso parameters 834 835 Logically Collective on Mat 836 837 Input Parameters: 838 + F - the factored matrix obtained by calling MatGetFactor() 839 . icntl - index of Mkl_Pardiso parameter 840 - ival - value of Mkl_Pardiso parameter 841 842 Options Database: 843 . -mat_mkl_cpardiso_<icntl> <ival> 844 845 Level: Intermediate 846 847 Notes: 848 This routine cannot be used if you are solving the linear system with TS, SNES, or KSP, only if you directly call MatGetFactor() so use the options 849 database approach when working with TS, SNES, or KSP. 850 851 References: 852 . Mkl_Pardiso Users' Guide 853 854 .seealso: MatGetFactor() 855 @*/ 856 PetscErrorCode MatMkl_CPardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival) 857 { 858 PetscErrorCode ierr; 859 860 PetscFunctionBegin; 861 ierr = PetscTryMethod(F,"MatMkl_CPardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 862 PetscFunctionReturn(0); 863 } 864 865 /*MC 866 MATSOLVERMKL_CPARDISO - A matrix type providing direct solvers (LU) for parallel matrices via the external package MKL_CPARDISO. 867 868 Works with MATMPIAIJ matrices 869 870 Use -pc_type lu -pc_factor_mat_solver_type mkl_cpardiso to use this direct solver 871 872 Options Database Keys: 873 + -mat_mkl_cpardiso_65 - Suggested number of threads to use within MKL_CPARDISO 874 . -mat_mkl_cpardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time 875 . -mat_mkl_cpardiso_67 - Indicates the actual matrix for the solution phase 876 . -mat_mkl_cpardiso_68 - Message level information, use 1 to get detailed information on the solver options 877 . -mat_mkl_cpardiso_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 878 . -mat_mkl_cpardiso_1 - Use default values 879 . -mat_mkl_cpardiso_2 - Fill-in reducing ordering for the input matrix 880 . -mat_mkl_cpardiso_4 - Preconditioned CGS/CG 881 . -mat_mkl_cpardiso_5 - User permutation 882 . -mat_mkl_cpardiso_6 - Write solution on x 883 . -mat_mkl_cpardiso_8 - Iterative refinement step 884 . -mat_mkl_cpardiso_10 - Pivoting perturbation 885 . -mat_mkl_cpardiso_11 - Scaling vectors 886 . -mat_mkl_cpardiso_12 - Solve with transposed or conjugate transposed matrix A 887 . -mat_mkl_cpardiso_13 - Improved accuracy using (non-) symmetric weighted matching 888 . -mat_mkl_cpardiso_18 - Numbers of non-zero elements 889 . -mat_mkl_cpardiso_19 - Report number of floating point operations 890 . -mat_mkl_cpardiso_21 - Pivoting for symmetric indefinite matrices 891 . -mat_mkl_cpardiso_24 - Parallel factorization control 892 . -mat_mkl_cpardiso_25 - Parallel forward/backward solve control 893 . -mat_mkl_cpardiso_27 - Matrix checker 894 . -mat_mkl_cpardiso_31 - Partial solve and computing selected components of the solution vectors 895 . -mat_mkl_cpardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode 896 - -mat_mkl_cpardiso_60 - Intel MKL_CPARDISO mode 897 898 Level: beginner 899 900 Notes: 901 Use -mat_mkl_cpardiso_68 1 to display the number of threads the solver is using. MKL does not provide a way to directly access this 902 information. 903 904 For more information on the options check the MKL_CPARDISO manual 905 906 .seealso: PCFactorSetMatSolverType(), MatSolverType 907 908 M*/ 909 910 static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type) 911 { 912 PetscFunctionBegin; 913 *type = MATSOLVERMKL_CPARDISO; 914 PetscFunctionReturn(0); 915 } 916 917 /* MatGetFactor for MPI AIJ matrices */ 918 static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A,MatFactorType ftype,Mat *F) 919 { 920 Mat B; 921 PetscErrorCode ierr; 922 Mat_MKL_CPARDISO *mat_mkl_cpardiso; 923 PetscBool isSeqAIJ,isMPIBAIJ,isMPISBAIJ; 924 925 PetscFunctionBegin; 926 /* Create the factorization matrix */ 927 928 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 929 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&isMPIBAIJ);CHKERRQ(ierr); 930 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMPISBAIJ);CHKERRQ(ierr); 931 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 932 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 933 ierr = PetscStrallocpy("mkl_cpardiso",&((PetscObject)B)->type_name);CHKERRQ(ierr); 934 ierr = MatSetUp(B);CHKERRQ(ierr); 935 936 ierr = PetscNewLog(B,&mat_mkl_cpardiso);CHKERRQ(ierr); 937 938 if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO; 939 else if (isMPIBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO; 940 else if (isMPISBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO; 941 else mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO; 942 943 if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO; 944 else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_CPARDISO; 945 B->ops->destroy = MatDestroy_MKL_CPARDISO; 946 947 B->ops->view = MatView_MKL_CPARDISO; 948 B->ops->getinfo = MatGetInfo_MKL_CPARDISO; 949 950 B->factortype = ftype; 951 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 952 953 B->data = mat_mkl_cpardiso; 954 955 /* set solvertype */ 956 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 957 ierr = PetscStrallocpy(MATSOLVERMKL_CPARDISO,&B->solvertype);CHKERRQ(ierr); 958 959 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_cpardiso);CHKERRQ(ierr); 960 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMkl_CPardisoSetCntl_C",MatMkl_CPardisoSetCntl_MKL_CPARDISO);CHKERRQ(ierr); 961 ierr = PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso);CHKERRQ(ierr); 962 963 *F = B; 964 PetscFunctionReturn(0); 965 } 966 967 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void) 968 { 969 PetscErrorCode ierr; 970 971 PetscFunctionBegin; 972 ierr = MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);CHKERRQ(ierr); 973 ierr = MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);CHKERRQ(ierr); 974 ierr = MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIBAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);CHKERRQ(ierr); 975 ierr = MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPISBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_mpiaij_mkl_cpardiso);CHKERRQ(ierr); 976 PetscFunctionReturn(0); 977 } 978