1 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 2 #include <../src/mat/impls/sbaij/seq/sbaij.h> 3 #include <../src/mat/impls/dense/seq/dense.h> 4 5 #if defined(PETSC_HAVE_MKL_INTEL_ILP64) 6 #define MKL_ILP64 7 #endif 8 #include <mkl_pardiso.h> 9 10 PETSC_EXTERN void PetscSetMKL_PARDISOThreads(int); 11 12 /* 13 * Possible mkl_pardiso phases that controls the execution of the solver. 14 * For more information check mkl_pardiso 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 30 #if defined(PETSC_USE_64BIT_INDICES) 31 #if defined(PETSC_HAVE_MKL_INTEL_ILP64) 32 #define INT_TYPE long long int 33 #define MKL_PARDISO pardiso 34 #define MKL_PARDISO_INIT pardisoinit 35 #else 36 /* this is the case where the MKL BLAS/LAPACK are 32 bit integers but the 64 bit integer version of 37 of Pardiso code is used; hence the need for the 64 below*/ 38 #define INT_TYPE long long int 39 #define MKL_PARDISO pardiso_64 40 #define MKL_PARDISO_INIT pardiso_64init 41 void pardiso_64init(void *pt, INT_TYPE *mtype, INT_TYPE iparm []) 42 { 43 int iparm_copy[IPARM_SIZE], mtype_copy, i; 44 45 mtype_copy = *mtype; 46 pardisoinit(pt, &mtype_copy, iparm_copy); 47 for(i = 0; i < IPARM_SIZE; i++){ 48 iparm[i] = iparm_copy[i]; 49 } 50 } 51 #endif 52 #else 53 #define INT_TYPE int 54 #define MKL_PARDISO pardiso 55 #define MKL_PARDISO_INIT pardisoinit 56 #endif 57 58 59 /* 60 * Internal data structure. 61 * For more information check mkl_pardiso manual. 62 */ 63 typedef struct { 64 65 /* Configuration vector*/ 66 INT_TYPE iparm[IPARM_SIZE]; 67 68 /* 69 * Internal mkl_pardiso memory location. 70 * After the first call to mkl_pardiso do not modify pt, as that could cause a serious memory leak. 71 */ 72 void *pt[IPARM_SIZE]; 73 74 /* Basic mkl_pardiso info*/ 75 INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err; 76 77 /* Matrix structure*/ 78 void *a; 79 INT_TYPE *ia, *ja; 80 81 /* Number of non-zero elements*/ 82 INT_TYPE nz; 83 84 /* Row permutaton vector*/ 85 INT_TYPE *perm; 86 87 /* Define if matrix preserves sparse structure.*/ 88 MatStructure matstruc; 89 90 PetscBool needsym; 91 PetscBool freeaij; 92 93 /* Schur complement */ 94 PetscScalar *schur; 95 PetscInt schur_size; 96 PetscInt *schur_idxs; 97 PetscScalar *schur_work; 98 PetscBLASInt schur_work_size; 99 PetscBool solve_interior; 100 101 /* True if mkl_pardiso function have been used.*/ 102 PetscBool CleanUp; 103 104 /* Conversion to a format suitable for MKL */ 105 PetscErrorCode (*Convert)(Mat, PetscBool, MatReuse, PetscBool*, INT_TYPE*, INT_TYPE**, INT_TYPE**, PetscScalar**); 106 } Mat_MKL_PARDISO; 107 108 PetscErrorCode MatMKLPardiso_Convert_seqsbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v) 109 { 110 Mat_SeqSBAIJ *aa = (Mat_SeqSBAIJ*)A->data; 111 PetscInt bs = A->rmap->bs; 112 113 PetscFunctionBegin; 114 if (!sym) { 115 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen"); 116 } 117 if (bs == 1) { /* already in the correct format */ 118 *v = aa->a; 119 /* though PetscInt and INT_TYPE are of the same size since they are defined differently the Intel compiler requires a cast */ 120 *r = (INT_TYPE*)aa->i; 121 *c = (INT_TYPE*)aa->j; 122 *nnz = (INT_TYPE)aa->nz; 123 *free = PETSC_FALSE; 124 } else { 125 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Conversion from SeqSBAIJ to MKL Pardiso format still need to be implemented"); 126 } 127 PetscFunctionReturn(0); 128 } 129 130 PetscErrorCode MatMKLPardiso_Convert_seqbaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v) 131 { 132 PetscFunctionBegin; 133 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Conversion from SeqBAIJ to MKL Pardiso format still need to be implemented"); 134 PetscFunctionReturn(0); 135 } 136 137 PetscErrorCode MatMKLPardiso_Convert_seqaij(Mat A,PetscBool sym,MatReuse reuse,PetscBool *free,INT_TYPE *nnz,INT_TYPE **r,INT_TYPE **c,PetscScalar **v) 138 { 139 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)A->data; 140 PetscErrorCode ierr; 141 142 PetscFunctionBegin; 143 if (!sym) { /* already in the correct format */ 144 *v = aa->a; 145 *r = (INT_TYPE*)aa->i; 146 *c = (INT_TYPE*)aa->j; 147 *nnz = (INT_TYPE)aa->nz; 148 *free = PETSC_FALSE; 149 PetscFunctionReturn(0); 150 } 151 /* need to get the triangular part */ 152 if (reuse == MAT_INITIAL_MATRIX) { 153 PetscScalar *vals,*vv; 154 PetscInt *row,*col,*jj; 155 PetscInt m = A->rmap->n,nz,i; 156 157 nz = 0; 158 for (i=0; i<m; i++) { 159 nz += aa->i[i+1] - aa->diag[i]; 160 } 161 ierr = PetscMalloc2(m+1,&row,nz,&col);CHKERRQ(ierr); 162 ierr = PetscMalloc1(nz,&vals);CHKERRQ(ierr); 163 jj = col; 164 vv = vals; 165 166 row[0] = 0; 167 for (i=0; i<m; i++) { 168 PetscInt *aj = aa->j + aa->diag[i]; 169 PetscScalar *av = aa->a + aa->diag[i]; 170 PetscInt rl = aa->i[i+1] - aa->diag[i],j; 171 for (j=0; j<rl; j++) { 172 *jj = *aj; jj++; aj++; 173 *vv = *av; vv++; av++; 174 } 175 row[i+1] = row[i] + rl; 176 } 177 *v = vals; 178 *r = (INT_TYPE*)row; 179 *c = (INT_TYPE*)col; 180 *nnz = (INT_TYPE)nz; 181 *free = PETSC_TRUE; 182 } else { 183 PetscScalar *vv; 184 PetscInt m = A->rmap->n,i; 185 186 vv = *v; 187 for (i=0; i<m; i++) { 188 PetscScalar *av = aa->a + aa->diag[i]; 189 PetscInt rl = aa->i[i+1] - aa->diag[i],j; 190 for (j=0; j<rl; j++) { 191 *vv = *av; vv++; av++; 192 } 193 } 194 *free = PETSC_TRUE; 195 } 196 PetscFunctionReturn(0); 197 } 198 199 200 static PetscErrorCode MatMKLPardisoSolveSchur_Private(Mat F, PetscScalar *B, PetscScalar *X) 201 { 202 Mat_MKL_PARDISO *mpardiso =(Mat_MKL_PARDISO*)F->data; 203 Mat S,Xmat,Bmat; 204 MatFactorSchurStatus schurstatus; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 ierr = MatFactorFactorizeSchurComplement(F);CHKERRQ(ierr); 209 ierr = MatFactorGetSchurComplement(F,&S,&schurstatus);CHKERRQ(ierr); 210 if (X == B && schurstatus == MAT_FACTOR_SCHUR_INVERTED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"X and B cannot point to the same address"); 211 ierr = MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,B,&Bmat);CHKERRQ(ierr); 212 ierr = MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->nrhs,X,&Xmat);CHKERRQ(ierr); 213 if (X != B) { /* using MatMatSolve */ 214 ierr = MatCopy(Bmat,Xmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 215 } 216 217 #if defined(PETSC_USE_COMPLEX) 218 if (mpardiso->iparm[12-1] == 1) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Hermitian solve not implemented yet"); 219 #endif 220 221 switch (schurstatus) { 222 case MAT_FACTOR_SCHUR_FACTORED: 223 if (!mpardiso->iparm[12-1]) { 224 ierr = MatMatSolve(S,Bmat,Xmat);CHKERRQ(ierr); 225 } else { /* transpose solve */ 226 ierr = MatMatSolveTranspose(S,Bmat,Xmat);CHKERRQ(ierr); 227 } 228 break; 229 case MAT_FACTOR_SCHUR_INVERTED: 230 if (!mpardiso->iparm[12-1]) { 231 ierr = MatMatMult(S,Bmat,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Xmat);CHKERRQ(ierr); 232 } else { /* transpose solve */ 233 ierr = MatTransposeMatMult(S,Bmat,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Xmat);CHKERRQ(ierr); 234 } 235 break; 236 default: 237 SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Unhandled MatFactorSchurStatus %D",F->schur_status); 238 break; 239 } 240 ierr = MatFactorRestoreSchurComplement(F,&S,schurstatus);CHKERRQ(ierr); 241 ierr = MatDestroy(&Bmat);CHKERRQ(ierr); 242 ierr = MatDestroy(&Xmat);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode MatFactorSetSchurIS_MKL_PARDISO(Mat F, IS is) 247 { 248 Mat_MKL_PARDISO *mpardiso =(Mat_MKL_PARDISO*)F->data; 249 const PetscInt *idxs; 250 PetscInt size,i; 251 PetscMPIInt csize; 252 PetscBool sorted; 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)F),&csize);CHKERRQ(ierr); 257 if (csize > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MKL_PARDISO parallel Schur complements not yet supported from PETSc"); 258 ierr = ISSorted(is,&sorted);CHKERRQ(ierr); 259 if (!sorted) { 260 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS for MKL_PARDISO Schur complements needs to be sorted"); 261 } 262 ierr = ISGetLocalSize(is,&size);CHKERRQ(ierr); 263 if (mpardiso->schur_size != size) { 264 mpardiso->schur_size = size; 265 ierr = PetscFree2(mpardiso->schur,mpardiso->schur_work);CHKERRQ(ierr); 266 ierr = PetscFree(mpardiso->schur_idxs);CHKERRQ(ierr); 267 ierr = PetscBLASIntCast(PetscMax(mpardiso->n,2*size),&mpardiso->schur_work_size);CHKERRQ(ierr); 268 ierr = PetscMalloc2(size*size,&mpardiso->schur,mpardiso->schur_work_size,&mpardiso->schur_work);CHKERRQ(ierr); 269 ierr = PetscMalloc1(size,&mpardiso->schur_idxs);CHKERRQ(ierr); 270 } 271 ierr = MatCreateSeqDense(PETSC_COMM_SELF,mpardiso->schur_size,mpardiso->schur_size,mpardiso->schur,&F->schur);CHKERRQ(ierr); 272 if (mpardiso->mtype == 2) { 273 ierr = MatSetOption(F->schur,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr); 274 } 275 276 ierr = PetscArrayzero(mpardiso->perm,mpardiso->n);CHKERRQ(ierr); 277 ierr = ISGetIndices(is,&idxs);CHKERRQ(ierr); 278 ierr = PetscArraycpy(mpardiso->schur_idxs,idxs,size);CHKERRQ(ierr); 279 for (i=0;i<size;i++) mpardiso->perm[idxs[i]] = 1; 280 ierr = ISRestoreIndices(is,&idxs);CHKERRQ(ierr); 281 if (size) { /* turn on Schur switch if the set of indices is not empty */ 282 mpardiso->iparm[36-1] = 2; 283 } 284 PetscFunctionReturn(0); 285 } 286 287 PetscErrorCode MatDestroy_MKL_PARDISO(Mat A) 288 { 289 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data; 290 PetscErrorCode ierr; 291 292 PetscFunctionBegin; 293 if (mat_mkl_pardiso->CleanUp) { 294 mat_mkl_pardiso->phase = JOB_RELEASE_OF_ALL_MEMORY; 295 296 MKL_PARDISO (mat_mkl_pardiso->pt, 297 &mat_mkl_pardiso->maxfct, 298 &mat_mkl_pardiso->mnum, 299 &mat_mkl_pardiso->mtype, 300 &mat_mkl_pardiso->phase, 301 &mat_mkl_pardiso->n, 302 NULL, 303 NULL, 304 NULL, 305 NULL, 306 &mat_mkl_pardiso->nrhs, 307 mat_mkl_pardiso->iparm, 308 &mat_mkl_pardiso->msglvl, 309 NULL, 310 NULL, 311 &mat_mkl_pardiso->err); 312 } 313 ierr = PetscFree(mat_mkl_pardiso->perm);CHKERRQ(ierr); 314 ierr = PetscFree2(mat_mkl_pardiso->schur,mat_mkl_pardiso->schur_work);CHKERRQ(ierr); 315 ierr = PetscFree(mat_mkl_pardiso->schur_idxs);CHKERRQ(ierr); 316 if (mat_mkl_pardiso->freeaij) { 317 ierr = PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);CHKERRQ(ierr); 318 ierr = PetscFree(mat_mkl_pardiso->a);CHKERRQ(ierr); 319 } 320 ierr = PetscFree(A->data);CHKERRQ(ierr); 321 322 /* clear composed functions */ 323 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);CHKERRQ(ierr); 324 ierr = PetscObjectComposeFunction((PetscObject)A,"MatFactorSetSchurIS_C",NULL);CHKERRQ(ierr); 325 ierr = PetscObjectComposeFunction((PetscObject)A,"MatMkl_PardisoSetCntl_C",NULL);CHKERRQ(ierr); 326 PetscFunctionReturn(0); 327 } 328 329 static PetscErrorCode MatMKLPardisoScatterSchur_Private(Mat_MKL_PARDISO *mpardiso, PetscScalar *whole, PetscScalar *schur, PetscBool reduce) 330 { 331 PetscFunctionBegin; 332 if (reduce) { /* data given for the whole matrix */ 333 PetscInt i,m=0,p=0; 334 for (i=0;i<mpardiso->nrhs;i++) { 335 PetscInt j; 336 for (j=0;j<mpardiso->schur_size;j++) { 337 schur[p+j] = whole[m+mpardiso->schur_idxs[j]]; 338 } 339 m += mpardiso->n; 340 p += mpardiso->schur_size; 341 } 342 } else { /* from Schur to whole */ 343 PetscInt i,m=0,p=0; 344 for (i=0;i<mpardiso->nrhs;i++) { 345 PetscInt j; 346 for (j=0;j<mpardiso->schur_size;j++) { 347 whole[m+mpardiso->schur_idxs[j]] = schur[p+j]; 348 } 349 m += mpardiso->n; 350 p += mpardiso->schur_size; 351 } 352 } 353 PetscFunctionReturn(0); 354 } 355 356 PetscErrorCode MatSolve_MKL_PARDISO(Mat A,Vec b,Vec x) 357 { 358 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data; 359 PetscErrorCode ierr; 360 PetscScalar *xarray; 361 const PetscScalar *barray; 362 363 PetscFunctionBegin; 364 mat_mkl_pardiso->nrhs = 1; 365 ierr = VecGetArray(x,&xarray);CHKERRQ(ierr); 366 ierr = VecGetArrayRead(b,&barray);CHKERRQ(ierr); 367 368 if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 369 else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION; 370 371 if (barray == xarray) { /* if the two vectors share the same memory */ 372 PetscScalar *work; 373 if (!mat_mkl_pardiso->schur_work) { 374 ierr = PetscMalloc1(mat_mkl_pardiso->n,&work);CHKERRQ(ierr); 375 } else { 376 work = mat_mkl_pardiso->schur_work; 377 } 378 mat_mkl_pardiso->iparm[6-1] = 1; 379 MKL_PARDISO (mat_mkl_pardiso->pt, 380 &mat_mkl_pardiso->maxfct, 381 &mat_mkl_pardiso->mnum, 382 &mat_mkl_pardiso->mtype, 383 &mat_mkl_pardiso->phase, 384 &mat_mkl_pardiso->n, 385 mat_mkl_pardiso->a, 386 mat_mkl_pardiso->ia, 387 mat_mkl_pardiso->ja, 388 NULL, 389 &mat_mkl_pardiso->nrhs, 390 mat_mkl_pardiso->iparm, 391 &mat_mkl_pardiso->msglvl, 392 (void*)xarray, 393 (void*)work, 394 &mat_mkl_pardiso->err); 395 if (!mat_mkl_pardiso->schur_work) { 396 ierr = PetscFree(work);CHKERRQ(ierr); 397 } 398 } else { 399 mat_mkl_pardiso->iparm[6-1] = 0; 400 MKL_PARDISO (mat_mkl_pardiso->pt, 401 &mat_mkl_pardiso->maxfct, 402 &mat_mkl_pardiso->mnum, 403 &mat_mkl_pardiso->mtype, 404 &mat_mkl_pardiso->phase, 405 &mat_mkl_pardiso->n, 406 mat_mkl_pardiso->a, 407 mat_mkl_pardiso->ia, 408 mat_mkl_pardiso->ja, 409 mat_mkl_pardiso->perm, 410 &mat_mkl_pardiso->nrhs, 411 mat_mkl_pardiso->iparm, 412 &mat_mkl_pardiso->msglvl, 413 (void*)barray, 414 (void*)xarray, 415 &mat_mkl_pardiso->err); 416 } 417 ierr = VecRestoreArrayRead(b,&barray);CHKERRQ(ierr); 418 419 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err); 420 421 if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */ 422 PetscInt shift = mat_mkl_pardiso->schur_size; 423 424 /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */ 425 if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) { 426 shift = 0; 427 } 428 429 if (!mat_mkl_pardiso->solve_interior) { 430 /* solve Schur complement */ 431 ierr = MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);CHKERRQ(ierr); 432 ierr = MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);CHKERRQ(ierr); 433 ierr = MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);CHKERRQ(ierr); 434 } else { /* if we are solving for the interior problem, any value in barray[schur] forward-substituted to xarray[schur] will be neglected */ 435 PetscInt i; 436 for (i=0;i<mat_mkl_pardiso->schur_size;i++) { 437 xarray[mat_mkl_pardiso->schur_idxs[i]] = 0.; 438 } 439 } 440 441 /* expansion phase */ 442 mat_mkl_pardiso->iparm[6-1] = 1; 443 mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION; 444 MKL_PARDISO (mat_mkl_pardiso->pt, 445 &mat_mkl_pardiso->maxfct, 446 &mat_mkl_pardiso->mnum, 447 &mat_mkl_pardiso->mtype, 448 &mat_mkl_pardiso->phase, 449 &mat_mkl_pardiso->n, 450 mat_mkl_pardiso->a, 451 mat_mkl_pardiso->ia, 452 mat_mkl_pardiso->ja, 453 mat_mkl_pardiso->perm, 454 &mat_mkl_pardiso->nrhs, 455 mat_mkl_pardiso->iparm, 456 &mat_mkl_pardiso->msglvl, 457 (void*)xarray, 458 (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */ 459 &mat_mkl_pardiso->err); 460 461 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err); 462 mat_mkl_pardiso->iparm[6-1] = 0; 463 } 464 ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr); 465 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 466 PetscFunctionReturn(0); 467 } 468 469 PetscErrorCode MatSolveTranspose_MKL_PARDISO(Mat A,Vec b,Vec x) 470 { 471 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data; 472 PetscInt oiparm12; 473 PetscErrorCode ierr; 474 475 PetscFunctionBegin; 476 oiparm12 = mat_mkl_pardiso->iparm[12 - 1]; 477 mat_mkl_pardiso->iparm[12 - 1] = 2; 478 ierr = MatSolve_MKL_PARDISO(A,b,x);CHKERRQ(ierr); 479 mat_mkl_pardiso->iparm[12 - 1] = oiparm12; 480 PetscFunctionReturn(0); 481 } 482 483 PetscErrorCode MatMatSolve_MKL_PARDISO(Mat A,Mat B,Mat X) 484 { 485 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(A)->data; 486 PetscErrorCode ierr; 487 const PetscScalar *barray; 488 PetscScalar *xarray; 489 PetscBool flg; 490 491 PetscFunctionBegin; 492 ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 493 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATSEQDENSE matrix"); 494 ierr = PetscObjectBaseTypeCompare((PetscObject)X,MATSEQDENSE,&flg);CHKERRQ(ierr); 495 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATSEQDENSE matrix"); 496 497 ierr = MatGetSize(B,NULL,(PetscInt*)&mat_mkl_pardiso->nrhs);CHKERRQ(ierr); 498 499 if (mat_mkl_pardiso->nrhs > 0) { 500 ierr = MatDenseGetArrayRead(B,&barray);CHKERRQ(ierr); 501 ierr = MatDenseGetArray(X,&xarray);CHKERRQ(ierr); 502 503 if (barray == xarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"B and X cannot share the same memory location"); 504 if (!mat_mkl_pardiso->schur) mat_mkl_pardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT; 505 else mat_mkl_pardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION; 506 mat_mkl_pardiso->iparm[6-1] = 0; 507 508 MKL_PARDISO (mat_mkl_pardiso->pt, 509 &mat_mkl_pardiso->maxfct, 510 &mat_mkl_pardiso->mnum, 511 &mat_mkl_pardiso->mtype, 512 &mat_mkl_pardiso->phase, 513 &mat_mkl_pardiso->n, 514 mat_mkl_pardiso->a, 515 mat_mkl_pardiso->ia, 516 mat_mkl_pardiso->ja, 517 mat_mkl_pardiso->perm, 518 &mat_mkl_pardiso->nrhs, 519 mat_mkl_pardiso->iparm, 520 &mat_mkl_pardiso->msglvl, 521 (void*)barray, 522 (void*)xarray, 523 &mat_mkl_pardiso->err); 524 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err); 525 526 ierr = MatDenseRestoreArrayRead(B,&barray);CHKERRQ(ierr); 527 if (mat_mkl_pardiso->schur) { /* solve Schur complement and expand solution */ 528 PetscScalar *o_schur_work = NULL; 529 PetscInt shift = mat_mkl_pardiso->schur_size*mat_mkl_pardiso->nrhs,scale; 530 PetscInt mem = mat_mkl_pardiso->n*mat_mkl_pardiso->nrhs; 531 532 /* allocate extra memory if it is needed */ 533 scale = 1; 534 if (A->schur_status == MAT_FACTOR_SCHUR_INVERTED) scale = 2; 535 536 mem *= scale; 537 if (mem > mat_mkl_pardiso->schur_work_size) { 538 o_schur_work = mat_mkl_pardiso->schur_work; 539 ierr = PetscMalloc1(mem,&mat_mkl_pardiso->schur_work);CHKERRQ(ierr); 540 } 541 542 /* if inverted, uses BLAS *MM subroutines, otherwise LAPACK *TRS */ 543 if (A->schur_status != MAT_FACTOR_SCHUR_INVERTED) shift = 0; 544 545 /* solve Schur complement */ 546 if (!mat_mkl_pardiso->solve_interior) { 547 ierr = MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work,PETSC_TRUE);CHKERRQ(ierr); 548 ierr = MatMKLPardisoSolveSchur_Private(A,mat_mkl_pardiso->schur_work,mat_mkl_pardiso->schur_work+shift);CHKERRQ(ierr); 549 ierr = MatMKLPardisoScatterSchur_Private(mat_mkl_pardiso,xarray,mat_mkl_pardiso->schur_work+shift,PETSC_FALSE);CHKERRQ(ierr); 550 } else { /* if we are solving for the interior problem, any value in barray[schur,n] forward-substituted to xarray[schur,n] will be neglected */ 551 PetscInt i,n,m=0; 552 for (n=0;n<mat_mkl_pardiso->nrhs;n++) { 553 for (i=0;i<mat_mkl_pardiso->schur_size;i++) { 554 xarray[mat_mkl_pardiso->schur_idxs[i]+m] = 0.; 555 } 556 m += mat_mkl_pardiso->n; 557 } 558 } 559 560 /* expansion phase */ 561 mat_mkl_pardiso->iparm[6-1] = 1; 562 mat_mkl_pardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION; 563 MKL_PARDISO (mat_mkl_pardiso->pt, 564 &mat_mkl_pardiso->maxfct, 565 &mat_mkl_pardiso->mnum, 566 &mat_mkl_pardiso->mtype, 567 &mat_mkl_pardiso->phase, 568 &mat_mkl_pardiso->n, 569 mat_mkl_pardiso->a, 570 mat_mkl_pardiso->ia, 571 mat_mkl_pardiso->ja, 572 mat_mkl_pardiso->perm, 573 &mat_mkl_pardiso->nrhs, 574 mat_mkl_pardiso->iparm, 575 &mat_mkl_pardiso->msglvl, 576 (void*)xarray, 577 (void*)mat_mkl_pardiso->schur_work, /* according to the specs, the solution vector is always used */ 578 &mat_mkl_pardiso->err); 579 if (o_schur_work) { /* restore original schur_work (minimal size) */ 580 ierr = PetscFree(mat_mkl_pardiso->schur_work);CHKERRQ(ierr); 581 mat_mkl_pardiso->schur_work = o_schur_work; 582 } 583 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err); 584 mat_mkl_pardiso->iparm[6-1] = 0; 585 } 586 } 587 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 588 PetscFunctionReturn(0); 589 } 590 591 PetscErrorCode MatFactorNumeric_MKL_PARDISO(Mat F,Mat A,const MatFactorInfo *info) 592 { 593 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)(F)->data; 594 PetscErrorCode ierr; 595 596 PetscFunctionBegin; 597 mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN; 598 ierr = (*mat_mkl_pardiso->Convert)(A,mat_mkl_pardiso->needsym,MAT_REUSE_MATRIX,&mat_mkl_pardiso->freeaij,&mat_mkl_pardiso->nz,&mat_mkl_pardiso->ia,&mat_mkl_pardiso->ja,(PetscScalar**)&mat_mkl_pardiso->a);CHKERRQ(ierr); 599 600 mat_mkl_pardiso->phase = JOB_NUMERICAL_FACTORIZATION; 601 MKL_PARDISO (mat_mkl_pardiso->pt, 602 &mat_mkl_pardiso->maxfct, 603 &mat_mkl_pardiso->mnum, 604 &mat_mkl_pardiso->mtype, 605 &mat_mkl_pardiso->phase, 606 &mat_mkl_pardiso->n, 607 mat_mkl_pardiso->a, 608 mat_mkl_pardiso->ia, 609 mat_mkl_pardiso->ja, 610 mat_mkl_pardiso->perm, 611 &mat_mkl_pardiso->nrhs, 612 mat_mkl_pardiso->iparm, 613 &mat_mkl_pardiso->msglvl, 614 NULL, 615 (void*)mat_mkl_pardiso->schur, 616 &mat_mkl_pardiso->err); 617 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err); 618 619 if (F->schur) { /* schur output from pardiso is in row major format */ 620 ierr = MatFactorRestoreSchurComplement(F,NULL,MAT_FACTOR_SCHUR_UNFACTORED);CHKERRQ(ierr); 621 ierr = MatTranspose(F->schur,MAT_INPLACE_MATRIX,&F->schur);CHKERRQ(ierr); 622 } 623 mat_mkl_pardiso->matstruc = SAME_NONZERO_PATTERN; 624 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 625 PetscFunctionReturn(0); 626 } 627 628 PetscErrorCode PetscSetMKL_PARDISOFromOptions(Mat F, Mat A) 629 { 630 Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data; 631 PetscErrorCode ierr; 632 PetscInt icntl,threads=1; 633 PetscBool flg; 634 635 PetscFunctionBegin; 636 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MKL_PARDISO Options","Mat");CHKERRQ(ierr); 637 638 ierr = PetscOptionsInt("-mat_mkl_pardiso_65","Number of threads to use within PARDISO","None",threads,&threads,&flg);CHKERRQ(ierr); 639 if (flg) PetscSetMKL_PARDISOThreads((int)threads); 640 641 ierr = PetscOptionsInt("-mat_mkl_pardiso_66","Maximum number of factors with identical sparsity structure that must be kept in memory at the same time","None",mat_mkl_pardiso->maxfct,&icntl,&flg);CHKERRQ(ierr); 642 if (flg) mat_mkl_pardiso->maxfct = icntl; 643 644 ierr = PetscOptionsInt("-mat_mkl_pardiso_67","Indicates the actual matrix for the solution phase","None",mat_mkl_pardiso->mnum,&icntl,&flg);CHKERRQ(ierr); 645 if (flg) mat_mkl_pardiso->mnum = icntl; 646 647 ierr = PetscOptionsInt("-mat_mkl_pardiso_68","Message level information","None",mat_mkl_pardiso->msglvl,&icntl,&flg);CHKERRQ(ierr); 648 if (flg) mat_mkl_pardiso->msglvl = icntl; 649 650 ierr = PetscOptionsInt("-mat_mkl_pardiso_69","Defines the matrix type","None",mat_mkl_pardiso->mtype,&icntl,&flg);CHKERRQ(ierr); 651 if(flg){ 652 void *pt[IPARM_SIZE]; 653 mat_mkl_pardiso->mtype = icntl; 654 MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm); 655 #if defined(PETSC_USE_REAL_SINGLE) 656 mat_mkl_pardiso->iparm[27] = 1; 657 #else 658 mat_mkl_pardiso->iparm[27] = 0; 659 #endif 660 mat_mkl_pardiso->iparm[34] = 1; /* use 0-based indexing */ 661 } 662 ierr = PetscOptionsInt("-mat_mkl_pardiso_1","Use default values","None",mat_mkl_pardiso->iparm[0],&icntl,&flg);CHKERRQ(ierr); 663 664 if (flg && icntl != 0) { 665 ierr = PetscOptionsInt("-mat_mkl_pardiso_2","Fill-in reducing ordering for the input matrix","None",mat_mkl_pardiso->iparm[1],&icntl,&flg);CHKERRQ(ierr); 666 if (flg) mat_mkl_pardiso->iparm[1] = icntl; 667 668 ierr = PetscOptionsInt("-mat_mkl_pardiso_4","Preconditioned CGS/CG","None",mat_mkl_pardiso->iparm[3],&icntl,&flg);CHKERRQ(ierr); 669 if (flg) mat_mkl_pardiso->iparm[3] = icntl; 670 671 ierr = PetscOptionsInt("-mat_mkl_pardiso_5","User permutation","None",mat_mkl_pardiso->iparm[4],&icntl,&flg);CHKERRQ(ierr); 672 if (flg) mat_mkl_pardiso->iparm[4] = icntl; 673 674 ierr = PetscOptionsInt("-mat_mkl_pardiso_6","Write solution on x","None",mat_mkl_pardiso->iparm[5],&icntl,&flg);CHKERRQ(ierr); 675 if (flg) mat_mkl_pardiso->iparm[5] = icntl; 676 677 ierr = PetscOptionsInt("-mat_mkl_pardiso_8","Iterative refinement step","None",mat_mkl_pardiso->iparm[7],&icntl,&flg);CHKERRQ(ierr); 678 if (flg) mat_mkl_pardiso->iparm[7] = icntl; 679 680 ierr = PetscOptionsInt("-mat_mkl_pardiso_10","Pivoting perturbation","None",mat_mkl_pardiso->iparm[9],&icntl,&flg);CHKERRQ(ierr); 681 if (flg) mat_mkl_pardiso->iparm[9] = icntl; 682 683 ierr = PetscOptionsInt("-mat_mkl_pardiso_11","Scaling vectors","None",mat_mkl_pardiso->iparm[10],&icntl,&flg);CHKERRQ(ierr); 684 if (flg) mat_mkl_pardiso->iparm[10] = icntl; 685 686 ierr = PetscOptionsInt("-mat_mkl_pardiso_12","Solve with transposed or conjugate transposed matrix A","None",mat_mkl_pardiso->iparm[11],&icntl,&flg);CHKERRQ(ierr); 687 if (flg) mat_mkl_pardiso->iparm[11] = icntl; 688 689 ierr = PetscOptionsInt("-mat_mkl_pardiso_13","Improved accuracy using (non-) symmetric weighted matching","None",mat_mkl_pardiso->iparm[12],&icntl,&flg);CHKERRQ(ierr); 690 if (flg) mat_mkl_pardiso->iparm[12] = icntl; 691 692 ierr = PetscOptionsInt("-mat_mkl_pardiso_18","Numbers of non-zero elements","None",mat_mkl_pardiso->iparm[17],&icntl,&flg);CHKERRQ(ierr); 693 if (flg) mat_mkl_pardiso->iparm[17] = icntl; 694 695 ierr = PetscOptionsInt("-mat_mkl_pardiso_19","Report number of floating point operations","None",mat_mkl_pardiso->iparm[18],&icntl,&flg);CHKERRQ(ierr); 696 if (flg) mat_mkl_pardiso->iparm[18] = icntl; 697 698 ierr = PetscOptionsInt("-mat_mkl_pardiso_21","Pivoting for symmetric indefinite matrices","None",mat_mkl_pardiso->iparm[20],&icntl,&flg);CHKERRQ(ierr); 699 if (flg) mat_mkl_pardiso->iparm[20] = icntl; 700 701 ierr = PetscOptionsInt("-mat_mkl_pardiso_24","Parallel factorization control","None",mat_mkl_pardiso->iparm[23],&icntl,&flg);CHKERRQ(ierr); 702 if (flg) mat_mkl_pardiso->iparm[23] = icntl; 703 704 ierr = PetscOptionsInt("-mat_mkl_pardiso_25","Parallel forward/backward solve control","None",mat_mkl_pardiso->iparm[24],&icntl,&flg);CHKERRQ(ierr); 705 if (flg) mat_mkl_pardiso->iparm[24] = icntl; 706 707 ierr = PetscOptionsInt("-mat_mkl_pardiso_27","Matrix checker","None",mat_mkl_pardiso->iparm[26],&icntl,&flg);CHKERRQ(ierr); 708 if (flg) mat_mkl_pardiso->iparm[26] = icntl; 709 710 ierr = PetscOptionsInt("-mat_mkl_pardiso_31","Partial solve and computing selected components of the solution vectors","None",mat_mkl_pardiso->iparm[30],&icntl,&flg);CHKERRQ(ierr); 711 if (flg) mat_mkl_pardiso->iparm[30] = icntl; 712 713 ierr = PetscOptionsInt("-mat_mkl_pardiso_34","Optimal number of threads for conditional numerical reproducibility (CNR) mode","None",mat_mkl_pardiso->iparm[33],&icntl,&flg);CHKERRQ(ierr); 714 if (flg) mat_mkl_pardiso->iparm[33] = icntl; 715 716 ierr = PetscOptionsInt("-mat_mkl_pardiso_60","Intel MKL_PARDISO mode","None",mat_mkl_pardiso->iparm[59],&icntl,&flg);CHKERRQ(ierr); 717 if (flg) mat_mkl_pardiso->iparm[59] = icntl; 718 } 719 PetscOptionsEnd(); 720 PetscFunctionReturn(0); 721 } 722 723 PetscErrorCode MatFactorMKL_PARDISOInitialize_Private(Mat A, MatFactorType ftype, Mat_MKL_PARDISO *mat_mkl_pardiso) 724 { 725 PetscErrorCode ierr; 726 PetscInt i; 727 728 PetscFunctionBegin; 729 for ( i = 0; i < IPARM_SIZE; i++ ){ 730 mat_mkl_pardiso->iparm[i] = 0; 731 } 732 for ( i = 0; i < IPARM_SIZE; i++ ){ 733 mat_mkl_pardiso->pt[i] = 0; 734 } 735 /* Default options for both sym and unsym */ 736 mat_mkl_pardiso->iparm[ 0] = 1; /* Solver default parameters overriden with provided by iparm */ 737 mat_mkl_pardiso->iparm[ 1] = 2; /* Metis reordering */ 738 mat_mkl_pardiso->iparm[ 5] = 0; /* Write solution into x */ 739 mat_mkl_pardiso->iparm[ 7] = 0; /* Max number of iterative refinement steps */ 740 mat_mkl_pardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */ 741 mat_mkl_pardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */ 742 #if 0 743 mat_mkl_pardiso->iparm[23] = 1; /* Parallel factorization control*/ 744 #endif 745 mat_mkl_pardiso->iparm[34] = 1; /* Cluster Sparse Solver use C-style indexing for ia and ja arrays */ 746 mat_mkl_pardiso->iparm[39] = 0; /* Input: matrix/rhs/solution stored on master */ 747 748 mat_mkl_pardiso->CleanUp = PETSC_FALSE; 749 mat_mkl_pardiso->maxfct = 1; /* Maximum number of numerical factorizations. */ 750 mat_mkl_pardiso->mnum = 1; /* Which factorization to use. */ 751 mat_mkl_pardiso->msglvl = 0; /* 0: do not print 1: Print statistical information in file */ 752 mat_mkl_pardiso->phase = -1; 753 mat_mkl_pardiso->err = 0; 754 755 mat_mkl_pardiso->n = A->rmap->N; 756 mat_mkl_pardiso->nrhs = 1; 757 mat_mkl_pardiso->err = 0; 758 mat_mkl_pardiso->phase = -1; 759 760 if(ftype == MAT_FACTOR_LU){ 761 mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */ 762 mat_mkl_pardiso->iparm[10] = 1; /* Use nonsymmetric permutation and scaling MPS */ 763 mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */ 764 765 } else { 766 mat_mkl_pardiso->iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */ 767 mat_mkl_pardiso->iparm[10] = 0; /* Use nonsymmetric permutation and scaling MPS */ 768 mat_mkl_pardiso->iparm[12] = 1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */ 769 /* mat_mkl_pardiso->iparm[20] = 1; */ /* Apply 1x1 and 2x2 Bunch-Kaufman pivoting during the factorization process */ 770 #if defined(PETSC_USE_DEBUG) 771 mat_mkl_pardiso->iparm[26] = 1; /* Matrix checker */ 772 #endif 773 } 774 ierr = PetscMalloc1(A->rmap->N*sizeof(INT_TYPE), &mat_mkl_pardiso->perm);CHKERRQ(ierr); 775 for(i = 0; i < A->rmap->N; i++){ 776 mat_mkl_pardiso->perm[i] = 0; 777 } 778 mat_mkl_pardiso->schur_size = 0; 779 PetscFunctionReturn(0); 780 } 781 782 PetscErrorCode MatFactorSymbolic_AIJMKL_PARDISO_Private(Mat F,Mat A,const MatFactorInfo *info) 783 { 784 Mat_MKL_PARDISO *mat_mkl_pardiso = (Mat_MKL_PARDISO*)F->data; 785 PetscErrorCode ierr; 786 787 PetscFunctionBegin; 788 mat_mkl_pardiso->matstruc = DIFFERENT_NONZERO_PATTERN; 789 ierr = PetscSetMKL_PARDISOFromOptions(F,A);CHKERRQ(ierr); 790 791 /* throw away any previously computed structure */ 792 if (mat_mkl_pardiso->freeaij) { 793 ierr = PetscFree2(mat_mkl_pardiso->ia,mat_mkl_pardiso->ja);CHKERRQ(ierr); 794 ierr = PetscFree(mat_mkl_pardiso->a);CHKERRQ(ierr); 795 } 796 ierr = (*mat_mkl_pardiso->Convert)(A,mat_mkl_pardiso->needsym,MAT_INITIAL_MATRIX,&mat_mkl_pardiso->freeaij,&mat_mkl_pardiso->nz,&mat_mkl_pardiso->ia,&mat_mkl_pardiso->ja,(PetscScalar**)&mat_mkl_pardiso->a);CHKERRQ(ierr); 797 mat_mkl_pardiso->n = A->rmap->N; 798 799 mat_mkl_pardiso->phase = JOB_ANALYSIS; 800 801 MKL_PARDISO (mat_mkl_pardiso->pt, 802 &mat_mkl_pardiso->maxfct, 803 &mat_mkl_pardiso->mnum, 804 &mat_mkl_pardiso->mtype, 805 &mat_mkl_pardiso->phase, 806 &mat_mkl_pardiso->n, 807 mat_mkl_pardiso->a, 808 mat_mkl_pardiso->ia, 809 mat_mkl_pardiso->ja, 810 mat_mkl_pardiso->perm, 811 &mat_mkl_pardiso->nrhs, 812 mat_mkl_pardiso->iparm, 813 &mat_mkl_pardiso->msglvl, 814 NULL, 815 NULL, 816 &mat_mkl_pardiso->err); 817 if (mat_mkl_pardiso->err < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error reported by MKL_PARDISO: err=%d. Please check manual",mat_mkl_pardiso->err); 818 819 mat_mkl_pardiso->CleanUp = PETSC_TRUE; 820 821 if (F->factortype == MAT_FACTOR_LU) F->ops->lufactornumeric = MatFactorNumeric_MKL_PARDISO; 822 else F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_PARDISO; 823 824 F->ops->solve = MatSolve_MKL_PARDISO; 825 F->ops->solvetranspose = MatSolveTranspose_MKL_PARDISO; 826 F->ops->matsolve = MatMatSolve_MKL_PARDISO; 827 PetscFunctionReturn(0); 828 } 829 830 PetscErrorCode MatLUFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 831 { 832 PetscErrorCode ierr; 833 834 PetscFunctionBegin; 835 ierr = MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);CHKERRQ(ierr); 836 PetscFunctionReturn(0); 837 } 838 839 #if !defined(PETSC_USE_COMPLEX) 840 PetscErrorCode MatGetInertia_MKL_PARDISO(Mat F,PetscInt *nneg,PetscInt *nzero,PetscInt *npos) 841 { 842 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)F->data; 843 844 PetscFunctionBegin; 845 if (nneg) *nneg = mat_mkl_pardiso->iparm[22]; 846 if (npos) *npos = mat_mkl_pardiso->iparm[21]; 847 if (nzero) *nzero = F->rmap->N -(mat_mkl_pardiso->iparm[22] + mat_mkl_pardiso->iparm[21]); 848 PetscFunctionReturn(0); 849 } 850 #endif 851 852 PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_PARDISO(Mat F,Mat A,IS r,const MatFactorInfo *info) 853 { 854 PetscErrorCode ierr; 855 856 PetscFunctionBegin; 857 ierr = MatFactorSymbolic_AIJMKL_PARDISO_Private(F, A, info);CHKERRQ(ierr); 858 #if defined(PETSC_USE_COMPLEX) 859 F->ops->getinertia = NULL; 860 #else 861 F->ops->getinertia = MatGetInertia_MKL_PARDISO; 862 #endif 863 PetscFunctionReturn(0); 864 } 865 866 PetscErrorCode MatView_MKL_PARDISO(Mat A, PetscViewer viewer) 867 { 868 PetscErrorCode ierr; 869 PetscBool iascii; 870 PetscViewerFormat format; 871 Mat_MKL_PARDISO *mat_mkl_pardiso=(Mat_MKL_PARDISO*)A->data; 872 PetscInt i; 873 874 PetscFunctionBegin; 875 if (A->ops->solve != MatSolve_MKL_PARDISO) PetscFunctionReturn(0); 876 877 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 878 if (iascii) { 879 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 880 if (format == PETSC_VIEWER_ASCII_INFO) { 881 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO run parameters:\n");CHKERRQ(ierr); 882 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO phase: %d \n",mat_mkl_pardiso->phase);CHKERRQ(ierr); 883 for(i = 1; i <= 64; i++){ 884 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO iparm[%d]: %d \n",i, mat_mkl_pardiso->iparm[i - 1]);CHKERRQ(ierr); 885 } 886 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO maxfct: %d \n", mat_mkl_pardiso->maxfct);CHKERRQ(ierr); 887 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mnum: %d \n", mat_mkl_pardiso->mnum);CHKERRQ(ierr); 888 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO mtype: %d \n", mat_mkl_pardiso->mtype);CHKERRQ(ierr); 889 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO n: %d \n", mat_mkl_pardiso->n);CHKERRQ(ierr); 890 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO nrhs: %d \n", mat_mkl_pardiso->nrhs);CHKERRQ(ierr); 891 ierr = PetscViewerASCIIPrintf(viewer,"MKL_PARDISO msglvl: %d \n", mat_mkl_pardiso->msglvl);CHKERRQ(ierr); 892 } 893 } 894 PetscFunctionReturn(0); 895 } 896 897 898 PetscErrorCode MatGetInfo_MKL_PARDISO(Mat A, MatInfoType flag, MatInfo *info) 899 { 900 Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)A->data; 901 902 PetscFunctionBegin; 903 info->block_size = 1.0; 904 info->nz_used = mat_mkl_pardiso->nz; 905 info->nz_allocated = mat_mkl_pardiso->nz; 906 info->nz_unneeded = 0.0; 907 info->assemblies = 0.0; 908 info->mallocs = 0.0; 909 info->memory = 0.0; 910 info->fill_ratio_given = 0; 911 info->fill_ratio_needed = 0; 912 info->factor_mallocs = 0; 913 PetscFunctionReturn(0); 914 } 915 916 PetscErrorCode MatMkl_PardisoSetCntl_MKL_PARDISO(Mat F,PetscInt icntl,PetscInt ival) 917 { 918 Mat_MKL_PARDISO *mat_mkl_pardiso =(Mat_MKL_PARDISO*)F->data; 919 920 PetscFunctionBegin; 921 if(icntl <= 64){ 922 mat_mkl_pardiso->iparm[icntl - 1] = ival; 923 } else { 924 if(icntl == 65) PetscSetMKL_PARDISOThreads(ival); 925 else if(icntl == 66) mat_mkl_pardiso->maxfct = ival; 926 else if(icntl == 67) mat_mkl_pardiso->mnum = ival; 927 else if(icntl == 68) mat_mkl_pardiso->msglvl = ival; 928 else if(icntl == 69){ 929 void *pt[IPARM_SIZE]; 930 mat_mkl_pardiso->mtype = ival; 931 MKL_PARDISO_INIT(pt, &mat_mkl_pardiso->mtype, mat_mkl_pardiso->iparm); 932 #if defined(PETSC_USE_REAL_SINGLE) 933 mat_mkl_pardiso->iparm[27] = 1; 934 #else 935 mat_mkl_pardiso->iparm[27] = 0; 936 #endif 937 mat_mkl_pardiso->iparm[34] = 1; 938 } else if(icntl==70) mat_mkl_pardiso->solve_interior = (PetscBool)!!ival; 939 } 940 PetscFunctionReturn(0); 941 } 942 943 /*@ 944 MatMkl_PardisoSetCntl - Set Mkl_Pardiso parameters 945 946 Logically Collective on Mat 947 948 Input Parameters: 949 + F - the factored matrix obtained by calling MatGetFactor() 950 . icntl - index of Mkl_Pardiso parameter 951 - ival - value of Mkl_Pardiso parameter 952 953 Options Database: 954 . -mat_mkl_pardiso_<icntl> <ival> 955 956 Level: beginner 957 958 References: 959 . Mkl_Pardiso Users' Guide 960 961 .seealso: MatGetFactor() 962 @*/ 963 PetscErrorCode MatMkl_PardisoSetCntl(Mat F,PetscInt icntl,PetscInt ival) 964 { 965 PetscErrorCode ierr; 966 967 PetscFunctionBegin; 968 ierr = PetscTryMethod(F,"MatMkl_PardisoSetCntl_C",(Mat,PetscInt,PetscInt),(F,icntl,ival));CHKERRQ(ierr); 969 PetscFunctionReturn(0); 970 } 971 972 /*MC 973 MATSOLVERMKL_PARDISO - A matrix type providing direct solvers (LU) for 974 sequential matrices via the external package MKL_PARDISO. 975 976 Works with MATSEQAIJ matrices 977 978 Use -pc_type lu -pc_factor_mat_solver_type mkl_pardiso to use this direct solver 979 980 Options Database Keys: 981 + -mat_mkl_pardiso_65 - Number of threads to use within MKL_PARDISO 982 . -mat_mkl_pardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time 983 . -mat_mkl_pardiso_67 - Indicates the actual matrix for the solution phase 984 . -mat_mkl_pardiso_68 - Message level information 985 . -mat_mkl_pardiso_69 - Defines the matrix type. IMPORTANT: When you set this flag, iparm parameters are going to be set to the default ones for the matrix type 986 . -mat_mkl_pardiso_1 - Use default values 987 . -mat_mkl_pardiso_2 - Fill-in reducing ordering for the input matrix 988 . -mat_mkl_pardiso_4 - Preconditioned CGS/CG 989 . -mat_mkl_pardiso_5 - User permutation 990 . -mat_mkl_pardiso_6 - Write solution on x 991 . -mat_mkl_pardiso_8 - Iterative refinement step 992 . -mat_mkl_pardiso_10 - Pivoting perturbation 993 . -mat_mkl_pardiso_11 - Scaling vectors 994 . -mat_mkl_pardiso_12 - Solve with transposed or conjugate transposed matrix A 995 . -mat_mkl_pardiso_13 - Improved accuracy using (non-) symmetric weighted matching 996 . -mat_mkl_pardiso_18 - Numbers of non-zero elements 997 . -mat_mkl_pardiso_19 - Report number of floating point operations 998 . -mat_mkl_pardiso_21 - Pivoting for symmetric indefinite matrices 999 . -mat_mkl_pardiso_24 - Parallel factorization control 1000 . -mat_mkl_pardiso_25 - Parallel forward/backward solve control 1001 . -mat_mkl_pardiso_27 - Matrix checker 1002 . -mat_mkl_pardiso_31 - Partial solve and computing selected components of the solution vectors 1003 . -mat_mkl_pardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode 1004 - -mat_mkl_pardiso_60 - Intel MKL_PARDISO mode 1005 1006 Level: beginner 1007 1008 For more information please check mkl_pardiso manual 1009 1010 .seealso: PCFactorSetMatSolverType(), MatSolverType 1011 1012 M*/ 1013 static PetscErrorCode MatFactorGetSolverType_mkl_pardiso(Mat A, MatSolverType *type) 1014 { 1015 PetscFunctionBegin; 1016 *type = MATSOLVERMKL_PARDISO; 1017 PetscFunctionReturn(0); 1018 } 1019 1020 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mkl_pardiso(Mat A,MatFactorType ftype,Mat *F) 1021 { 1022 Mat B; 1023 PetscErrorCode ierr; 1024 Mat_MKL_PARDISO *mat_mkl_pardiso; 1025 PetscBool isSeqAIJ,isSeqBAIJ,isSeqSBAIJ; 1026 1027 PetscFunctionBegin; 1028 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isSeqAIJ);CHKERRQ(ierr); 1029 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&isSeqBAIJ);CHKERRQ(ierr); 1030 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&isSeqSBAIJ);CHKERRQ(ierr); 1031 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 1032 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1033 ierr = PetscStrallocpy("mkl_pardiso",&((PetscObject)B)->type_name);CHKERRQ(ierr); 1034 ierr = MatSetUp(B);CHKERRQ(ierr); 1035 1036 ierr = PetscNewLog(B,&mat_mkl_pardiso);CHKERRQ(ierr); 1037 B->data = mat_mkl_pardiso; 1038 1039 ierr = MatFactorMKL_PARDISOInitialize_Private(A, ftype, mat_mkl_pardiso);CHKERRQ(ierr); 1040 if (ftype == MAT_FACTOR_LU) { 1041 B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_PARDISO; 1042 B->factortype = MAT_FACTOR_LU; 1043 mat_mkl_pardiso->needsym = PETSC_FALSE; 1044 if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij; 1045 else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij; 1046 else if (isSeqSBAIJ) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU factor with SEQSBAIJ format! Use MAT_FACTOR_CHOLESKY instead"); 1047 else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO LU with %s format",((PetscObject)A)->type_name); 1048 #if defined(PETSC_USE_COMPLEX) 1049 mat_mkl_pardiso->mtype = 13; 1050 #else 1051 mat_mkl_pardiso->mtype = 11; 1052 #endif 1053 } else { 1054 #if defined(PETSC_USE_COMPLEX) 1055 SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead",((PetscObject)A)->type_name); 1056 #endif 1057 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_PARDISO; 1058 B->factortype = MAT_FACTOR_CHOLESKY; 1059 if (isSeqAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqaij; 1060 else if (isSeqBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqbaij; 1061 else if (isSeqSBAIJ) mat_mkl_pardiso->Convert = MatMKLPardiso_Convert_seqsbaij; 1062 else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for PARDISO CHOLESKY with %s format",((PetscObject)A)->type_name); 1063 1064 mat_mkl_pardiso->needsym = PETSC_TRUE; 1065 if (A->spd_set && A->spd) mat_mkl_pardiso->mtype = 2; 1066 else mat_mkl_pardiso->mtype = -2; 1067 } 1068 B->ops->destroy = MatDestroy_MKL_PARDISO; 1069 B->ops->view = MatView_MKL_PARDISO; 1070 B->factortype = ftype; 1071 B->ops->getinfo = MatGetInfo_MKL_PARDISO; 1072 B->assembled = PETSC_TRUE; 1073 1074 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 1075 ierr = PetscStrallocpy(MATSOLVERMKL_PARDISO,&B->solvertype);CHKERRQ(ierr); 1076 1077 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_mkl_pardiso);CHKERRQ(ierr); 1078 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorSetSchurIS_C",MatFactorSetSchurIS_MKL_PARDISO);CHKERRQ(ierr); 1079 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMkl_PardisoSetCntl_C",MatMkl_PardisoSetCntl_MKL_PARDISO);CHKERRQ(ierr); 1080 1081 *F = B; 1082 PetscFunctionReturn(0); 1083 } 1084 1085 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_MKL_Pardiso(void) 1086 { 1087 PetscErrorCode ierr; 1088 1089 PetscFunctionBegin; 1090 ierr = MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_mkl_pardiso);CHKERRQ(ierr); 1091 ierr = MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);CHKERRQ(ierr); 1092 ierr = MatSolverTypeRegister(MATSOLVERMKL_PARDISO,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_mkl_pardiso);CHKERRQ(ierr); 1093 PetscFunctionReturn(0); 1094 } 1095 1096