// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010-2011 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/.

#include "common.h"
#include <Eigen/LU>

// computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges
EIGEN_LAPACK_FUNC(getrf, (int *m, int *n, RealScalar *pa, int *lda, int *ipiv, int *info)) {
  *info = 0;
  if (*m < 0)
    *info = -1;
  else if (*n < 0)
    *info = -2;
  else if (*lda < std::max(1, *m))
    *info = -4;
  if (*info != 0) {
    int e = -*info;
    return xerbla_(SCALAR_SUFFIX_UP "GETRF", &e, 6);
  }

  if (*m == 0 || *n == 0) return 0;

  Scalar *a = reinterpret_cast<Scalar *>(pa);
  int nb_transpositions;
  int ret = int(
      Eigen::internal::partial_lu_impl<Scalar, ColMajor, int>::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions));

  for (int i = 0; i < std::min(*m, *n); ++i) ipiv[i]++;

  if (ret >= 0) *info = ret + 1;

  return 0;
}

// GETRS solves a system of linear equations
//     A * X = B  or  A' * X = B
//   with a general N-by-N matrix A using the LU factorization computed  by GETRF
EIGEN_LAPACK_FUNC(getrs, (char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb,
                          int *info)) {
  *info = 0;
  if (OP(*trans) == INVALID)
    *info = -1;
  else if (*n < 0)
    *info = -2;
  else if (*nrhs < 0)
    *info = -3;
  else if (*lda < std::max(1, *n))
    *info = -5;
  else if (*ldb < std::max(1, *n))
    *info = -8;
  if (*info != 0) {
    int e = -*info;
    return xerbla_(SCALAR_SUFFIX_UP "GETRS", &e, 6);
  }

  Scalar *a = reinterpret_cast<Scalar *>(pa);
  Scalar *b = reinterpret_cast<Scalar *>(pb);
  MatrixType lu(a, *n, *n, *lda);
  MatrixType B(b, *n, *nrhs, *ldb);

  for (int i = 0; i < *n; ++i) ipiv[i]--;
  if (OP(*trans) == NOTR) {
    B = PivotsType(ipiv, *n) * B;
    lu.triangularView<UnitLower>().solveInPlace(B);
    lu.triangularView<Upper>().solveInPlace(B);
  } else if (OP(*trans) == TR) {
    lu.triangularView<Upper>().transpose().solveInPlace(B);
    lu.triangularView<UnitLower>().transpose().solveInPlace(B);
    B = PivotsType(ipiv, *n).transpose() * B;
  } else if (OP(*trans) == ADJ) {
    lu.triangularView<Upper>().adjoint().solveInPlace(B);
    lu.triangularView<UnitLower>().adjoint().solveInPlace(B);
    B = PivotsType(ipiv, *n).transpose() * B;
  }
  for (int i = 0; i < *n; ++i) ipiv[i]++;

  return 0;
}