blob: d30c8cec650ebe97aef7806a8aeb26dbf6897994 [file] [log] [blame]
// 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);
}
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)
(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);
}
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]++;
}