extern "C" {

#include <slu_ddefs.h></slu_ddefs.h>

#include <slu_util.h></slu_util.h>

}

#include <algorithm></algorithm>

#include "tlin/tlin_cblas_wrap.h"

#include "tlin/tlin_superlu_wrap.h"

//*************************************************************************

//    Preliminaries

//*************************************************************************

struct tlin::SuperMatrix final : public ::SuperMatrix {};

struct tlin::superlu_options_t final : public ::superlu_options_t {};

//=======================================================================

namespace {

static const tlin::spmat::HashMap::size_t neg = -1;

bool initialized = false;

static tlin::superlu_options_t defaultOpt;

struct DefaultOptsInitializer {

  DefaultOptsInitializer() {

    set_default_options(&defaultOpt);

    defaultOpt.PrintStat = NO;

  }

} _instance;

inline bool rowLess(const tlin::spmat::HashMap::BucketNode *a,

                    const tlin::spmat::HashMap::BucketNode *b) {

  return a->m_key.first < b->m_key.first;

}

}

//*************************************************************************

//    SuperLU-specific Functions

//*************************************************************************

void tlin::allocS(SuperMatrix *&A, int rows, int cols, int nnz) {

  A              = (SuperMatrix *)SUPERLU_MALLOC(sizeof(SuperMatrix));

  double *values = doubleMalloc(nnz);

  int *rowind    = intMalloc(nnz);

  int *colptr    = intMalloc(cols + 1);

  dCreate_CompCol_Matrix(A, rows, cols, nnz, values, rowind, colptr, SLU_NC,

                         SLU_D, SLU_GE);

}

//---------------------------------------------------------------

void tlin::allocD(SuperMatrix *&A, int rows, int cols) {

  A = (SuperMatrix *)SUPERLU_MALLOC(sizeof(SuperMatrix));

  double *values = doubleMalloc(rows * cols * sizeof(double));

  dCreate_Dense_Matrix(A, rows, cols, values, rows, SLU_DN, SLU_D, SLU_GE);

}

//---------------------------------------------------------------

void tlin::allocS(SuperMatrix *&A, int rows, int cols, int nnz, int *colptr,

                  int *rowind, double *values) {

  A = (SuperMatrix *)SUPERLU_MALLOC(sizeof(SuperMatrix));

  dCreate_CompCol_Matrix(A, rows, cols, nnz, values, rowind, colptr, SLU_NC,

                         SLU_D, SLU_GE);

}

//---------------------------------------------------------------

void tlin::allocD(SuperMatrix *&A, int rows, int cols, int lda,

                  double *values) {

  A = (SuperMatrix *)SUPERLU_MALLOC(sizeof(SuperMatrix));

  dCreate_Dense_Matrix(A, rows, cols, values, lda, SLU_DN, SLU_D, SLU_GE);

}

//---------------------------------------------------------------

void tlin::freeS(SuperMatrix *m) {

  if (!m) return;

  Destroy_CompCol_Matrix(m);

  SUPERLU_FREE(m);

}

//---------------------------------------------------------------

void tlin::freeD(SuperMatrix *m) {

  if (!m) return;

  Destroy_Dense_Matrix(m);

  SUPERLU_FREE(m);

}

//---------------------------------------------------------------

void tlin::createS(SuperMatrix &A, int rows, int cols, int nnz) {

  double *values = doubleMalloc(nnz);

  int *rowind    = intMalloc(nnz);

  int *colptr    = intMalloc(cols + 1);

  dCreate_CompCol_Matrix(&A, rows, cols, nnz, values, rowind, colptr, SLU_NC,

                         SLU_D, SLU_GE);

}

//---------------------------------------------------------------

void tlin::createD(SuperMatrix &A, int rows, int cols) {

  double *values = doubleMalloc(rows * cols * sizeof(double));

  dCreate_Dense_Matrix(&A, rows, cols, values, rows, SLU_DN, SLU_D, SLU_GE);

}

//---------------------------------------------------------------

void tlin::createS(SuperMatrix &A, int rows, int cols, int nnz, int *colptr,

                   int *rowind, double *values) {

  dCreate_CompCol_Matrix(&A, rows, cols, nnz, values, rowind, colptr, SLU_NC,

                         SLU_D, SLU_GE);

}

//---------------------------------------------------------------

void tlin::createD(SuperMatrix &A, int rows, int cols, int lda,

                   double *values) {

  dCreate_Dense_Matrix(&A, rows, cols, values, lda, SLU_DN, SLU_D, SLU_GE);

}

//---------------------------------------------------------------

void tlin::destroyS(SuperMatrix &A, bool destroyData) {

  if (destroyData)

    Destroy_CompCol_Matrix(&A);

  else

    SUPERLU_FREE(A.Store);

}

//---------------------------------------------------------------

void tlin::destroyD(SuperMatrix &A, bool destroyData) {

  if (destroyData)

    Destroy_Dense_Matrix(&A);

  else

    SUPERLU_FREE(A.Store);

}

//---------------------------------------------------------------

void tlin::freeF(SuperFactors *F) {

  if (!F) return;

  Destroy_SuperNode_Matrix(F->L);

  Destroy_CompCol_Matrix(F->U);

  SUPERLU_FREE(F->L);

  SUPERLU_FREE(F->U);

  SUPERLU_FREE(F->perm_r);

  SUPERLU_FREE(F->perm_c);

  SUPERLU_FREE(F);

}

//---------------------------------------------------------------

void tlin::readDN(SuperMatrix *A, int &lda, double *&values) {

  assert(A->Stype == SLU_DN);

  DNformat *storage = (DNformat *)A->Store;

  lda    = storage->lda;

  values = (double *)storage->nzval;

}

//---------------------------------------------------------------

void tlin::readNC(SuperMatrix *A, int &nnz, int *&colptr, int *&rowind,

                  double *&values) {

  assert(A->Stype == SLU_NC);  // Only SLU_NC (CCS) format is supported here

  NCformat *storage = (NCformat *)A->Store;

  nnz    = storage->nnz;

  values = (double *)storage->nzval;

  rowind = (int *)storage->rowind;

  colptr = (int *)storage->colptr;

}

//---------------------------------------------------------------

void tlin::traduceD(const mat &m, SuperMatrix *&A) {

  int rows = m.rows(), cols = m.cols();

  if (!A) allocD(A, rows, cols);

  int lda;

  double *Avalues = 0;

  readDN(A, lda, Avalues);

  assert(A->nrow == rows && A->ncol == cols && lda == rows);

  memcpy(Avalues, m.values(), rows * cols * sizeof(double));

}

//---------------------------------------------------------------

void tlin::traduceS(spmat &m, SuperMatrix *&A) {

  int rows = m.rows(), cols = m.cols(), nnz = (int)m.entries().size();

  spmat::HashMap &entries = m.entries();

  // Build or extract pointers to out's data

  double *values;

  int *rowind, *colptr;

  if (!A) allocS(A, rows, cols, nnz);

  // Retrieve NC arrays from A

  int Annz;

  readNC(A, Annz, colptr, rowind, values);

  assert(A->nrow == rows && A->ncol == cols && Annz == nnz);

  // Rehash to cols buckets

  if (entries.hashFunctor().m_cols != cols) entries.hashFunctor().m_cols = cols;

  entries.rehash(cols);

  // Copy each bucket to the corresponding col

  const std::vector<spmat::hashmap::size_t> &buckets = m.entries().buckets();</spmat::hashmap::size_t>

  const tcg::list<spmat::hashmap::bucketnode> &nodes = m.entries().items();</spmat::hashmap::bucketnode>

  std::vector<const *="" spmat::hashmap::bucketnode=""> colEntries;</const>

  std::vector<const *="" spmat::hashmap::bucketnode="">::size_type j, size;</const>

  double *currVal = values;

  int *currRowind = rowind;

  for (int i = 0; i < cols; ++i) {

    colptr[i]                      = (int)(currVal - values);

    spmat::HashMap::size_t nodeIdx = buckets[i];

    // Retrieve all column entry pointers.

    colEntries.clear();

    while (nodeIdx != neg) {

      const spmat::HashMap::BucketNode &node = nodes[nodeIdx];

      colEntries.push_back(&node);

      nodeIdx = nodes[nodeIdx].m_next;

    }

    // Sort them by row

    std::sort(colEntries.begin(), colEntries.end(), rowLess);

    // Finally, write them in the SuperMatrix

    size = colEntries.size();

    for (j = 0; j < size; ++j) {

      const spmat::HashMap::BucketNode &node = *colEntries[j];

      *currRowind = node.m_key.first;

      *currVal    = node.m_val;

      ++currVal, ++currRowind;

    }

  }

  colptr[cols] = nnz;

}

//---------------------------------------------------------------

void tlin::traduceD(const tlin::sparse_matrix<double> &m, SuperMatrix *&A) {</double>

  int rows = m.rows(), cols = m.cols();

  const spmat::HashMap &entries = m.entries();

  // Build or extract pointers to out's data

  double *values;

  if (!A) allocD(A, rows, cols);

  // Retrieve DN arrays from A

  int lda;

  readDN(A, lda, values);

  assert(A->nrow == rows && A->ncol == cols && lda == rows);

  // Copy each value in entries to A

  spmat::HashMap::const_iterator it;

  for (it = entries.begin(); it != entries.end(); ++it)

    values[it->m_key.second * rows + it->m_key.first] = it->m_val;

}

//---------------------------------------------------------------

void tlin::factorize(SuperMatrix *A, SuperFactors *&F, superlu_options_t *opt) {

  assert(A->nrow == A->ncol);

  int n = A->nrow;

  if (!F) F = (SuperFactors *)SUPERLU_MALLOC(sizeof(SuperFactors));

  if (!opt) opt = &defaultOpt;

  F->perm_c = intMalloc(n);

  get_perm_c(3, A, F->perm_c);

  SuperMatrix AC;

  int *etree = intMalloc(n);

  sp_preorder(opt, A, F->perm_c, etree, &AC);

  F->L      = (SuperMatrix *)SUPERLU_MALLOC(sizeof(SuperMatrix));

  F->U      = (SuperMatrix *)SUPERLU_MALLOC(sizeof(SuperMatrix));

  F->perm_r = intMalloc(n);

  SuperLUStat_t stat;

  StatInit(&stat);

  int result;

#if defined(SUPERLU_MAJOR_VERSION) && (SUPERLU_MAJOR_VERSION >= 5)

  GlobalLU_t Glu;

  dgstrf(opt, &AC, sp_ienv(1), sp_ienv(2), etree, NULL, 0, F->perm_c, F->perm_r,

         F->L, F->U, &Glu, &stat, &result);

#else

  dgstrf(opt, &AC, sp_ienv(1), sp_ienv(2), etree, NULL, 0, F->perm_c, F->perm_r,

         F->L, F->U, &stat, &result);

#endif

  StatFree(&stat);

  Destroy_CompCol_Permuted(&AC);

  SUPERLU_FREE(etree);

  if (result != 0) freeF(F), F = 0;

}

//---------------------------------------------------------------

void tlin::solve(SuperFactors *F, SuperMatrix *BX, superlu_options_t *opt) {

  assert(F);

  if (!opt) opt = &defaultOpt;

  SuperLUStat_t stat;

  StatInit(&stat);

  int result;

  dgstrs(NOTRANS, F->L, F->U, F->perm_c, F->perm_r, BX, &stat, &result);

  StatFree(&stat);

}

//---------------------------------------------------------------

void tlin::solve(SuperFactors *F, SuperMatrix *B, SuperMatrix *&X,

                 superlu_options_t *opt) {

  if (!X) allocD(X, B->nrow, B->ncol);

  double *Bvalues = 0, *Xvalues = 0;

  int lda;

  readDN(B, lda, Bvalues);

  readDN(X, lda, Xvalues);

  memcpy(Xvalues, Bvalues, B->nrow * B->ncol * sizeof(double));

  solve(F, X, opt);

}

//---------------------------------------------------------------

void tlin::solve(SuperMatrix *A, SuperMatrix *BX, superlu_options_t *opt) {

  assert(A->nrow == A->ncol);

  int n = A->nrow;

  if (!opt) opt = &defaultOpt;

  SuperMatrix L, U;

  int *perm_c, *perm_r;

  perm_c = intMalloc(n);

  perm_r = intMalloc(n);

  SuperLUStat_t stat;

  StatInit(&stat);

  int result;

  dgssv(opt, A, perm_c, perm_r, &L, &U, BX, &stat, &result);

  Destroy_SuperNode_Matrix(&L);

  Destroy_CompCol_Matrix(&U);

  SUPERLU_FREE(perm_r);

  SUPERLU_FREE(perm_c);

  StatFree(&stat);

}

//---------------------------------------------------------------

void tlin::solve(SuperMatrix *A, SuperMatrix *B, SuperMatrix *&X,

                 superlu_options_t *opt) {

  if (!X) allocD(X, B->nrow, B->ncol);

  double *Bvalues = 0, *Xvalues = 0;

  int lda;

  readDN(B, lda, Bvalues);

  readDN(X, lda, Xvalues);

  memcpy(Xvalues, Bvalues, B->nrow * B->ncol * sizeof(double));

  solve(A, X, opt);

}

//---------------------------------------------------------------

void tlin::solve(SuperFactors *F, double *bx, superlu_options_t *opt) {

  SuperMatrix BX;

  int rows = F->L->nrow;

  createD(BX, rows, 1, rows, bx);

  tlin::solve(F, &BX, opt);

  SUPERLU_FREE(BX.Store);

}

//---------------------------------------------------------------

void tlin::solve(SuperFactors *F, double *b, double *&x,

                 superlu_options_t *opt) {

  SuperMatrix B, X;

  int rows = F->L->nrow;

  if (!x) x = (double *)malloc(rows * sizeof(double));

  createD(B, rows, 1, rows, b);

  createD(X, rows, 1, rows, x);

  SuperMatrix *Xptr = &X;  //&X is const

  tlin::solve(F, &B, Xptr, opt);

  SUPERLU_FREE(B.Store);

  SUPERLU_FREE(X.Store);

}

//---------------------------------------------------------------

void tlin::solve(SuperMatrix *A, double *bx, superlu_options_t *opt) {

  SuperMatrix BX;

  int rows = A->nrow;

  createD(BX, rows, 1, rows, bx);

  tlin::solve(A, &BX, opt);

  SUPERLU_FREE(BX.Store);

}

//---------------------------------------------------------------

void tlin::solve(SuperMatrix *A, double *b, double *&x,

                 superlu_options_t *opt) {

  SuperMatrix B, X;

  int rows = A->nrow;

  if (!x) x = (double *)malloc(rows * sizeof(double));

  createD(B, rows, 1, rows, b);

  createD(X, rows, 1, rows, x);

  SuperMatrix *Xptr = &X;  //&X is const

  tlin::solve(A, &B, Xptr, opt);

  SUPERLU_FREE(B.Store);

  SUPERLU_FREE(X.Store);

}

//*************************************************************************

//    BLAS-related Functions

//*************************************************************************

void tlin::multiplyS(const SuperMatrix *A, const double *x, double *&y) {

  /*

int sp_dgemv (char *, double, SuperMatrix *, double *,

                  int, double, double *, int);

*/

  if (!y) {

    y = (double *)malloc(A->nrow * sizeof(double));

    memset(y, 0, A->nrow * sizeof(double));

  }

  SuperMatrix *_A = const_cast<supermatrix *="">(A);</supermatrix>

  double *_x      = const_cast<double *="">(x);</double>

  sp_dgemv((char *)"N", 1.0, _A, _x, 1, 1.0, y, 1);

}

//---------------------------------------------------------------

void tlin::multiplyD(const SuperMatrix *A, const double *x, double *&y) {

  int lda;

  double *values;

  tlin::readDN(const_cast<supermatrix *="">(A), lda, values);</supermatrix>

  tlin::multiply(A->nrow, A->ncol, values, x, y);

}