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