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