mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-04-12 19:20:36 +08:00
83 lines
2.7 KiB
C++
83 lines
2.7 KiB
C++
// 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);
|
|
|
|
using Eigen::UnitLower;
|
|
using Eigen::Upper;
|
|
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]++;
|
|
}
|