#include <petscsys.h>
#include <../src/mat/impls/aij/mpi/mpiaij.h> /*I  "petscmat.h"  I*/
#include <../src/mat/impls/sbaij/mpi/mpisbaij.h>

#if defined(PETSC_HAVE_MKL_INTEL_ILP64)
  #define MKL_ILP64
#endif
#include <mkl.h>
#include <mkl_cluster_sparse_solver.h>

/*
 *  Possible mkl_cpardiso phases that controls the execution of the solver.
 *  For more information check mkl_cpardiso manual.
 */
#define JOB_ANALYSIS                                                    11
#define JOB_ANALYSIS_NUMERICAL_FACTORIZATION                            12
#define JOB_ANALYSIS_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT 13
#define JOB_NUMERICAL_FACTORIZATION                                     22
#define JOB_NUMERICAL_FACTORIZATION_SOLVE_ITERATIVE_REFINEMENT          23
#define JOB_SOLVE_ITERATIVE_REFINEMENT                                  33
#define JOB_SOLVE_FORWARD_SUBSTITUTION                                  331
#define JOB_SOLVE_DIAGONAL_SUBSTITUTION                                 332
#define JOB_SOLVE_BACKWARD_SUBSTITUTION                                 333
#define JOB_RELEASE_OF_LU_MEMORY                                        0
#define JOB_RELEASE_OF_ALL_MEMORY                                       -1

#define IPARM_SIZE 64
#define INT_TYPE   MKL_INT

static const char *Err_MSG_CPardiso(int errNo)
{
  switch (errNo) {
  case -1:
    return "input inconsistent";
    break;
  case -2:
    return "not enough memory";
    break;
  case -3:
    return "reordering problem";
    break;
  case -4:
    return "zero pivot, numerical factorization or iterative refinement problem";
    break;
  case -5:
    return "unclassified (internal) error";
    break;
  case -6:
    return "preordering failed (matrix types 11, 13 only)";
    break;
  case -7:
    return "diagonal matrix problem";
    break;
  case -8:
    return "32-bit integer overflow problem";
    break;
  case -9:
    return "not enough memory for OOC";
    break;
  case -10:
    return "problems with opening OOC temporary files";
    break;
  case -11:
    return "read/write problems with the OOC data file";
    break;
  default:
    return "unknown error";
  }
}

#define PetscCallCluster(f) PetscStackCallExternalVoid("cluster_sparse_solver", f);

/*
 *  Internal data structure.
 *  For more information check mkl_cpardiso manual.
 */

typedef struct {
  /* Configuration vector */
  INT_TYPE iparm[IPARM_SIZE];

  /*
   * Internal mkl_cpardiso memory location.
   * After the first call to mkl_cpardiso do not modify pt, as that could cause a serious memory leak.
   */
  void *pt[IPARM_SIZE];

  MPI_Fint comm_mkl_cpardiso;

  /* Basic mkl_cpardiso info*/
  INT_TYPE phase, maxfct, mnum, mtype, n, nrhs, msglvl, err;

  /* Matrix values and matrix nonzero structure */
  PetscScalar *a;

  INT_TYPE *ia, *ja;

  /* Number of non-zero elements */
  INT_TYPE nz;

  /* Row permutaton vector*/
  INT_TYPE *perm;

  /* Define is matrix preserve sparse structure. */
  MatStructure matstruc;

  PetscErrorCode (*ConvertToTriples)(Mat, MatReuse, PetscInt *, PetscInt **, PetscInt **, PetscScalar **);

  /* True if mkl_cpardiso function have been used. */
  PetscBool CleanUp;
} Mat_MKL_CPARDISO;

/*
 * Copy the elements of matrix A.
 * Input:
 *   - Mat A: MATSEQAIJ matrix
 *   - int shift: matrix index.
 *     - 0 for c representation
 *     - 1 for fortran representation
 *   - MatReuse reuse:
 *     - MAT_INITIAL_MATRIX: Create a new aij representation
 *     - MAT_REUSE_MATRIX: Reuse all aij representation and just change values
 * Output:
 *   - int *nnz: Number of nonzero-elements.
 *   - int **r pointer to i index
 *   - int **c pointer to j elements
 *   - MATRIXTYPE **v: Non-zero elements
 */
static PetscErrorCode MatCopy_seqaij_seqaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
{
  Mat_SeqAIJ *aa = (Mat_SeqAIJ *)A->data;

  PetscFunctionBegin;
  *v = aa->a;
  if (reuse == MAT_INITIAL_MATRIX) {
    *r   = (INT_TYPE *)aa->i;
    *c   = (INT_TYPE *)aa->j;
    *nnz = aa->nz;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
{
  const PetscInt    *ai, *aj, *bi, *bj, *garray, m = A->rmap->n, *ajj, *bjj;
  PetscInt           rstart, nz, i, j, countA, countB;
  PetscInt          *row, *col;
  const PetscScalar *av, *bv;
  PetscScalar       *val;
  Mat_MPIAIJ        *mat = (Mat_MPIAIJ *)A->data;
  Mat_SeqAIJ        *aa  = (Mat_SeqAIJ *)mat->A->data;
  Mat_SeqAIJ        *bb  = (Mat_SeqAIJ *)mat->B->data;
  PetscInt           colA_start, jB, jcol;

  PetscFunctionBegin;
  ai     = aa->i;
  aj     = aa->j;
  bi     = bb->i;
  bj     = bb->j;
  rstart = A->rmap->rstart;
  av     = aa->a;
  bv     = bb->a;

  garray = mat->garray;

  if (reuse == MAT_INITIAL_MATRIX) {
    nz   = aa->nz + bb->nz;
    *nnz = nz;
    PetscCall(PetscMalloc3(m + 1, &row, nz, &col, nz, &val));
    *r = row;
    *c = col;
    *v = val;
  } else {
    row = *r;
    col = *c;
    val = *v;
  }

  nz = 0;
  for (i = 0; i < m; i++) {
    row[i] = nz;
    countA = ai[i + 1] - ai[i];
    countB = bi[i + 1] - bi[i];
    ajj    = aj + ai[i]; /* ptr to the beginning of this row */
    bjj    = bj + bi[i];

    /* B part, smaller col index */
    colA_start = rstart + ajj[0]; /* the smallest global col index of A */
    jB         = 0;
    for (j = 0; j < countB; j++) {
      jcol = garray[bjj[j]];
      if (jcol > colA_start) break;
      col[nz]   = jcol;
      val[nz++] = *bv++;
    }
    jB = j;

    /* A part */
    for (j = 0; j < countA; j++) {
      col[nz]   = rstart + ajj[j];
      val[nz++] = *av++;
    }

    /* B part, larger col index */
    for (j = jB; j < countB; j++) {
      col[nz]   = garray[bjj[j]];
      val[nz++] = *bv++;
    }
  }
  row[m] = nz;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
{
  const PetscInt    *ai, *aj, *bi, *bj, *garray, bs = A->rmap->bs, bs2 = bs * bs, m = A->rmap->n / bs, *ajj, *bjj;
  PetscInt           rstart, nz, i, j, countA, countB;
  PetscInt          *row, *col;
  const PetscScalar *av, *bv;
  PetscScalar       *val;
  Mat_MPIBAIJ       *mat = (Mat_MPIBAIJ *)A->data;
  Mat_SeqBAIJ       *aa  = (Mat_SeqBAIJ *)mat->A->data;
  Mat_SeqBAIJ       *bb  = (Mat_SeqBAIJ *)mat->B->data;
  PetscInt           colA_start, jB, jcol;

  PetscFunctionBegin;
  ai     = aa->i;
  aj     = aa->j;
  bi     = bb->i;
  bj     = bb->j;
  rstart = A->rmap->rstart / bs;
  av     = aa->a;
  bv     = bb->a;

  garray = mat->garray;

  if (reuse == MAT_INITIAL_MATRIX) {
    nz   = aa->nz + bb->nz;
    *nnz = nz;
    PetscCall(PetscMalloc3(m + 1, &row, nz, &col, nz * bs2, &val));
    *r = row;
    *c = col;
    *v = val;
  } else {
    row = *r;
    col = *c;
    val = *v;
  }

  nz = 0;
  for (i = 0; i < m; i++) {
    row[i] = nz + 1;
    countA = ai[i + 1] - ai[i];
    countB = bi[i + 1] - bi[i];
    ajj    = aj + ai[i]; /* ptr to the beginning of this row */
    bjj    = bj + bi[i];

    /* B part, smaller col index */
    colA_start = rstart + (countA > 0 ? ajj[0] : 0); /* the smallest global col index of A */
    jB         = 0;
    for (j = 0; j < countB; j++) {
      jcol = garray[bjj[j]];
      if (jcol > colA_start) break;
      col[nz++] = jcol + 1;
    }
    jB = j;
    PetscCall(PetscArraycpy(val, bv, jB * bs2));
    val += jB * bs2;
    bv += jB * bs2;

    /* A part */
    for (j = 0; j < countA; j++) col[nz++] = rstart + ajj[j] + 1;
    PetscCall(PetscArraycpy(val, av, countA * bs2));
    val += countA * bs2;
    av += countA * bs2;

    /* B part, larger col index */
    for (j = jB; j < countB; j++) col[nz++] = garray[bjj[j]] + 1;
    PetscCall(PetscArraycpy(val, bv, (countB - jB) * bs2));
    val += (countB - jB) * bs2;
    bv += (countB - jB) * bs2;
  }
  row[m] = nz + 1;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO(Mat A, MatReuse reuse, PetscInt *nnz, PetscInt **r, PetscInt **c, PetscScalar **v)
{
  const PetscInt    *ai, *aj, *bi, *bj, *garray, bs = A->rmap->bs, bs2 = bs * bs, m = A->rmap->n / bs, *ajj, *bjj;
  PetscInt           rstart, nz, i, j, countA, countB;
  PetscInt          *row, *col;
  const PetscScalar *av, *bv;
  PetscScalar       *val;
  Mat_MPISBAIJ      *mat = (Mat_MPISBAIJ *)A->data;
  Mat_SeqSBAIJ      *aa  = (Mat_SeqSBAIJ *)mat->A->data;
  Mat_SeqBAIJ       *bb  = (Mat_SeqBAIJ *)mat->B->data;

  PetscFunctionBegin;
  ai     = aa->i;
  aj     = aa->j;
  bi     = bb->i;
  bj     = bb->j;
  rstart = A->rmap->rstart / bs;
  av     = aa->a;
  bv     = bb->a;

  garray = mat->garray;

  if (reuse == MAT_INITIAL_MATRIX) {
    nz   = aa->nz + bb->nz;
    *nnz = nz;
    PetscCall(PetscMalloc3(m + 1, &row, nz, &col, nz * bs2, &val));
    *r = row;
    *c = col;
    *v = val;
  } else {
    row = *r;
    col = *c;
    val = *v;
  }

  nz = 0;
  for (i = 0; i < m; i++) {
    row[i] = nz + 1;
    countA = ai[i + 1] - ai[i];
    countB = bi[i + 1] - bi[i];
    ajj    = aj + ai[i]; /* ptr to the beginning of this row */
    bjj    = bj + bi[i];

    /* A part */
    for (j = 0; j < countA; j++) col[nz++] = rstart + ajj[j] + 1;
    PetscCall(PetscArraycpy(val, av, countA * bs2));
    val += countA * bs2;
    av += countA * bs2;

    /* B part, larger col index */
    for (j = 0; j < countB; j++) col[nz++] = garray[bjj[j]] + 1;
    PetscCall(PetscArraycpy(val, bv, countB * bs2));
    val += countB * bs2;
    bv += countB * bs2;
  }
  row[m] = nz + 1;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
 * Free memory for Mat_MKL_CPARDISO structure and pointers to objects.
 */
static PetscErrorCode MatDestroy_MKL_CPARDISO(Mat A)
{
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
  MPI_Comm          comm;

  PetscFunctionBegin;
  /* Terminate instance, deallocate memories */
  if (mat_mkl_cpardiso->CleanUp) {
    mat_mkl_cpardiso->phase = JOB_RELEASE_OF_ALL_MEMORY;

    PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, NULL, NULL, NULL, mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs,
                                           mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
  }
  if (mat_mkl_cpardiso->ConvertToTriples != MatCopy_seqaij_seqaij_MKL_CPARDISO) PetscCall(PetscFree3(mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja, mat_mkl_cpardiso->a));
  comm = MPI_Comm_f2c(mat_mkl_cpardiso->comm_mkl_cpardiso);
  PetscCallMPI(MPI_Comm_free(&comm));
  PetscCall(PetscFree(A->data));

  /* clear composed functions */
  PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatMkl_CPardisoSetCntl_C", NULL));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
 * Computes Ax = b
 */
static PetscErrorCode MatSolve_MKL_CPARDISO(Mat A, Vec b, Vec x)
{
  Mat_MKL_CPARDISO  *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
  PetscScalar       *xarray;
  const PetscScalar *barray;

  PetscFunctionBegin;
  mat_mkl_cpardiso->nrhs = 1;
  PetscCall(VecGetArray(x, &xarray));
  PetscCall(VecGetArrayRead(b, &barray));

  /* solve phase */
  mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
  PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
                                         mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
  PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));

  PetscCall(VecRestoreArray(x, &xarray));
  PetscCall(VecRestoreArrayRead(b, &barray));
  mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatForwardSolve_MKL_CPARDISO(Mat A, Vec b, Vec x)
{
  Mat_MKL_CPARDISO  *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
  PetscScalar       *xarray;
  const PetscScalar *barray;

  PetscFunctionBegin;
  mat_mkl_cpardiso->nrhs = 1;
  PetscCall(VecGetArray(x, &xarray));
  PetscCall(VecGetArrayRead(b, &barray));

  /* solve phase */
  mat_mkl_cpardiso->phase = JOB_SOLVE_FORWARD_SUBSTITUTION;
  PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
                                         mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
  PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));

  PetscCall(VecRestoreArray(x, &xarray));
  PetscCall(VecRestoreArrayRead(b, &barray));
  mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatBackwardSolve_MKL_CPARDISO(Mat A, Vec b, Vec x)
{
  Mat_MKL_CPARDISO  *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
  PetscScalar       *xarray;
  const PetscScalar *barray;

  PetscFunctionBegin;
  mat_mkl_cpardiso->nrhs = 1;
  PetscCall(VecGetArray(x, &xarray));
  PetscCall(VecGetArrayRead(b, &barray));

  /* solve phase */
  mat_mkl_cpardiso->phase = JOB_SOLVE_BACKWARD_SUBSTITUTION;
  PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
                                         mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
  PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));

  PetscCall(VecRestoreArray(x, &xarray));
  PetscCall(VecRestoreArrayRead(b, &barray));
  mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatSolveTranspose_MKL_CPARDISO(Mat A, Vec b, Vec x)
{
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;

  PetscFunctionBegin;
  mat_mkl_cpardiso->iparm[12 - 1] = PetscDefined(USE_COMPLEX) ? 1 : 2;
  PetscCall(MatSolve_MKL_CPARDISO(A, b, x));
  mat_mkl_cpardiso->iparm[12 - 1] = 0;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatMatSolve_MKL_CPARDISO(Mat A, Mat B, Mat X)
{
  Mat_MKL_CPARDISO  *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
  PetscScalar       *xarray;
  const PetscScalar *barray;

  PetscFunctionBegin;
  PetscCall(MatGetSize(B, NULL, (PetscInt *)&mat_mkl_cpardiso->nrhs));

  if (mat_mkl_cpardiso->nrhs > 0) {
    PetscCall(MatDenseGetArrayRead(B, &barray));
    PetscCall(MatDenseGetArray(X, &xarray));

    PetscCheck(barray != xarray, PETSC_COMM_SELF, PETSC_ERR_SUP, "B and X cannot share the same memory location");

    /* solve phase */
    mat_mkl_cpardiso->phase = JOB_SOLVE_ITERATIVE_REFINEMENT;
    PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
                                           mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, (void *)barray, (void *)xarray, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
    PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));
    PetscCall(MatDenseRestoreArrayRead(B, &barray));
    PetscCall(MatDenseRestoreArray(X, &xarray));
  }
  mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
 * LU Decomposition
 */
static PetscErrorCode MatFactorNumeric_MKL_CPARDISO(Mat F, Mat A, const MatFactorInfo *info)
{
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;

  PetscFunctionBegin;
  mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
  PetscCall((*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_REUSE_MATRIX, &mat_mkl_cpardiso->nz, &mat_mkl_cpardiso->ia, &mat_mkl_cpardiso->ja, &mat_mkl_cpardiso->a));

  mat_mkl_cpardiso->phase = JOB_NUMERICAL_FACTORIZATION;
  PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
                                         mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, &mat_mkl_cpardiso->err));
  PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\". Please check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));

  mat_mkl_cpardiso->matstruc = SAME_NONZERO_PATTERN;
  mat_mkl_cpardiso->CleanUp  = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* Sets mkl_cpardiso options from the options database */
static PetscErrorCode MatSetFromOptions_MKL_CPARDISO(Mat F, Mat A)
{
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;
  PetscInt          icntl, threads;
  PetscBool         flg;

  PetscFunctionBegin;
  PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "MKL Cluster PARDISO Options", "Mat");
  PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_65", "Suggested number of threads to use within MKL Cluster PARDISO", "None", threads, &threads, &flg));
  if (flg) mkl_set_num_threads((int)threads);

  PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_66", "Maximum number of factors with identical sparsity structure that must be kept in memory at the same time", "None", mat_mkl_cpardiso->maxfct, &icntl, &flg));
  if (flg) mat_mkl_cpardiso->maxfct = icntl;

  PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_67", "Indicates the actual matrix for the solution phase", "None", mat_mkl_cpardiso->mnum, &icntl, &flg));
  if (flg) mat_mkl_cpardiso->mnum = icntl;

  PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_68", "Message level information", "None", mat_mkl_cpardiso->msglvl, &icntl, &flg));
  if (flg) mat_mkl_cpardiso->msglvl = icntl;

  PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_69", "Defines the matrix type", "None", mat_mkl_cpardiso->mtype, &icntl, &flg));
  if (flg) mat_mkl_cpardiso->mtype = icntl;
  PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_1", "Use default values", "None", mat_mkl_cpardiso->iparm[0], &icntl, &flg));

  if (flg && icntl != 0) {
    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_2", "Fill-in reducing ordering for the input matrix", "None", mat_mkl_cpardiso->iparm[1], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[1] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_4", "Preconditioned CGS/CG", "None", mat_mkl_cpardiso->iparm[3], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[3] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_5", "User permutation", "None", mat_mkl_cpardiso->iparm[4], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[4] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_6", "Write solution on x", "None", mat_mkl_cpardiso->iparm[5], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[5] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_8", "Iterative refinement step", "None", mat_mkl_cpardiso->iparm[7], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[7] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_10", "Pivoting perturbation", "None", mat_mkl_cpardiso->iparm[9], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[9] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_11", "Scaling vectors", "None", mat_mkl_cpardiso->iparm[10], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[10] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_12", "Solve with transposed or conjugate transposed matrix A", "None", mat_mkl_cpardiso->iparm[11], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[11] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_13", "Improved accuracy using (non-) symmetric weighted matching", "None", mat_mkl_cpardiso->iparm[12], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[12] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_18", "Numbers of non-zero elements", "None", mat_mkl_cpardiso->iparm[17], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[17] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_19", "Report number of floating point operations", "None", mat_mkl_cpardiso->iparm[18], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[18] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_21", "Pivoting for symmetric indefinite matrices", "None", mat_mkl_cpardiso->iparm[20], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[20] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_24", "Parallel factorization control", "None", mat_mkl_cpardiso->iparm[23], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[23] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_25", "Parallel forward/backward solve control", "None", mat_mkl_cpardiso->iparm[24], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[24] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_27", "Matrix checker", "None", mat_mkl_cpardiso->iparm[26], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[26] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_31", "Partial solve and computing selected components of the solution vectors", "None", mat_mkl_cpardiso->iparm[30], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[30] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_34", "Optimal number of threads for conditional numerical reproducibility (CNR) mode", "None", mat_mkl_cpardiso->iparm[33], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[33] = icntl;

    PetscCall(PetscOptionsInt("-mat_mkl_cpardiso_60", "Intel MKL Cluster PARDISO mode", "None", mat_mkl_cpardiso->iparm[59], &icntl, &flg));
    if (flg) mat_mkl_cpardiso->iparm[59] = icntl;
  }

  PetscOptionsEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PetscInitialize_MKL_CPARDISO(Mat A, Mat_MKL_CPARDISO *mat_mkl_cpardiso)
{
  PetscInt    bs;
  PetscBool   match;
  PetscMPIInt size;
  MPI_Comm    comm;

  PetscFunctionBegin;
  PetscCallMPI(MPI_Comm_dup(PetscObjectComm((PetscObject)A), &comm));
  PetscCallMPI(MPI_Comm_size(comm, &size));
  mat_mkl_cpardiso->comm_mkl_cpardiso = MPI_Comm_c2f(comm);

  mat_mkl_cpardiso->CleanUp = PETSC_FALSE;
  mat_mkl_cpardiso->maxfct  = 1;
  mat_mkl_cpardiso->mnum    = 1;
  mat_mkl_cpardiso->n       = A->rmap->N;
  if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
  mat_mkl_cpardiso->msglvl = 0;
  mat_mkl_cpardiso->nrhs   = 1;
  mat_mkl_cpardiso->err    = 0;
  mat_mkl_cpardiso->phase  = -1;
  mat_mkl_cpardiso->mtype  = PetscDefined(USE_COMPLEX) ? 13 : 11;

  mat_mkl_cpardiso->iparm[27] = PetscDefined(USE_REAL_SINGLE) ? 1 : 0;

  mat_mkl_cpardiso->iparm[0]  = 1;  /* Solver default parameters overridden with provided by iparm */
  mat_mkl_cpardiso->iparm[1]  = 2;  /* Use METIS for fill-in reordering */
  mat_mkl_cpardiso->iparm[5]  = 0;  /* Write solution into x */
  mat_mkl_cpardiso->iparm[7]  = 2;  /* Max number of iterative refinement steps */
  mat_mkl_cpardiso->iparm[9]  = 13; /* Perturb the pivot elements with 1E-13 */
  mat_mkl_cpardiso->iparm[10] = 1;  /* Use nonsymmetric permutation and scaling MPS */
  mat_mkl_cpardiso->iparm[12] = 1;  /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
  mat_mkl_cpardiso->iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
  mat_mkl_cpardiso->iparm[18] = -1; /* Output: Mflops for LU factorization */
  mat_mkl_cpardiso->iparm[26] = 1;  /* Check input data for correctness */

  mat_mkl_cpardiso->iparm[39] = 0;
  if (size > 1) {
    mat_mkl_cpardiso->iparm[39] = 2;
    mat_mkl_cpardiso->iparm[40] = A->rmap->rstart;
    mat_mkl_cpardiso->iparm[41] = A->rmap->rend - 1;
  }
  PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &match, MATMPIBAIJ, MATMPISBAIJ, ""));
  if (match) {
    PetscCall(MatGetBlockSize(A, &bs));
    mat_mkl_cpardiso->iparm[36] = bs;
    mat_mkl_cpardiso->iparm[40] /= bs;
    mat_mkl_cpardiso->iparm[41] /= bs;
    mat_mkl_cpardiso->iparm[40]++;
    mat_mkl_cpardiso->iparm[41]++;
    mat_mkl_cpardiso->iparm[34] = 0; /* Fortran style */
  } else {
    mat_mkl_cpardiso->iparm[34] = 1; /* C style */
  }

  mat_mkl_cpardiso->perm = 0;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
 * Symbolic decomposition. Mkl_Pardiso analysis phase.
 */
static PetscErrorCode MatLUFactorSymbolic_AIJMKL_CPARDISO(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
{
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;

  PetscFunctionBegin;
  mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;

  /* Set MKL_CPARDISO options from the options database */
  PetscCall(MatSetFromOptions_MKL_CPARDISO(F, A));
  PetscCall((*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_INITIAL_MATRIX, &mat_mkl_cpardiso->nz, &mat_mkl_cpardiso->ia, &mat_mkl_cpardiso->ja, &mat_mkl_cpardiso->a));

  mat_mkl_cpardiso->n = A->rmap->N;
  if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];

  /* analysis phase */
  mat_mkl_cpardiso->phase = JOB_ANALYSIS;

  PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
                                         mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
  PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\".Check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));

  mat_mkl_cpardiso->CleanUp = PETSC_TRUE;
  F->ops->lufactornumeric   = MatFactorNumeric_MKL_CPARDISO;
  F->ops->solve             = MatSolve_MKL_CPARDISO;
  F->ops->forwardsolve      = MatForwardSolve_MKL_CPARDISO;
  F->ops->backwardsolve     = MatBackwardSolve_MKL_CPARDISO;
  F->ops->solvetranspose    = MatSolveTranspose_MKL_CPARDISO;
  F->ops->matsolve          = MatMatSolve_MKL_CPARDISO;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatCholeskyFactorSymbolic_AIJMKL_CPARDISO(Mat F, Mat A, IS perm, const MatFactorInfo *info)
{
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;

  PetscFunctionBegin;
  mat_mkl_cpardiso->matstruc = DIFFERENT_NONZERO_PATTERN;

  /* Set MKL_CPARDISO options from the options database */
  PetscCall(MatSetFromOptions_MKL_CPARDISO(F, A));
  PetscCall((*mat_mkl_cpardiso->ConvertToTriples)(A, MAT_INITIAL_MATRIX, &mat_mkl_cpardiso->nz, &mat_mkl_cpardiso->ia, &mat_mkl_cpardiso->ja, &mat_mkl_cpardiso->a));

  mat_mkl_cpardiso->n = A->rmap->N;
  if (mat_mkl_cpardiso->iparm[36]) mat_mkl_cpardiso->n /= mat_mkl_cpardiso->iparm[36];
  PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for PARDISO CHOLESKY with complex scalars! Use MAT_FACTOR_LU instead");
  if (A->spd == PETSC_BOOL3_TRUE) mat_mkl_cpardiso->mtype = 2;
  else mat_mkl_cpardiso->mtype = -2;

  /* analysis phase */
  mat_mkl_cpardiso->phase = JOB_ANALYSIS;

  PetscCallCluster(cluster_sparse_solver(mat_mkl_cpardiso->pt, &mat_mkl_cpardiso->maxfct, &mat_mkl_cpardiso->mnum, &mat_mkl_cpardiso->mtype, &mat_mkl_cpardiso->phase, &mat_mkl_cpardiso->n, mat_mkl_cpardiso->a, mat_mkl_cpardiso->ia, mat_mkl_cpardiso->ja,
                                         mat_mkl_cpardiso->perm, &mat_mkl_cpardiso->nrhs, mat_mkl_cpardiso->iparm, &mat_mkl_cpardiso->msglvl, NULL, NULL, &mat_mkl_cpardiso->comm_mkl_cpardiso, (PetscInt *)&mat_mkl_cpardiso->err));
  PetscCheck(mat_mkl_cpardiso->err >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error reported by MKL Cluster PARDISO: err=%d, msg = \"%s\".Check manual", mat_mkl_cpardiso->err, Err_MSG_CPardiso(mat_mkl_cpardiso->err));

  mat_mkl_cpardiso->CleanUp     = PETSC_TRUE;
  F->ops->choleskyfactornumeric = MatFactorNumeric_MKL_CPARDISO;
  F->ops->solve                 = MatSolve_MKL_CPARDISO;
  F->ops->solvetranspose        = MatSolveTranspose_MKL_CPARDISO;
  F->ops->matsolve              = MatMatSolve_MKL_CPARDISO;
  if (A->spd == PETSC_BOOL3_TRUE) {
    F->ops->forwardsolve  = MatForwardSolve_MKL_CPARDISO;
    F->ops->backwardsolve = MatBackwardSolve_MKL_CPARDISO;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatView_MKL_CPARDISO(Mat A, PetscViewer viewer)
{
  PetscBool         isascii;
  PetscViewerFormat format;
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;
  PetscInt          i;

  PetscFunctionBegin;
  /* check if matrix is mkl_cpardiso type */
  if (A->ops->solve != MatSolve_MKL_CPARDISO) PetscFunctionReturn(PETSC_SUCCESS);

  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  if (isascii) {
    PetscCall(PetscViewerGetFormat(viewer, &format));
    if (format == PETSC_VIEWER_ASCII_INFO) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO run parameters:\n"));
      PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO phase:             %d \n", mat_mkl_cpardiso->phase));
      for (i = 1; i <= 64; i++) PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO iparm[%d]:     %d \n", i, mat_mkl_cpardiso->iparm[i - 1]));
      PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO maxfct:     %d \n", mat_mkl_cpardiso->maxfct));
      PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO mnum:     %d \n", mat_mkl_cpardiso->mnum));
      PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO mtype:     %d \n", mat_mkl_cpardiso->mtype));
      PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO n:     %d \n", mat_mkl_cpardiso->n));
      PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO nrhs:     %d \n", mat_mkl_cpardiso->nrhs));
      PetscCall(PetscViewerASCIIPrintf(viewer, "MKL Cluster PARDISO msglvl:     %d \n", mat_mkl_cpardiso->msglvl));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatGetInfo_MKL_CPARDISO(Mat A, MatInfoType flag, MatInfo *info)
{
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)A->data;

  PetscFunctionBegin;
  info->block_size        = 1.0;
  info->nz_allocated      = mat_mkl_cpardiso->nz + 0.0;
  info->nz_unneeded       = 0.0;
  info->assemblies        = 0.0;
  info->mallocs           = 0.0;
  info->memory            = 0.0;
  info->fill_ratio_given  = 0;
  info->fill_ratio_needed = 0;
  info->factor_mallocs    = 0;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode MatMkl_CPardisoSetCntl_MKL_CPARDISO(Mat F, PetscInt icntl, PetscInt ival)
{
  Mat_MKL_CPARDISO *mat_mkl_cpardiso = (Mat_MKL_CPARDISO *)F->data;

  PetscFunctionBegin;
  if (icntl <= 64) {
    mat_mkl_cpardiso->iparm[icntl - 1] = ival;
  } else {
    if (icntl == 65) mkl_set_num_threads((int)ival);
    else if (icntl == 66) mat_mkl_cpardiso->maxfct = ival;
    else if (icntl == 67) mat_mkl_cpardiso->mnum = ival;
    else if (icntl == 68) mat_mkl_cpardiso->msglvl = ival;
    else if (icntl == 69) mat_mkl_cpardiso->mtype = ival;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  MatMkl_CPardisoSetCntl - Set MKL Cluster PARDISO parameters
  <https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/onemkl-pardiso-parallel-direct-sparse-solver-iface.html>

  Logically Collective

  Input Parameters:
+ F     - the factored matrix obtained by calling `MatGetFactor()`
. icntl - index of MKL Cluster PARDISO parameter
- ival  - value of MKL Cluster PARDISO parameter

  Options Database Key:
. -mat_mkl_cpardiso_<icntl> <ival> - set the option numbered icntl to ival

  Level: intermediate

  Note:
  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
  database approach when working with `TS`, `SNES`, or `KSP`. See `MATSOLVERMKL_CPARDISO` for the options

.seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MATMPIAIJ`, `MATSOLVERMKL_CPARDISO`
@*/
PetscErrorCode MatMkl_CPardisoSetCntl(Mat F, PetscInt icntl, PetscInt ival)
{
  PetscFunctionBegin;
  PetscTryMethod(F, "MatMkl_CPardisoSetCntl_C", (Mat, PetscInt, PetscInt), (F, icntl, ival));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
  MATSOLVERMKL_CPARDISO -  A matrix type providing direct solvers (LU) for parallel matrices via the external package MKL Cluster PARDISO
  <https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/onemkl-pardiso-parallel-direct-sparse-solver-iface.html>

  Works with `MATMPIAIJ` matrices

  Use `-pc_type lu` `-pc_factor_mat_solver_type mkl_cpardiso` to use this direct solver

  Options Database Keys:
+ -mat_mkl_cpardiso_65 - Suggested number of threads to use within MKL Cluster PARDISO
. -mat_mkl_cpardiso_66 - Maximum number of factors with identical sparsity structure that must be kept in memory at the same time
. -mat_mkl_cpardiso_67 - Indicates the actual matrix for the solution phase
. -mat_mkl_cpardiso_68 - Message level information, use 1 to get detailed information on the solver options
. -mat_mkl_cpardiso_69 - Defines the matrix type. IMPORTANT: When you set this flag, iparm parameters are going to be set to the default ones for the matrix type
. -mat_mkl_cpardiso_1  - Use default values
. -mat_mkl_cpardiso_2  - Fill-in reducing ordering for the input matrix
. -mat_mkl_cpardiso_4  - Preconditioned CGS/CG
. -mat_mkl_cpardiso_5  - User permutation
. -mat_mkl_cpardiso_6  - Write solution on x
. -mat_mkl_cpardiso_8  - Iterative refinement step
. -mat_mkl_cpardiso_10 - Pivoting perturbation
. -mat_mkl_cpardiso_11 - Scaling vectors
. -mat_mkl_cpardiso_12 - Solve with transposed or conjugate transposed matrix A
. -mat_mkl_cpardiso_13 - Improved accuracy using (non-) symmetric weighted matching
. -mat_mkl_cpardiso_18 - Numbers of non-zero elements
. -mat_mkl_cpardiso_19 - Report number of floating point operations
. -mat_mkl_cpardiso_21 - Pivoting for symmetric indefinite matrices
. -mat_mkl_cpardiso_24 - Parallel factorization control
. -mat_mkl_cpardiso_25 - Parallel forward/backward solve control
. -mat_mkl_cpardiso_27 - Matrix checker
. -mat_mkl_cpardiso_31 - Partial solve and computing selected components of the solution vectors
. -mat_mkl_cpardiso_34 - Optimal number of threads for conditional numerical reproducibility (CNR) mode
- -mat_mkl_cpardiso_60 - Intel MKL Cluster PARDISO mode

  Level: beginner

  Notes:
  Use `-mat_mkl_cpardiso_68 1` to display the number of threads the solver is using. MKL does not provide a way to directly access this
  information.

  For more information on the options check
  <https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/onemkl-pardiso-parallel-direct-sparse-solver-iface.html>

.seealso: [](ch_matrices), `Mat`, `PCFactorSetMatSolverType()`, `MatSolverType`, `MatMkl_CPardisoSetCntl()`, `MatGetFactor()`, `MATSOLVERMKL_PARDISO`
M*/

static PetscErrorCode MatFactorGetSolverType_mkl_cpardiso(Mat A, MatSolverType *type)
{
  PetscFunctionBegin;
  *type = MATSOLVERMKL_CPARDISO;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* MatGetFactor for MPI AIJ matrices */
static PetscErrorCode MatGetFactor_mpiaij_mkl_cpardiso(Mat A, MatFactorType ftype, Mat *F)
{
  Mat               B;
  Mat_MKL_CPARDISO *mat_mkl_cpardiso;
  PetscBool         isSeqAIJ, isMPIBAIJ, isMPISBAIJ;

  PetscFunctionBegin;
  /* Create the factorization matrix */

  PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isSeqAIJ));
  PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIBAIJ, &isMPIBAIJ));
  PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPISBAIJ, &isMPISBAIJ));
  PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
  PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
  PetscCall(PetscStrallocpy("mkl_cpardiso", &((PetscObject)B)->type_name));
  PetscCall(MatSetUp(B));

  PetscCall(PetscNew(&mat_mkl_cpardiso));

  if (isSeqAIJ) mat_mkl_cpardiso->ConvertToTriples = MatCopy_seqaij_seqaij_MKL_CPARDISO;
  else if (isMPIBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpibaij_mpibaij_MKL_CPARDISO;
  else if (isMPISBAIJ) mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpisbaij_mpisbaij_MKL_CPARDISO;
  else mat_mkl_cpardiso->ConvertToTriples = MatConvertToTriples_mpiaij_mpiaij_MKL_CPARDISO;

  if (ftype == MAT_FACTOR_LU) B->ops->lufactorsymbolic = MatLUFactorSymbolic_AIJMKL_CPARDISO;
  else B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_AIJMKL_CPARDISO;
  B->ops->destroy = MatDestroy_MKL_CPARDISO;

  B->ops->view    = MatView_MKL_CPARDISO;
  B->ops->getinfo = MatGetInfo_MKL_CPARDISO;

  B->factortype = ftype;
  B->assembled  = PETSC_TRUE; /* required by -ksp_view */

  B->data = mat_mkl_cpardiso;

  /* set solvertype */
  PetscCall(PetscFree(B->solvertype));
  PetscCall(PetscStrallocpy(MATSOLVERMKL_CPARDISO, &B->solvertype));

  PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_mkl_cpardiso));
  PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMkl_CPardisoSetCntl_C", MatMkl_CPardisoSetCntl_MKL_CPARDISO));
  PetscCall(PetscInitialize_MKL_CPARDISO(A, mat_mkl_cpardiso));

  *F = B;
  PetscFunctionReturn(PETSC_SUCCESS);
}

PETSC_INTERN PetscErrorCode MatSolverTypeRegister_MKL_CPardiso(void)
{
  PetscFunctionBegin;
  PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
  PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
  PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPIBAIJ, MAT_FACTOR_LU, MatGetFactor_mpiaij_mkl_cpardiso));
  PetscCall(MatSolverTypeRegister(MATSOLVERMKL_CPARDISO, MATMPISBAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_mpiaij_mkl_cpardiso));
  PetscFunctionReturn(PETSC_SUCCESS);
}