tahoma2d/toonz/sources/tnzext/tlin/tlin_superlu_wrap.cpp

extern "C" {
#include "slu_ddefs.h"
#include "slu_util.h"
}

#include <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();
  const tcg::list<spmat::HashMap::BucketNode> &nodes = m.entries().items();

  std::vector<const spmat::HashMap::BucketNode *> colEntries;
  std::vector<const spmat::HashMap::BucketNode *>::size_type j, size;

  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) {
  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;
  dgstrf(opt, &AC, sp_ienv(1), sp_ienv(2), etree, NULL, 0, F->perm_c, F->perm_r,
         F->L, F->U, &stat, &result);

  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);
  double *_x      = const_cast<double *>(x);

  sp_dgemv("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);

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