41 #ifndef BPP_NUMERIC_MATRIX_LUDECOMPOSITION_H
42 #define BPP_NUMERIC_MATRIX_LUDECOMPOSITION_H
45 #include "../../Exceptions.h"
46 #include "../NumTools.h"
81 std::vector<size_t>
piv;
86 size_t piv_length =
piv.size();
88 X.
resize(piv_length, j1 - j0 + 1);
90 for (
size_t i = 0; i < piv_length; i++)
92 for (
size_t j = j0; j <= j1; j++)
94 X(i, j - j0) = A(
piv[i], j);
99 static void permuteCopy(
const std::vector<Real>& A,
const std::vector<size_t>&
piv, std::vector<Real>& X)
101 size_t piv_length =
piv.size();
102 if (piv_length != A.size())
105 X.resize(piv_length);
107 for (
size_t i = 0; i < piv_length; i++)
199 L_(A.getNumberOfRows(), A.getNumberOfColumns()),
200 U_(A.getNumberOfColumns(), A.getNumberOfColumns()),
201 m(A.getNumberOfRows()),
202 n(A.getNumberOfColumns()),
204 piv(A.getNumberOfRows())
206 for (
size_t i = 0; i <
m; i++)
211 for (
size_t k = 0; k <
n; k++)
215 for (
size_t i = k + 1; i <
m; i++)
217 if (NumTools::abs<Real>(
LU(i, k)) > NumTools::abs<Real>(
LU(p, k)))
225 for (
size_t j = 0; j <
n; j++)
227 double t =
LU(p, j);
LU(p, j) =
LU(k, j);
LU(k, j) = t;
235 for (
size_t i = k + 1; i <
m; i++)
237 LU(i, k) /=
LU(k, k);
238 for (
size_t j = k + 1; j <
n; j++)
240 LU(i, j) -=
LU(i, k) *
LU(k, j);
254 for (
size_t i = 0; i <
m; i++)
256 for (
size_t j = 0; j <
n; j++)
282 for (
size_t i = 0; i <
n; i++)
284 for (
size_t j = 0; j <
n; j++)
322 for (
size_t j = 0; j <
n; j++)
349 Real minD = NumTools::abs<Real>(
LU(0, 0));
350 for (
size_t i = 1; i <
m; i++)
352 Real currentValue = NumTools::abs<Real>(
LU(i, i));
353 if (currentValue < minD)
368 for (
size_t k = 0; k <
n; k++)
370 for (
size_t i = k + 1; i <
n; i++)
372 for (
size_t j = 0; j < nx; j++)
374 X(i, j) -= X(k, j) *
LU(i, k);
386 for (
size_t j = 0; j < nx; j++)
390 for (
size_t i = 0; i < k; i++)
392 for (
size_t j = 0; j < nx; j++)
394 X(i, j) -= X(k, j) *
LU(i, k);
414 Real
solve (
const std::vector<Real>& b, std::vector<Real>& x)
const
423 Real minD = NumTools::abs<Real>(
LU(0, 0));
424 for (
size_t i = 1; i <
m; i++)
426 Real currentValue = NumTools::abs<Real>(
LU(i, i));
427 if (currentValue < minD)
439 for (
size_t k = 0; k <
n; k++)
441 for (
size_t i = k + 1; i <
n; i++)
443 x[i] -= x[k] *
LU(i, k);
454 for (
size_t i = 0; i < k; i++)
456 x[i] -= x[k] *
LU(i, k);
Number exception: integers.
const RowMatrix< Real > & getU()
Return upper triangular factor.
Real det() const
Compute determinant using LU factors.
LUDecomposition(const Matrix< Real > &A)
LU Decomposition.
std::vector< size_t > getPivot() const
Return pivot permutation vector.
std::vector< size_t > piv
Real solve(const std::vector< Real > &b, std::vector< Real > &x) const
Solve A*x = b, where x and b are vectors of length equal to the number of rows in A.
Real solve(const Matrix< Real > &B, Matrix< Real > &X) const
Solve A*X = B.
const RowMatrix< Real > & getL()
Return lower triangular factor.
static void permuteCopy(const Matrix< Real > &A, const std::vector< size_t > &piv, size_t j0, size_t j1, Matrix< Real > &X)
static void permuteCopy(const std::vector< Real > &A, const std::vector< size_t > &piv, std::vector< Real > &X)
The matrix template interface.
virtual size_t getNumberOfColumns() const =0
virtual void resize(size_t nRows, size_t nCols)=0
Resize the matrix.
virtual size_t getNumberOfRows() const =0
The base class exception for zero division error.