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_Comm 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 PetscErrorCode ierr; 320 321 PetscFunctionBegin; 322 /* Terminate instance, deallocate memories */ 323 if (mat_mkl_cpardiso->CleanUp) { 324 mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY; 325 326 cluster_sparse_solver ( 327 mat_mkl_cpardiso->pt, 328 &mat_mkl_cpardiso->maxfct, 329 &mat_mkl_cpardiso->mnum, 330 &mat_mkl_cpardiso->mtype, 331 &mat_mkl_cpardiso->phase, 332 &mat_mkl_cpardiso->n, 333 NULL, 334 NULL, 335 NULL, 336 mat_mkl_cpardiso->perm, 337 &mat_mkl_cpardiso->nrhs, 338 mat_mkl_cpardiso->iparm, 339 &mat_mkl_cpardiso->msglvl, 340 NULL, 341 NULL, 342 &mat_mkl_cpardiso->comm_mkl_cpardiso, 343 (PetscInt*)&mat_mkl_cpardiso->err); 344 } 345 346 if (mat_mkl_cpardiso->ConvertToTriples != MatCopy_seqaij_seqaij_MKL_CPARDISO) { 347 ierr = PetscFree3(mat_mkl_cpardiso->ia,mat_mkl_cpardiso->ja,mat_mkl_cpardiso->a);CHKERRQ(ierr); 348 } 349 ierr = MPI_Comm_free(&(mat_mkl_cpardiso->comm_mkl_cpardiso));CHKERRQ(ierr); 350 ierr = PetscFree(A->data);CHKERRQ(ierr); 351 352 /* clear composed functions */ 353 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);CHKERRQ(ierr); 354 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMkl_CPardisoSetCntl_C",NULL);CHKERRQ(ierr); 355 PetscFunctionReturn(0); 356 } 357 358 /* 359 * Computes Ax = b 360 */ 361 PetscErrorCode MatSolve_MKL_CPARDISO(Mat A,Vec b,Vec x) 362 { 363 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data; 364 PetscErrorCode ierr; 365 PetscScalar *xarray; 366 const PetscScalar *barray; 367 368 PetscFunctionBegin; 369 mat_mkl_cpardiso->nrhs = 1; 370 ierr = VecGetArray(x,&xarray);CHKERRQ(ierr); 371 ierr = VecGetArrayRead(b,&barray);CHKERRQ(ierr); 372 373 /* solve phase */ 374 /*-------------*/ 375 mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 376 cluster_sparse_solver ( 377 mat_mkl_cpardiso->pt, 378 &mat_mkl_cpardiso->maxfct, 379 &mat_mkl_cpardiso->mnum, 380 &mat_mkl_cpardiso->mtype, 381 &mat_mkl_cpardiso->phase, 382 &mat_mkl_cpardiso->n, 383 mat_mkl_cpardiso->a, 384 mat_mkl_cpardiso->ia, 385 mat_mkl_cpardiso->ja, 386 mat_mkl_cpardiso->perm, 387 &mat_mkl_cpardiso->nrhs, 388 mat_mkl_cpardiso->iparm, 389 &mat_mkl_cpardiso->msglvl, 390 (void*)barray, 391 (void*)xarray, 392 &mat_mkl_cpardiso->comm_mkl_cpardiso, 393 (PetscInt*)&mat_mkl_cpardiso->err); 394 395 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)); 396 397 ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr); 398 ierr = VecRestoreArrayRead(b,&barray);CHKERRQ(ierr); 399 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 400 PetscFunctionReturn(0); 401 } 402 403 PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A,Vec b,Vec x) 404 { 405 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data; 406 PetscErrorCode ierr; 407 408 PetscFunctionBegin; 409 #if defined(PETSC_USE_COMPLEX) 410 mat_mkl_cpardiso->iparm[12 - 1] = 1; 411 #else 412 mat_mkl_cpardiso->iparm[12 - 1] = 2; 413 #endif 414 ierr = MatSolve_MKL_CPARDISO(A,b,x);CHKERRQ(ierr); 415 mat_mkl_cpardiso->iparm[12 - 1] = 0; 416 PetscFunctionReturn(0); 417 } 418 419 PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A,Mat B,Mat X) 420 { 421 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(A)->data; 422 PetscErrorCode ierr; 423 PetscScalar *xarray; 424 const PetscScalar *barray; 425 426 PetscFunctionBegin; 427 ierr = MatGetSize(B,NULL,(PetscInt*)&mat_mkl_cpardiso->nrhs);CHKERRQ(ierr); 428 429 if (mat_mkl_cpardiso->nrhs > 0) { 430 ierr = MatDenseGetArrayRead(B,&barray); 431 ierr = MatDenseGetArray(X,&xarray); 432 433 /* solve phase */ 434 /*-------------*/ 435 mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 436 cluster_sparse_solver ( 437 mat_mkl_cpardiso->pt, 438 &mat_mkl_cpardiso->maxfct, 439 &mat_mkl_cpardiso->mnum, 440 &mat_mkl_cpardiso->mtype, 441 &mat_mkl_cpardiso->phase, 442 &mat_mkl_cpardiso->n, 443 mat_mkl_cpardiso->a, 444 mat_mkl_cpardiso->ia, 445 mat_mkl_cpardiso->ja, 446 mat_mkl_cpardiso->perm, 447 &mat_mkl_cpardiso->nrhs, 448 mat_mkl_cpardiso->iparm, 449 &mat_mkl_cpardiso->msglvl, 450 (void*)barray, 451 (void*)xarray, 452 &mat_mkl_cpardiso->comm_mkl_cpardiso, 453 (PetscInt*)&mat_mkl_cpardiso->err); 454 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)); 455 ierr = MatDenseRestoreArrayRead(B,&barray); 456 ierr = MatDenseRestoreArray(X,&xarray); 457 458 } 459 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 460 PetscFunctionReturn(0); 461 } 462 463 /* 464 * LU Decomposition 465 */ 466 PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F,Mat A,const MatFactorInfo *info) 467 { 468 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)(F)->data; 469 PetscErrorCode ierr; 470 471 PetscFunctionBegin; 472 mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN; 473 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); 474 475 mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION; 476 cluster_sparse_solver ( 477 mat_mkl_cpardiso->pt, 478 &mat_mkl_cpardiso->maxfct, 479 &mat_mkl_cpardiso->mnum, 480 &mat_mkl_cpardiso->mtype, 481 &mat_mkl_cpardiso->phase, 482 &mat_mkl_cpardiso->n, 483 mat_mkl_cpardiso->a, 484 mat_mkl_cpardiso->ia, 485 mat_mkl_cpardiso->ja, 486 mat_mkl_cpardiso->perm, 487 &mat_mkl_cpardiso->nrhs, 488 mat_mkl_cpardiso->iparm, 489 &mat_mkl_cpardiso->msglvl, 490 NULL, 491 NULL, 492 &mat_mkl_cpardiso->comm_mkl_cpardiso, 493 &mat_mkl_cpardiso->err); 494 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)); 495 496 mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN; 497 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 498 PetscFunctionReturn(0); 499 } 500 501 /* Sets mkl_cpardiso options from the options database */ 502 PetscErrorCode PetscSetMKL_CPARDISOFromOptions(Mat F, Mat A) 503 { 504 Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data; 505 PetscErrorCode ierr; 506 PetscInt icntl,threads; 507 PetscBool flg; 508 509 PetscFunctionBegin; 510 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_CPARDISO Options","Mat");CHKERRQ(ierr); 511 ierr = PetscOptionsInt("-mat_mkl_cpardiso_65","Number of threads to use","None",threads,&threads,&flg);CHKERRQ(ierr); 512 if (flg) mkl_set_num_threads((int)threads); 513 514 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); 515 if (flg) mat_mkl_cpardiso->maxfct = icntl; 516 517 ierr = PetscOptionsInt("-mat_mkl_cpardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_cpardiso->mnum,&icntl,&flg);CHKERRQ(ierr); 518 if (flg) mat_mkl_cpardiso->mnum = icntl; 519 520 ierr = PetscOptionsInt("-mat_mkl_cpardiso_68","Message level information","None",mat_mkl_cpardiso->msglvl,&icntl,&flg);CHKERRQ(ierr); 521 if (flg) mat_mkl_cpardiso->msglvl = icntl; 522 523 ierr = PetscOptionsInt("-mat_mkl_cpardiso_69","Defines the matrix type","None",mat_mkl_cpardiso->mtype,&icntl,&flg);CHKERRQ(ierr); 524 if (flg) { 525 mat_mkl_cpardiso->mtype = icntl; 526 #if defined(PETSC_USE_REAL_SINGLE) 527 mat_mkl_cpardiso->iparm[27] = 1; 528 #else 529 mat_mkl_cpardiso->iparm[27] = 0; 530 #endif 531 mat_mkl_cpardiso->iparm[34] = 1; 532 } 533 ierr = PetscOptionsInt("-mat_mkl_cpardiso_1","Use default values","None",mat_mkl_cpardiso->iparm[0],&icntl,&flg);CHKERRQ(ierr); 534 535 if (flg && icntl != 0) { 536 ierr = PetscOptionsInt("-mat_mkl_cpardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_cpardiso->iparm[1],&icntl,&flg);CHKERRQ(ierr); 537 if (flg) mat_mkl_cpardiso->iparm[1] = icntl; 538 539 ierr = PetscOptionsInt("-mat_mkl_cpardiso_4","Preconditioned CGS/CG","None",mat_mkl_cpardiso->iparm[3],&icntl,&flg);CHKERRQ(ierr); 540 if (flg) mat_mkl_cpardiso->iparm[3] = icntl; 541 542 ierr = PetscOptionsInt("-mat_mkl_cpardiso_5","User permutation","None",mat_mkl_cpardiso->iparm[4],&icntl,&flg);CHKERRQ(ierr); 543 if (flg) mat_mkl_cpardiso->iparm[4] = icntl; 544 545 ierr = PetscOptionsInt("-mat_mkl_cpardiso_6","Write solution on x","None",mat_mkl_cpardiso->iparm[5],&icntl,&flg);CHKERRQ(ierr); 546 if (flg) mat_mkl_cpardiso->iparm[5] = icntl; 547 548 ierr = PetscOptionsInt("-mat_mkl_cpardiso_8","Iterative refinement step","None",mat_mkl_cpardiso->iparm[7],&icntl,&flg);CHKERRQ(ierr); 549 if (flg) mat_mkl_cpardiso->iparm[7] = icntl; 550 551 ierr = PetscOptionsInt("-mat_mkl_cpardiso_10","Pivoting perturbation","None",mat_mkl_cpardiso->iparm[9],&icntl,&flg);CHKERRQ(ierr); 552 if (flg) mat_mkl_cpardiso->iparm[9] = icntl; 553 554 ierr = PetscOptionsInt("-mat_mkl_cpardiso_11","Scaling vectors","None",mat_mkl_cpardiso->iparm[10],&icntl,&flg);CHKERRQ(ierr); 555 if (flg) mat_mkl_cpardiso->iparm[10] = icntl; 556 557 ierr = PetscOptionsInt("-mat_mkl_cpardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_cpardiso->iparm[11],&icntl,&flg);CHKERRQ(ierr); 558 if (flg) mat_mkl_cpardiso->iparm[11] = icntl; 559 560 ierr = PetscOptionsInt("-mat_mkl_cpardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_cpardiso->iparm[12],&icntl, 561 &flg);CHKERRQ(ierr); 562 if (flg) mat_mkl_cpardiso->iparm[12] = icntl; 563 564 ierr = PetscOptionsInt("-mat_mkl_cpardiso_18","Numbers of non-zero elements","None",mat_mkl_cpardiso->iparm[17],&icntl, 565 &flg);CHKERRQ(ierr); 566 if (flg) mat_mkl_cpardiso->iparm[17] = icntl; 567 568 ierr = PetscOptionsInt("-mat_mkl_cpardiso_19","Report number of floating point operations","None",mat_mkl_cpardiso->iparm[18],&icntl,&flg);CHKERRQ(ierr); 569 if (flg) mat_mkl_cpardiso->iparm[18] = icntl; 570 571 ierr = PetscOptionsInt("-mat_mkl_cpardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_cpardiso->iparm[20],&icntl,&flg);CHKERRQ(ierr); 572 if (flg) mat_mkl_cpardiso->iparm[20] = icntl; 573 574 ierr = PetscOptionsInt("-mat_mkl_cpardiso_24","Parallel factorization control","None",mat_mkl_cpardiso->iparm[23],&icntl,&flg);CHKERRQ(ierr); 575 if (flg) mat_mkl_cpardiso->iparm[23] = icntl; 576 577 ierr = PetscOptionsInt("-mat_mkl_cpardiso_25","Parallel forward/backward solve control","None",mat_mkl_cpardiso->iparm[24],&icntl,&flg);CHKERRQ(ierr); 578 if (flg) mat_mkl_cpardiso->iparm[24] = icntl; 579 580 ierr = PetscOptionsInt("-mat_mkl_cpardiso_27","Matrix checker","None",mat_mkl_cpardiso->iparm[26],&icntl,&flg);CHKERRQ(ierr); 581 if (flg) mat_mkl_cpardiso->iparm[26] = icntl; 582 583 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); 584 if (flg) mat_mkl_cpardiso->iparm[30] = icntl; 585 586 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); 587 if (flg) mat_mkl_cpardiso->iparm[33] = icntl; 588 589 ierr = PetscOptionsInt("-mat_mkl_cpardiso_60","Intel MKL_CPARDISO mode","None",mat_mkl_cpardiso->iparm[59],&icntl,&flg);CHKERRQ(ierr); 590 if (flg) mat_mkl_cpardiso->iparm[59] = icntl; 591 } 592 593 ierr = PetscOptionsEnd();CHKERRQ(ierr); 594 PetscFunctionReturn(0); 595 } 596 597 PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso) 598 { 599 PetscErrorCode ierr; 600 PetscInt bs; 601 PetscBool match; 602 PetscMPIInt size; 603 604 PetscFunctionBegin; 605 606 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(mat_mkl_cpardiso->comm_mkl_cpardiso));CHKERRQ(ierr); 607 ierr = MPI_Comm_size(mat_mkl_cpardiso->comm_mkl_cpardiso, &size);CHKERRQ(ierr); 608 609 mat_mkl_cpardiso->CleanUp = PETSC_FALSE; 610 mat_mkl_cpardiso->maxfct = 1; 611 mat_mkl_cpardiso->mnum = 1; 612 mat_mkl_cpardiso->n = A->rmap->N; 613 if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36]; 614 mat_mkl_cpardiso->msglvl = 0; 615 mat_mkl_cpardiso->nrhs = 1; 616 mat_mkl_cpardiso->err = 0; 617 mat_mkl_cpardiso->phase = -1; 618 #if defined(PETSC_USE_COMPLEX) 619 mat_mkl_cpardiso->mtype = 13; 620 #else 621 mat_mkl_cpardiso->mtype = 11; 622 #endif 623 624 #if defined(PETSC_USE_REAL_SINGLE) 625 mat_mkl_cpardiso->iparm[27] = 1; 626 #else 627 mat_mkl_cpardiso->iparm[27] = 0; 628 #endif 629 630 mat_mkl_cpardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */ 631 mat_mkl_cpardiso->iparm[ 1] = 2; /* Use METIS for fill-in reordering */ 632 mat_mkl_cpardiso->iparm[ 5] = 0; /* Write solution into x */ 633 mat_mkl_cpardiso->iparm[ 7] = 2; /* Max number of iterative refinement steps */ 634 mat_mkl_cpardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */ 635 mat_mkl_cpardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ 636 mat_mkl_cpardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */ 637 mat_mkl_cpardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */ 638 mat_mkl_cpardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */ 639 mat_mkl_cpardiso->iparm[26] = 1; /* Check input data for correctness */ 640 641 mat_mkl_cpardiso->iparm[39] = 0; 642 if (size > 1) { 643 mat_mkl_cpardiso->iparm[39] = 2; 644 mat_mkl_cpardiso->iparm[40] = A->rmap->rstart; 645 mat_mkl_cpardiso->iparm[41] = A->rmap->rend-1; 646 } 647 ierr = PetscObjectTypeCompareAny((PetscObject)A,&match,MATMPIBAIJ,MATMPISBAIJ,"");CHKERRQ(ierr); 648 if (match) { 649 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 650 mat_mkl_cpardiso->iparm[36] = bs; 651 mat_mkl_cpardiso->iparm[40] /= bs; 652 mat_mkl_cpardiso->iparm[41] /= bs; 653 mat_mkl_cpardiso->iparm[40]++; 654 mat_mkl_cpardiso->iparm[41]++; 655 mat_mkl_cpardiso->iparm[34] = 0; /* Fortran style */ 656 } else { 657 mat_mkl_cpardiso->iparm[34] = 1; /* C style */ 658 } 659 660 mat_mkl_cpardiso->perm = 0; 661 PetscFunctionReturn(0); 662 } 663 664 /* 665 * Symbolic decomposition. Mkl_Pardiso analysis phase. 666 */ 667 PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 668 { 669 Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data; 670 PetscErrorCode ierr; 671 672 PetscFunctionBegin; 673 mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN; 674 675 /* Set MKL_CPARDISO options from the options database */ 676 ierr = PetscSetMKL_CPARDISOFromOptions(F,A);CHKERRQ(ierr); 677 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); 678 679 mat_mkl_cpardiso->n = A->rmap->N; 680 if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36]; 681 682 /* analysis phase */ 683 /*----------------*/ 684 mat_mkl_cpardiso->phase = JOB_ANALYSIS; 685 686 cluster_sparse_solver ( 687 mat_mkl_cpardiso->pt, 688 &mat_mkl_cpardiso->maxfct, 689 &mat_mkl_cpardiso->mnum, 690 &mat_mkl_cpardiso->mtype, 691 &mat_mkl_cpardiso->phase, 692 &mat_mkl_cpardiso->n, 693 mat_mkl_cpardiso->a, 694 mat_mkl_cpardiso->ia, 695 mat_mkl_cpardiso->ja, 696 mat_mkl_cpardiso->perm, 697 &mat_mkl_cpardiso->nrhs, 698 mat_mkl_cpardiso->iparm, 699 &mat_mkl_cpardiso->msglvl, 700 NULL, 701 NULL, 702 &mat_mkl_cpardiso->comm_mkl_cpardiso, 703 (PetscInt*)&mat_mkl_cpardiso->err); 704 705 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)); 706 707 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 708 F->ops->lufactornumeric = MatFactorNumeric_MKL_CPARDISO; 709 F->ops->solve = MatSolve_MKL_CPARDISO; 710 F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO; 711 F->ops->matsolve = MatMatSolve_MKL_CPARDISO; 712 PetscFunctionReturn(0); 713 } 714 715 PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_CPARDISO(Mat F,Mat A,IS perm,const MatFactorInfo *info) 716 { 717 Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO*)F->data; 718 PetscErrorCode ierr; 719 720 PetscFunctionBegin; 721 mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN; 722 723 724 /* Set MKL_CPARDISO options from the options database */ 725 ierr = PetscSetMKL_CPARDISOFromOptions(F,A);CHKERRQ(ierr); 726 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); 727 728 mat_mkl_cpardiso->n = A->rmap->N; 729 if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36]; 730 #if defined(PETSC_USE_COMPLEX) 731 SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead",((PetscObject)A)->type_name); 732 #endif 733 if (A->spd_set && A->spd) mat_mkl_cpardiso->mtype = 2; 734 else mat_mkl_cpardiso->mtype = -2; 735 736 /* analysis phase */ 737 /*----------------*/ 738 mat_mkl_cpardiso->phase = JOB_ANALYSIS; 739 740 cluster_sparse_solver ( 741 mat_mkl_cpardiso->pt, 742 &mat_mkl_cpardiso->maxfct, 743 &mat_mkl_cpardiso->mnum, 744 &mat_mkl_cpardiso->mtype, 745 &mat_mkl_cpardiso->phase, 746 &mat_mkl_cpardiso->n, 747 mat_mkl_cpardiso->a, 748 mat_mkl_cpardiso->ia, 749 mat_mkl_cpardiso->ja, 750 mat_mkl_cpardiso->perm, 751 &mat_mkl_cpardiso->nrhs, 752 mat_mkl_cpardiso->iparm, 753 &mat_mkl_cpardiso->msglvl, 754 NULL, 755 NULL, 756 &mat_mkl_cpardiso->comm_mkl_cpardiso, 757 (PetscInt*)&mat_mkl_cpardiso->err); 758 759 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)); 760 761 mat_mkl_cpardiso->CleanUp = PETSC_TRUE; 762 F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_CPARDISO; 763 F->ops->solve = MatSolve_MKL_CPARDISO; 764 F->ops->solvetranspose = MatSolveTranspose_MKL_CPARDISO; 765 F->ops->matsolve = MatMatSolve_MKL_CPARDISO; 766 PetscFunctionReturn(0); 767 } 768 769 PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer) 770 { 771 PetscErrorCode ierr; 772 PetscBool iascii; 773 PetscViewerFormat format; 774 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data; 775 PetscInt i; 776 777 PetscFunctionBegin; 778 /* check if matrix is mkl_cpardiso type */ 779 if (A->ops->solve != MatSolve_MKL_CPARDISO) PetscFunctionReturn(0); 780 781 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 782 if (iascii) { 783 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 784 if (format == PETSC_VIEWER_ASCII_INFO) { 785 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO run parameters:\n");CHKERRQ(ierr); 786 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO phase: %d \n",mat_mkl_cpardiso->phase);CHKERRQ(ierr); 787 for (i = 1; i <= 64; i++) { 788 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO iparm[%d]: %d \n",i, mat_mkl_cpardiso->iparm[i - 1]);CHKERRQ(ierr); 789 } 790 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO maxfct: %d \n", mat_mkl_cpardiso->maxfct);CHKERRQ(ierr); 791 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mnum: %d \n", mat_mkl_cpardiso->mnum);CHKERRQ(ierr); 792 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO mtype: %d \n", mat_mkl_cpardiso->mtype);CHKERRQ(ierr); 793 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO n: %d \n", mat_mkl_cpardiso->n);CHKERRQ(ierr); 794 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO nrhs: %d \n", mat_mkl_cpardiso->nrhs);CHKERRQ(ierr); 795 ierr = PetscViewerASCIIPrintf(viewer,"MKL_CPARDISO msglvl: %d \n", mat_mkl_cpardiso->msglvl);CHKERRQ(ierr); 796 } 797 } 798 PetscFunctionReturn(0); 799 } 800 801 PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info) 802 { 803 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)A->data; 804 805 PetscFunctionBegin; 806 info->block_size = 1.0; 807 info->nz_allocated = mat_mkl_cpardiso->nz + 0.0; 808 info->nz_unneeded = 0.0; 809 info->assemblies = 0.0; 810 info->mallocs = 0.0; 811 info->memory = 0.0; 812 info->fill_ratio_given = 0; 813 info->fill_ratio_needed = 0; 814 info->factor_mallocs = 0; 815 PetscFunctionReturn(0); 816 } 817 818 PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F,PetscInt icntl,PetscInt ival) 819 { 820 Mat_MKL_CPARDISO *mat_mkl_cpardiso=(Mat_MKL_CPARDISO*)F->data; 821 822 PetscFunctionBegin; 823 if (icntl <= 64) { 824 mat_mkl_cpardiso->iparm[icntl - 1] = ival; 825 } else { 826 if (icntl == 65) mkl_set_num_threads((int)ival); 827 else if (icntl == 66) mat_mkl_cpardiso->maxfct = ival; 828 else if (icntl == 67) mat_mkl_cpardiso->mnum = ival; 829 else if (icntl == 68) mat_mkl_cpardiso->msglvl = ival; 830 else if (icntl == 69) { 831 mat_mkl_cpardiso->mtype = ival; 832 #if defined(PETSC_USE_REAL_SINGLE) 833 mat_mkl_cpardiso->iparm[27] = 1; 834 #else 835 mat_mkl_cpardiso->iparm[27] = 0; 836 #endif 837 mat_mkl_cpardiso->iparm[34] = 1; 838 } 839 } 840 PetscFunctionReturn(0); 841 } 842 843 /*@ 844 MatMkl_CPardisoSetCntl - Set Mkl_Pardiso parameters 845 846 Logically Collective on Mat 847 848 Input Parameters: 849 + F - the factored matrix obtained by calling MatGetFactor() 850 . icntl - index of Mkl_Pardiso parameter 851 - ival - value of Mkl_Pardiso parameter 852 853 Options Database: 854 . -mat_mkl_cpardiso_<icntl> <ival> 855 856 Level: Intermediate 857 858 Notes: 859 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 860 database approach when working with TS, SNES, or KSP. 861 862 References: 863 . Mkl_Pardiso Users' Guide 864 865 .seealso: MatGetFactor() 866 @*/ 867 PetscErrorCode MatMkl_CPardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival) 868 { 869 PetscErrorCode ierr; 870 871 PetscFunctionBegin; 872 ierr = PetscTryMethod(F,"MatMkl_CPardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 873 PetscFunctionReturn(0); 874 } 875 876 static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type) 877 { 878 PetscFunctionBegin; 879 *type = MATSOLVERMKL_CPARDISO; 880 PetscFunctionReturn(0); 881 } 882 883 /* MatGetFactor for MPI AIJ matrices */ 884 static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A,MatFactorType ftype,Mat *F) 885 { 886 Mat B; 887 PetscErrorCode ierr; 888 Mat_MKL_CPARDISO *mat_mkl_cpardiso; 889 PetscBool isSeqAIJ,isMPIBAIJ,isMPISBAIJ; 890 891 PetscFunctionBegin; 892 /* Create the factorization matrix */ 893 894 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 895 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&isMPIBAIJ);CHKERRQ(ierr); 896 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPISBAIJ,&isMPISBAIJ);CHKERRQ(ierr); 897 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 898 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 899 ierr = PetscStrallocpy("mkl_cpardiso",&((PetscObject)B)->type_name);CHKERRQ(ierr); 900 ierr = MatSetUp(B);CHKERRQ(ierr); 901 902 ierr = PetscNewLog(B,&mat_mkl_cpardiso);CHKERRQ(ierr); 903 904 if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO; 905 else if (isMPIBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO; 906 else if (isMPISBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO; 907 else mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO; 908 909 if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO; 910 else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_CPARDISO; 911 B->ops->destroy = MatDestroy_MKL_CPARDISO; 912 913 B->ops->view = MatView_MKL_CPARDISO; 914 B->ops->getinfo = MatGetInfo_MKL_CPARDISO; 915 916 B->factortype = ftype; 917 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 918 919 B->data = mat_mkl_cpardiso; 920 921 /* set solvertype */ 922 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 923 ierr = PetscStrallocpy(MATSOLVERMKL_CPARDISO,&B->solvertype);CHKERRQ(ierr); 924 925 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_cpardiso);CHKERRQ(ierr); 926 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMkl_CPardisoSetCntl_C",MatMkl_CPardisoSetCntl_MKL_CPARDISO);CHKERRQ(ierr); 927 ierr = PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso);CHKERRQ(ierr); 928 929 *F = B; 930 PetscFunctionReturn(0); 931 } 932 933 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void) 934 { 935 PetscErrorCode ierr; 936 937 PetscFunctionBegin; 938 ierr = MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);CHKERRQ(ierr); 939 ierr = MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);CHKERRQ(ierr); 940 ierr = MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPIBAIJ,MAT_FACTOR_LU,MatGetFactor_mpiaij_mkl_cpardiso);CHKERRQ(ierr); 941 ierr = MatSolverTypeRegister(MATSOLVERMKL_CPARDISO,MATMPISBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_mpiaij_mkl_cpardiso);CHKERRQ(ierr); 942 PetscFunctionReturn(0); 943 } 944