lapack/lu.inc

// 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/.
// SPDX-License-Identifier: 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, kBlasNameLength);
  }

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

  Scalar *a = reinterpret_cast<Scalar *>(pa);
  int nb_transpositions;
  int ret = int(Eigen::internal::partial_lu_impl<Scalar, Eigen::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;
}

// 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)
(const 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, kBlasNameLength);
  }

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

  using Eigen::UnitLower;
  using Eigen::Upper;
  // IPIV is input-only (and may be shared by concurrent solves): convert to 0-based indices in a copy.
  Eigen::Transpositions<Eigen::Dynamic, Eigen::Dynamic, int> pivots(*n);
  for (int i = 0; i < *n; ++i) pivots.indices()(i) = ipiv[i] - 1;
  if (OP(*trans) == NOTR) {
    B = pivots * 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 = pivots.transpose() * B;
  } else if (OP(*trans) == ADJ) {
    lu.triangularView<Upper>().adjoint().solveInPlace(B);
    lu.triangularView<UnitLower>().adjoint().solveInPlace(B);
    B = pivots.transpose() * B;
  }
}