// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0

#include "common.h"

template <typename Index, typename Scalar, int StorageOrder, bool ConjugateLhs, bool ConjugateRhs>
struct general_matrix_vector_product_wrapper {
  static void run(Index rows, Index cols, const Scalar *lhs, Index lhsStride, const Scalar *rhs, Index rhsIncr,
                  Scalar *res, Index resIncr, Scalar alpha) {
    typedef Eigen::internal::const_blas_data_mapper<Scalar, Index, StorageOrder> LhsMapper;
    typedef Eigen::internal::const_blas_data_mapper<Scalar, Index, Eigen::RowMajor> RhsMapper;

    Eigen::internal::general_matrix_vector_product<Index, Scalar, LhsMapper, StorageOrder, ConjugateLhs, Scalar,
                                                   RhsMapper, ConjugateRhs>::run(rows, cols, LhsMapper(lhs, lhsStride),
                                                                                 RhsMapper(rhs, rhsIncr), res, resIncr,
                                                                                 alpha);
  }
};

EIGEN_BLAS_FUNC(gemv)
(const char *opa, const int *m, const int *n, const RealScalar *palpha, const RealScalar *pa, const int *lda,
 const RealScalar *pb, const int *incb, const RealScalar *pbeta, RealScalar *pc, const int *incc) {
  typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
  static const functype func[4] = {
      // array index: NOTR
      (general_matrix_vector_product_wrapper<int, Scalar, Eigen::ColMajor, false, false>::run),
      // array index: TR
      (general_matrix_vector_product_wrapper<int, Scalar, Eigen::RowMajor, false, false>::run),
      // array index: ADJ
      (general_matrix_vector_product_wrapper<int, Scalar, Eigen::RowMajor, Conj, false>::run), 0};

  const Scalar *a = reinterpret_cast<const Scalar *>(pa);
  const Scalar *b = reinterpret_cast<const Scalar *>(pb);
  Scalar *c = reinterpret_cast<Scalar *>(pc);
  Scalar alpha = *reinterpret_cast<const Scalar *>(palpha);
  Scalar beta = *reinterpret_cast<const Scalar *>(pbeta);

  // check arguments
  int info = 0;
  if (OP(*opa) == INVALID)
    info = 1;
  else if (*m < 0)
    info = 2;
  else if (*n < 0)
    info = 3;
  else if (*lda < std::max(1, *m))
    info = 6;
  else if (*incb == 0)
    info = 8;
  else if (*incc == 0)
    info = 11;
  if (info) return xerbla_(SCALAR_SUFFIX_UP "GEMV ", &info, kBlasNameLength);

  if (*m == 0 || *n == 0 || (alpha == Scalar(0) && beta == Scalar(1))) return;

  int actual_m = *m;
  int actual_n = *n;
  int code = OP(*opa);
  if (code != NOTR) std::swap(actual_m, actual_n);

  const Scalar *actual_b = get_compact_vector(b, actual_n, *incb);
  Scalar *actual_c = get_compact_vector(c, actual_m, *incc);

  if (beta != Scalar(1)) {
    if (beta == Scalar(0))
      make_vector(actual_c, actual_m).setZero();
    else
      make_vector(actual_c, actual_m) *= beta;
  }

  if (code >= 4 || func[code] == 0) return;

  func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha);

  if (actual_b != b) delete[] actual_b;
  if (actual_c != c) delete[] copy_back(actual_c, c, actual_m, *incc);
}

EIGEN_BLAS_FUNC(trsv)
(const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda,
 RealScalar *pb, const int *incb) {
  typedef void (*functype)(int, const Scalar *, int, Scalar *);
  using Eigen::ColMajor;
  using Eigen::Lower;
  using Eigen::OnTheLeft;
  using Eigen::RowMajor;
  using Eigen::UnitDiag;
  using Eigen::Upper;
  static const functype func[16] = {
      // array index: NOTR  | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, false, ColMajor>::run),
      // array index: TR    | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, false, RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, Conj, RowMajor>::run), 0,
      // array index: NOTR  | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, false, ColMajor>::run),
      // array index: TR    | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, false, RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, Conj, RowMajor>::run), 0,
      // array index: NOTR  | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, false,
                                                ColMajor>::run),
      // array index: TR    | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, false,
                                                RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, Conj, RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, false,
                                                ColMajor>::run),
      // array index: TR    | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, false,
                                                RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, Conj, RowMajor>::run),
      0};

  const Scalar *a = reinterpret_cast<const Scalar *>(pa);
  Scalar *b = reinterpret_cast<Scalar *>(pb);

  int info = 0;
  if (UPLO(*uplo) == INVALID)
    info = 1;
  else if (OP(*opa) == INVALID)
    info = 2;
  else if (DIAG(*diag) == INVALID)
    info = 3;
  else if (*n < 0)
    info = 4;
  else if (*lda < std::max(1, *n))
    info = 6;
  else if (*incb == 0)
    info = 8;
  if (info) return xerbla_(SCALAR_SUFFIX_UP "TRSV ", &info, kBlasNameLength);

  Scalar *actual_b = get_compact_vector(b, *n, *incb);

  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
  func[code](*n, a, *lda, actual_b);

  if (actual_b != b) delete[] copy_back(actual_b, b, *n, *incb);
}

EIGEN_BLAS_FUNC(trmv)
(const char *uplo, const char *opa, const char *diag, const int *n, const RealScalar *pa, const int *lda,
 RealScalar *pb, const int *incb) {
  typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, const Scalar &);
  using Eigen::ColMajor;
  using Eigen::Lower;
  using Eigen::OnTheLeft;
  using Eigen::RowMajor;
  using Eigen::UnitDiag;
  using Eigen::Upper;
  static const functype func[16] = {
      // array index: NOTR  | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Upper | 0, Scalar, false, Scalar, false, ColMajor>::run),
      // array index: TR    | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Lower | 0, Scalar, false, Scalar, false, RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Lower | 0, Scalar, Conj, Scalar, false, RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Lower | 0, Scalar, false, Scalar, false, ColMajor>::run),
      // array index: TR    | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Upper | 0, Scalar, false, Scalar, false, RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Upper | 0, Scalar, Conj, Scalar, false, RowMajor>::run),
      0,
      // array index: NOTR  | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, false, Scalar, false,
                                                         ColMajor>::run),
      // array index: TR    | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, false, Scalar, false,
                                                         RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, Conj, Scalar, false,
                                                         RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, false, Scalar, false,
                                                         ColMajor>::run),
      // array index: TR    | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, false, Scalar, false,
                                                         RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, Conj, Scalar, false,
                                                         RowMajor>::run),
      0};

  const Scalar *a = reinterpret_cast<const Scalar *>(pa);
  Scalar *b = reinterpret_cast<Scalar *>(pb);

  int info = 0;
  if (UPLO(*uplo) == INVALID)
    info = 1;
  else if (OP(*opa) == INVALID)
    info = 2;
  else if (DIAG(*diag) == INVALID)
    info = 3;
  else if (*n < 0)
    info = 4;
  else if (*lda < std::max(1, *n))
    info = 6;
  else if (*incb == 0)
    info = 8;
  if (info) return xerbla_(SCALAR_SUFFIX_UP "TRMV ", &info, kBlasNameLength);

  if (*n == 0) return;

  Scalar *actual_b = get_compact_vector(b, *n, *incb);
  Eigen::Matrix<Scalar, Eigen::Dynamic, 1> res(*n);
  res.setZero();

  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
  if (code >= 16 || func[code] == 0) return;

  func[code](*n, *n, a, *lda, actual_b, 1, res.data(), 1, Scalar(1));

  copy_back(res.data(), b, *n, *incb);
  if (actual_b != b) delete[] actual_b;
}

/**  GBMV  performs one of the matrix-vector operations
 *
 *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
 *
 *  where alpha and beta are scalars, x and y are vectors and A is an
 *  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
 */
EIGEN_BLAS_FUNC(gbmv)
(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx,
 RealScalar *pbeta, RealScalar *py, int *incy) {
  const Scalar *a = reinterpret_cast<const Scalar *>(pa);
  const Scalar *x = reinterpret_cast<const Scalar *>(px);
  Scalar *y = reinterpret_cast<Scalar *>(py);
  Scalar alpha = *reinterpret_cast<const Scalar *>(palpha);
  Scalar beta = *reinterpret_cast<const Scalar *>(pbeta);
  int coeff_rows = *kl + *ku + 1;

  int info = 0;
  if (OP(*trans) == INVALID)
    info = 1;
  else if (*m < 0)
    info = 2;
  else if (*n < 0)
    info = 3;
  else if (*kl < 0)
    info = 4;
  else if (*ku < 0)
    info = 5;
  else if (*lda < coeff_rows)
    info = 8;
  else if (*incx == 0)
    info = 10;
  else if (*incy == 0)
    info = 13;
  if (info) return xerbla_(SCALAR_SUFFIX_UP "GBMV ", &info, kBlasNameLength);

  if (*m == 0 || *n == 0 || (alpha == Scalar(0) && beta == Scalar(1))) return;

  int actual_m = *m;
  int actual_n = *n;
  if (OP(*trans) != NOTR) std::swap(actual_m, actual_n);

  const Scalar *actual_x = get_compact_vector(x, actual_n, *incx);
  Scalar *actual_y = get_compact_vector(y, actual_m, *incy);

  if (beta != Scalar(1)) {
    if (beta == Scalar(0))
      make_vector(actual_y, actual_m).setZero();
    else
      make_vector(actual_y, actual_m) *= beta;
  }

  ConstMatrixType mat_coeffs(a, coeff_rows, *n, *lda);

  int nb = std::min(*n, (*m) + (*ku));
  for (int j = 0; j < nb; ++j) {
    int start = std::max(0, j - *ku);
    int end = std::min((*m) - 1, j + *kl);
    int len = end - start + 1;
    int offset = (*ku) - j + start;
    if (OP(*trans) == NOTR)
      make_vector(actual_y + start, len) += (alpha * actual_x[j]) * mat_coeffs.col(j).segment(offset, len);
    else if (OP(*trans) == TR)
      actual_y[j] +=
          alpha * (mat_coeffs.col(j).segment(offset, len).transpose() * make_vector(actual_x + start, len)).value();
    else
      actual_y[j] +=
          alpha * (mat_coeffs.col(j).segment(offset, len).adjoint() * make_vector(actual_x + start, len)).value();
  }

  if (actual_x != x) delete[] actual_x;
  if (actual_y != y) delete[] copy_back(actual_y, y, actual_m, *incy);
}

/**  TBMV  performs one of the matrix-vector operations
 *
 *     x := A*x,   or   x := A'*x,   or   x := conjg(A')*x,
 *
 *  where x is an n element vector and A is an n by n unit, or non-unit,
 *  upper or lower triangular band matrix, with ( k + 1 ) diagonals.
 *
 *  Band storage: upper triangle stores A[i,j] at a[(k+i-j) + j*lda],
 *  lower triangle stores A[i,j] at a[(i-j) + j*lda].
 */
EIGEN_BLAS_FUNC(tbmv)
(char *uplo, char *opa, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx) {
  Scalar *a = reinterpret_cast<Scalar *>(pa);
  Scalar *x = reinterpret_cast<Scalar *>(px);

  int info = 0;
  if (UPLO(*uplo) == INVALID)
    info = 1;
  else if (OP(*opa) == INVALID)
    info = 2;
  else if (DIAG(*diag) == INVALID)
    info = 3;
  else if (*n < 0)
    info = 4;
  else if (*k < 0)
    info = 5;
  else if (*lda < *k + 1)
    info = 7;
  else if (*incx == 0)
    info = 9;
  if (info) return xerbla_(SCALAR_SUFFIX_UP "TBMV ", &info, kBlasNameLength);

  if (*n == 0) return;

  Scalar *actual_x = get_compact_vector(x, *n, *incx);

  bool upper = (UPLO(*uplo) == UP);
  int op = OP(*opa);
  bool unit = (DIAG(*diag) == UNIT);

  if (*k >= 8) {
    // Vectorized path: use Eigen Map segments for the inner band operations.
    ConstMatrixType band(a, *k + 1, *n, *lda);
    if (op == NOTR) {
      if (upper) {
        for (int j = 0; j < *n; ++j) {
          if (actual_x[j] != Scalar(0)) {
            int start = std::max(0, j - *k);
            int len = j - start;
            int offset = *k - (j - start);
            Scalar temp = actual_x[j];
            if (len > 0) make_vector(actual_x + start, len) += temp * band.col(j).segment(offset, len);
            if (!unit) actual_x[j] = temp * band(*k, j);
          }
        }
      } else {
        for (int j = *n - 1; j >= 0; --j) {
          if (actual_x[j] != Scalar(0)) {
            int len = std::min(*n - 1, j + *k) - j;
            Scalar temp = actual_x[j];
            if (len > 0) make_vector(actual_x + j + 1, len) += temp * band.col(j).segment(1, len);
            if (!unit) actual_x[j] = temp * band(0, j);
          }
        }
      }
    } else if (op == TR) {
      if (upper) {
        for (int j = *n - 1; j >= 0; --j) {
          int start = std::max(0, j - *k);
          int len = j - start;
          int offset = *k - (j - start);
          Scalar temp = actual_x[j];
          if (!unit) temp *= band(*k, j);
          if (len > 0)
            temp += (band.col(j).segment(offset, len).cwiseProduct(make_vector(actual_x + start, len))).sum();
          actual_x[j] = temp;
        }
      } else {
        for (int j = 0; j < *n; ++j) {
          int len = std::min(*n - 1, j + *k) - j;
          Scalar temp = actual_x[j];
          if (!unit) temp *= band(0, j);
          if (len > 0) temp += (band.col(j).segment(1, len).cwiseProduct(make_vector(actual_x + j + 1, len))).sum();
          actual_x[j] = temp;
        }
      }
    } else {
      // Conjugate transpose: .dot() computes conj(lhs) . rhs.
      if (upper) {
        for (int j = *n - 1; j >= 0; --j) {
          int start = std::max(0, j - *k);
          int len = j - start;
          int offset = *k - (j - start);
          Scalar temp = actual_x[j];
          if (!unit) temp *= Eigen::numext::conj(band(*k, j));
          if (len > 0) temp += band.col(j).segment(offset, len).dot(make_vector(actual_x + start, len));
          actual_x[j] = temp;
        }
      } else {
        for (int j = 0; j < *n; ++j) {
          int len = std::min(*n - 1, j + *k) - j;
          Scalar temp = actual_x[j];
          if (!unit) temp *= Eigen::numext::conj(band(0, j));
          if (len > 0) temp += band.col(j).segment(1, len).dot(make_vector(actual_x + j + 1, len));
          actual_x[j] = temp;
        }
      }
    }
  } else {
    // Scalar path: for narrow bandwidth, avoid Map overhead.
    if (op == NOTR) {
      if (upper) {
        for (int j = 0; j < *n; ++j) {
          if (actual_x[j] != Scalar(0)) {
            Scalar temp = actual_x[j];
            for (int i = std::max(0, j - *k); i < j; ++i) actual_x[i] += temp * a[(*k + i - j) + j * *lda];
            if (!unit) actual_x[j] = temp * a[*k + j * *lda];
          }
        }
      } else {
        for (int j = *n - 1; j >= 0; --j) {
          if (actual_x[j] != Scalar(0)) {
            Scalar temp = actual_x[j];
            for (int i = j + 1; i <= std::min(*n - 1, j + *k); ++i) actual_x[i] += temp * a[(i - j) + j * *lda];
            if (!unit) actual_x[j] = temp * a[j * *lda];
          }
        }
      }
    } else {
      // Transpose or conjugate transpose.
      auto maybe_conj = [op](Scalar val) -> Scalar { return op == ADJ ? Eigen::numext::conj(val) : val; };
      if (upper) {
        for (int j = *n - 1; j >= 0; --j) {
          Scalar temp = actual_x[j];
          if (!unit) temp *= maybe_conj(a[*k + j * *lda]);
          for (int i = std::max(0, j - *k); i < j; ++i) temp += maybe_conj(a[(*k + i - j) + j * *lda]) * actual_x[i];
          actual_x[j] = temp;
        }
      } else {
        for (int j = 0; j < *n; ++j) {
          Scalar temp = actual_x[j];
          if (!unit) temp *= maybe_conj(a[j * *lda]);
          for (int i = j + 1; i <= std::min(*n - 1, j + *k); ++i)
            temp += maybe_conj(a[(i - j) + j * *lda]) * actual_x[i];
          actual_x[j] = temp;
        }
      }
    }
  }

  if (actual_x != x) delete[] copy_back(actual_x, x, *n, *incx);
}

/**  DTBSV  solves one of the systems of equations
 *
 *     A*x = b,   or   A'*x = b,
 *
 *  where b and x are n element vectors and A is an n by n unit, or
 *  non-unit, upper or lower triangular band matrix, with ( k + 1 )
 *  diagonals.
 *
 *  No test for singularity or near-singularity is included in this
 *  routine. Such tests must be performed before calling this routine.
 */
EIGEN_BLAS_FUNC(tbsv)
(char *uplo, char *op, char *diag, int *n, int *k, RealScalar *pa, int *lda, RealScalar *px, int *incx) {
  typedef void (*functype)(int, int, const Scalar *, int, Scalar *);
  using Eigen::ColMajor;
  using Eigen::Lower;
  using Eigen::OnTheLeft;
  using Eigen::RowMajor;
  using Eigen::UnitDiag;
  using Eigen::Upper;
  static const functype func[16] = {
      // array index: NOTR  | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Upper | 0, Scalar, false, Scalar, ColMajor>::run),
      // array index: TR    | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Lower | 0, Scalar, false, Scalar, RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Lower | 0, Scalar, Conj, Scalar, RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Lower | 0, Scalar, false, Scalar, ColMajor>::run),
      // array index: TR    | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Upper | 0, Scalar, false, Scalar, RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Upper | 0, Scalar, Conj, Scalar, RowMajor>::run),
      0,
      // array index: NOTR  | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Upper | UnitDiag, Scalar, false, Scalar, ColMajor>::run),
      // array index: TR    | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Lower | UnitDiag, Scalar, false, Scalar, RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Lower | UnitDiag, Scalar, Conj, Scalar, RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Lower | UnitDiag, Scalar, false, Scalar, ColMajor>::run),
      // array index: TR    | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Upper | UnitDiag, Scalar, false, Scalar, RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::band_solve_triangular_selector<int, Upper | UnitDiag, Scalar, Conj, Scalar, RowMajor>::run),
      0,
  };

  Scalar *a = reinterpret_cast<Scalar *>(pa);
  Scalar *x = reinterpret_cast<Scalar *>(px);
  int coeff_rows = *k + 1;

  int info = 0;
  if (UPLO(*uplo) == INVALID)
    info = 1;
  else if (OP(*op) == INVALID)
    info = 2;
  else if (DIAG(*diag) == INVALID)
    info = 3;
  else if (*n < 0)
    info = 4;
  else if (*k < 0)
    info = 5;
  else if (*lda < coeff_rows)
    info = 7;
  else if (*incx == 0)
    info = 9;
  if (info) return xerbla_(SCALAR_SUFFIX_UP "TBSV ", &info, kBlasNameLength);

  if (*n == 0 || (*k == 0 && DIAG(*diag) == UNIT)) return;

  int actual_n = *n;

  Scalar *actual_x = get_compact_vector(x, actual_n, *incx);

  int code = OP(*op) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
  if (code >= 16 || func[code] == 0) return;

  func[code](*n, *k, a, *lda, actual_x);

  if (actual_x != x) delete[] copy_back(actual_x, x, actual_n, *incx);
}

/**  DTPMV  performs one of the matrix-vector operations
 *
 *     x := A*x,   or   x := A'*x,
 *
 *  where x is an n element vector and  A is an n by n unit, or non-unit,
 *  upper or lower triangular matrix, supplied in packed form.
 */
EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) {
  typedef void (*functype)(int, const Scalar *, const Scalar *, Scalar *, Scalar);
  using Eigen::ColMajor;
  using Eigen::Lower;
  using Eigen::OnTheLeft;
  using Eigen::RowMajor;
  using Eigen::UnitDiag;
  using Eigen::Upper;
  static const functype func[16] = {
      // array index: NOTR  | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | 0, Scalar, false, Scalar, false,
                                                                ColMajor>::run),
      // array index: TR    | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | 0, Scalar, false, Scalar, false,
                                                                RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | 0, Scalar, Conj, Scalar, false,
                                                                RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | 0, Scalar, false, Scalar, false,
                                                                ColMajor>::run),
      // array index: TR    | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | 0, Scalar, false, Scalar, false,
                                                                RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | 0, Scalar, Conj, Scalar, false,
                                                                RowMajor>::run),
      0,
      // array index: NOTR  | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, false, Scalar, false,
                                                                ColMajor>::run),
      // array index: TR    | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, false, Scalar, false,
                                                                RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, Conj, Scalar, false,
                                                                RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Lower | UnitDiag, Scalar, false, Scalar, false,
                                                                ColMajor>::run),
      // array index: TR    | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, false, Scalar, false,
                                                                RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_matrix_vector_product<int, Upper | UnitDiag, Scalar, Conj, Scalar, false,
                                                                RowMajor>::run),
      0};

  Scalar *ap = reinterpret_cast<Scalar *>(pap);
  Scalar *x = reinterpret_cast<Scalar *>(px);

  int info = 0;
  if (UPLO(*uplo) == INVALID)
    info = 1;
  else if (OP(*opa) == INVALID)
    info = 2;
  else if (DIAG(*diag) == INVALID)
    info = 3;
  else if (*n < 0)
    info = 4;
  else if (*incx == 0)
    info = 7;
  if (info) return xerbla_(SCALAR_SUFFIX_UP "TPMV ", &info, kBlasNameLength);

  if (*n == 0) return;

  Scalar *actual_x = get_compact_vector(x, *n, *incx);
  Eigen::Matrix<Scalar, Eigen::Dynamic, 1> res(*n);
  res.setZero();

  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
  if (code >= 16 || func[code] == 0) return;

  func[code](*n, ap, actual_x, res.data(), Scalar(1));

  copy_back(res.data(), x, *n, *incx);
  if (actual_x != x) delete[] actual_x;
}

/**  DTPSV  solves one of the systems of equations
 *
 *     A*x = b,   or   A'*x = b,
 *
 *  where b and x are n element vectors and A is an n by n unit, or
 *  non-unit, upper or lower triangular matrix, supplied in packed form.
 *
 *  No test for singularity or near-singularity is included in this
 *  routine. Such tests must be performed before calling this routine.
 */
EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pap, RealScalar *px, int *incx) {
  typedef void (*functype)(int, const Scalar *, Scalar *);
  using Eigen::ColMajor;
  using Eigen::Lower;
  using Eigen::OnTheLeft;
  using Eigen::RowMajor;
  using Eigen::UnitDiag;
  using Eigen::Upper;
  static const functype func[16] = {
      // array index: NOTR  | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, false,
                                                       ColMajor>::run),
      // array index: TR    | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, false,
                                                       RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, Conj, RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | 0, false,
                                                       ColMajor>::run),
      // array index: TR    | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, false,
                                                       RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (NUNIT << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | 0, Conj, RowMajor>::run),
      0,
      // array index: NOTR  | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, false,
                                                       ColMajor>::run),
      // array index: TR    | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, false,
                                                       RowMajor>::run),
      // array index: ADJ   | (UP << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, Conj,
                                                       RowMajor>::run),
      0,
      // array index: NOTR  | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Lower | UnitDiag, false,
                                                       ColMajor>::run),
      // array index: TR    | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, false,
                                                       RowMajor>::run),
      // array index: ADJ   | (LO << 2) | (UNIT  << 3)
      (Eigen::internal::packed_triangular_solve_vector<Scalar, Scalar, int, OnTheLeft, Upper | UnitDiag, Conj,
                                                       RowMajor>::run),
      0};

  Scalar *ap = reinterpret_cast<Scalar *>(pap);
  Scalar *x = reinterpret_cast<Scalar *>(px);

  int info = 0;
  if (UPLO(*uplo) == INVALID)
    info = 1;
  else if (OP(*opa) == INVALID)
    info = 2;
  else if (DIAG(*diag) == INVALID)
    info = 3;
  else if (*n < 0)
    info = 4;
  else if (*incx == 0)
    info = 7;
  if (info) return xerbla_(SCALAR_SUFFIX_UP "TPSV ", &info, kBlasNameLength);

  Scalar *actual_x = get_compact_vector(x, *n, *incx);

  int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
  func[code](*n, ap, actual_x);

  if (actual_x != x) delete[] copy_back(actual_x, x, *n, *incx);
}