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