LU and LU inverse is done
This commit is contained in:
+143
-26
@@ -7,55 +7,172 @@ namespace decomp{
|
||||
|
||||
// Stores PA = LU with partial pivoting (row permutations).
|
||||
template <typename T>
|
||||
struct LU{
|
||||
utils::Matrix<T> LUmat; // combined L (unit diagonal implied) and U
|
||||
std::vector<uint64_t> piv; // pivot indices (row permutations)
|
||||
bool singular= false; // set true if zero (or near-zero) pivots encountered
|
||||
struct LUdcmp{
|
||||
uint64_t rows; // Stores number of rows.
|
||||
utils::Matrix<T> lu; // Stores the decomposition.
|
||||
std::vector<uint64_t> indx; // Stores the permutation.
|
||||
T d; // Used by det.
|
||||
|
||||
LU() = default;
|
||||
// Default Constructor
|
||||
LUdcmp() = default;
|
||||
|
||||
explicit LU(const utils::Matrix<T>& A) {
|
||||
factor(A);
|
||||
// Constructor
|
||||
LUdcmp(const utils::Matrix<T>& A){
|
||||
decomp(A);
|
||||
}
|
||||
|
||||
void factor(const utils::Matrix<T>&A){
|
||||
void decomp(const utils::Matrix<T>&A){
|
||||
|
||||
rows = A.rows();
|
||||
|
||||
const uint64_t rows = A.rows();
|
||||
if (rows != A.cols()){
|
||||
throw std::runtime_error("LU: factor non-square");
|
||||
throw std::runtime_error("LUdcmp: decomp non-square");
|
||||
}
|
||||
|
||||
if (rows == 0){
|
||||
LUmat = A;
|
||||
piv.clear();
|
||||
singular = false;
|
||||
return;
|
||||
}
|
||||
uint64_t imax{0};
|
||||
lu = A;
|
||||
indx.resize(rows);
|
||||
std::vector<T> vv(rows); // vv stores the implicit scaling of each row.
|
||||
T big{T{0}}, tmp{T{0}};// TINY{T{1.0e-40}};
|
||||
|
||||
T big, tmp;
|
||||
d = T{1}; // No row interchanges yet.
|
||||
|
||||
LUmat = A;
|
||||
piv.resize(rows); // piv stores the implicit scaling of each row.
|
||||
//double d = 1.0; // No row interchanges yet.
|
||||
|
||||
for (uint64_t i = 0; i < rows; ++i){ // Loop over rows to get the implicit scaling information.
|
||||
// Loop over rows to get the implicit scaling information.
|
||||
for (uint64_t i = 0; i < rows; ++i){
|
||||
big = T{0};
|
||||
for (uint64_t j = 0; j < rows; ++j){
|
||||
tmp=std::abs(LUmat(i,j));
|
||||
tmp = std::abs(lu(i,j));
|
||||
if (tmp > big){
|
||||
big = tmp;
|
||||
}
|
||||
}
|
||||
if (big == T{0}){
|
||||
throw std::runtime_error("Singular matrix in LU.factor()");
|
||||
throw std::runtime_error("LUdcmp: Singular matrix");
|
||||
}
|
||||
// No nonzero largest element.
|
||||
vv[i] = T{1}/big; // Save the scaling.
|
||||
}
|
||||
// This is the outermost kij loop. Initialize for the search for largest pivot element.
|
||||
for (uint64_t k = 0; k < rows; ++k){
|
||||
big = T{0};
|
||||
imax = k;
|
||||
for (uint64_t i = k; i < rows; ++i){
|
||||
tmp = vv[i] * static_cast<T>(std::fabs(static_cast<double>(lu(i,k))));
|
||||
if (tmp > big){ // Is the figure of merit for the pivot better than the best so far?
|
||||
big = tmp;
|
||||
imax = i;
|
||||
}
|
||||
}
|
||||
if (k != imax){ // Do we need to interchange rows?
|
||||
lu.swap_rows(imax, k); // Yes, do so...
|
||||
d = -d; // ...and change the parity of d.
|
||||
vv[imax] = vv[k]; // Also interchange the scale factor.
|
||||
}
|
||||
indx[k] = imax;
|
||||
if (lu(k,k) == T{0.0}){ // if the pivot element is zero, the matrix is singular (at least to the precision of thealgorithm).
|
||||
throw std::runtime_error("LUdcmp: Singular matrix");
|
||||
//lu(k,k) = TINY; // For some applications on singular matrices, it is desirable to substitute TINY for zero.
|
||||
}
|
||||
for (uint64_t i = k+1; i < rows; ++i){
|
||||
tmp = lu(i,k) /= lu(k,k); // Divide by the pivot element.
|
||||
for (uint64_t j = k+1; j < rows; ++j){ // Innermost loop: reduce remaining submatrix.
|
||||
lu(i,j) -= tmp*lu(k,j);
|
||||
}
|
||||
}
|
||||
}
|
||||
} // end void decomp(const utils::Matrix<T>&A)
|
||||
|
||||
// Solves the set of n linear equations A*x=b using the stored LU decomposition of A.
|
||||
void inplace_solve(const utils::Vector<T>& b, utils::Vector<T>& x){
|
||||
T sum{T{0}};
|
||||
|
||||
uint64_t ii{0}, ip{0};
|
||||
|
||||
if (b.size() != rows || x.size() != rows){
|
||||
throw std::runtime_error("LUdcmp: inplace_solve bad sizes");
|
||||
}
|
||||
x = b;
|
||||
|
||||
for (uint64_t i = 0; i < rows; ++i){ // When ii is set to a positive value, it will become the index of the first nonvanishing element of b.
|
||||
ip = indx[i];
|
||||
sum = x[ip];
|
||||
x[ip] = x[i];
|
||||
if (ii >= 0){
|
||||
for (uint64_t j = ii; j < i; ++j){
|
||||
sum -= lu(i,j)*x[j];
|
||||
}
|
||||
}else if (sum != T{0}) { // A nonzero element was encountered, so from now on we will have to do the sums in the loop above.
|
||||
ii = i+1;
|
||||
}
|
||||
x[i] = sum;
|
||||
}
|
||||
for (int64_t i = static_cast<int64_t>(rows)-1; i >= 0; --i){ // Now we do the backsubstitution.
|
||||
sum=x[i];
|
||||
for (uint64_t j = static_cast<uint64_t>(i)+1; j < rows; ++j){
|
||||
sum -= lu(static_cast<uint64_t>(i),j)*x[j];
|
||||
}
|
||||
x[static_cast<uint64_t>(i)] = sum/lu(static_cast<uint64_t>(i),static_cast<uint64_t>(i)); // Store a component of the solution vector x.
|
||||
}
|
||||
} // end inplace_solve(utils::Vector<T>& b, utils::Vector<T>& x)
|
||||
|
||||
// SSolves m sets of n linear equations A*X=B using the stored LU decomposition of A.
|
||||
void inplace_solve(const utils::Matrix<T>& B, utils::Matrix<T>& X) {
|
||||
|
||||
uint64_t m{B.cols()};
|
||||
|
||||
if (B.rows() != rows || X.rows() != rows || B.cols() != X.cols()){
|
||||
throw std::runtime_error("LUdcmp: inplace_solve bad sizes");
|
||||
}
|
||||
|
||||
utils::Vector<T> xx(rows);
|
||||
|
||||
for (uint64_t j = 0; j < m; ++j){ // Copy and solve each column in turn.
|
||||
|
||||
xx = B.get_col(j);
|
||||
inplace_solve(xx,xx);
|
||||
X.set_col(j, xx);
|
||||
}
|
||||
|
||||
} // end inplace_solve(utils::Matrix<T>& B, utils::Matrix<T>& X)
|
||||
|
||||
// Solves the set of n linear equations A*x=b using the stored LU decomposition of A.
|
||||
utils::Vector<T> solve(const utils::Vector<T>& b) {
|
||||
utils::Vector<T> x(rows,T{0});
|
||||
inplace_solve(b, x);
|
||||
return x;
|
||||
}
|
||||
|
||||
// Solves the set of n linear equations A*X=B using the stored LU decomposition of A.
|
||||
utils::Matrix<T> solve(const utils::Matrix<T>& B) {
|
||||
utils::Matrix<T> X(rows,B.cols(),T{0});
|
||||
inplace_solve(B, X);
|
||||
return X;
|
||||
}
|
||||
|
||||
void inplace_inverse(utils::Matrix<T>& Ainv){
|
||||
Ainv.eye(rows);
|
||||
inplace_solve(Ainv, Ainv);
|
||||
}
|
||||
|
||||
utils::Matrix<T> inverse(){
|
||||
utils::Matrix<T> Ainv;
|
||||
inplace_inverse(Ainv);
|
||||
return Ainv;
|
||||
}
|
||||
|
||||
T det(){
|
||||
T dd = d;
|
||||
for (uint64_t i = 0; i < rows; ++i){
|
||||
dd *= lu(i,i);
|
||||
}
|
||||
return dd;
|
||||
}
|
||||
|
||||
|
||||
}; // struct LU
|
||||
|
||||
typedef LU<float> LUf;
|
||||
typedef LU<double> LUd;
|
||||
typedef LUdcmp<float> LUdcmpf;
|
||||
typedef LUdcmp<double> LUdcmpd;
|
||||
|
||||
|
||||
} // namespace decomp
|
||||
@@ -23,8 +23,11 @@ namespace numerics{
|
||||
if (method == "Gauss-Jordan"){
|
||||
inverse_gj(A);
|
||||
}
|
||||
else if(method == "LU"){
|
||||
inplace_inverse_lu(A);
|
||||
}
|
||||
else{
|
||||
throw std::runtime_error("numerics::inplace_inverse(" + method + ") - Not implemented yet \r \nImplemented: 'Gauss-Jordan',");
|
||||
throw std::runtime_error("numerics::inplace_inverse(" + method + ") - Not implemented yet \r \nImplemented: 'Gauss-Jordan', 'LU'");
|
||||
}
|
||||
}
|
||||
|
||||
@@ -32,21 +35,11 @@ namespace numerics{
|
||||
|
||||
template <typename T>
|
||||
utils::Matrix<T> inverse(utils::Matrix<T>& A, std::string method = "Gauss-Jordan"){
|
||||
|
||||
|
||||
utils::Matrix<T> B = A;
|
||||
|
||||
inplace_inverse(B, method);
|
||||
|
||||
return B;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
} // namespace numerics
|
||||
|
||||
#endif // _inverse_n_
|
||||
@@ -7,12 +7,13 @@
|
||||
|
||||
#include <omp.h>
|
||||
|
||||
|
||||
namespace numerics{
|
||||
|
||||
template <typename T>
|
||||
void inverse_gj(utils::Matrix<T>& A){
|
||||
utils::Matrix<T> B(A.rows(),A.cols(), T{0});
|
||||
//utils::Matrix<T> B(A.rows(),A.cols(), T{0});
|
||||
utils::Matrix<T> B;
|
||||
B.eye(A.rows());
|
||||
|
||||
|
||||
uint64_t icol{0}, irow{0}, rows{A.rows()}, cols{A.cols()};
|
||||
@@ -36,6 +37,9 @@ namespace numerics{
|
||||
}
|
||||
}
|
||||
}
|
||||
if (big <= T{1e-14}){
|
||||
throw std::runtime_error("utill:inplace_inverse('Gauss-Jordan' - Singular Matrix");
|
||||
}
|
||||
ipiv[icol]++;
|
||||
if (irow != icol){
|
||||
for (uint64_t l = 0; l < rows; l++){ // SWAP
|
||||
|
||||
@@ -0,0 +1,20 @@
|
||||
#pragma once
|
||||
|
||||
#include "./decomp/lu.h"
|
||||
|
||||
|
||||
|
||||
|
||||
namespace numerics{
|
||||
|
||||
template <typename T>
|
||||
void inplace_inverse_lu(utils::Matrix<T>& A){
|
||||
if (A.rows() != A.cols()){
|
||||
throw std::runtime_error("numerics inverse_lu: non-square matrix");
|
||||
}
|
||||
|
||||
decomp::LUdcmp<T> lu(A);
|
||||
lu.inplace_inverse(A);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
+17
-1
@@ -45,7 +45,7 @@ public:
|
||||
return !(*this == A);
|
||||
}
|
||||
bool nearly_equal(const Matrix<T>& A, T tol = static_cast<T>(1e-9)) const {
|
||||
if (rows_ != A.rows_ || cols_ != A.cols_) return false;
|
||||
if (rows_ != A.rows() || cols_ != A.cols()) return false;
|
||||
for (uint64_t i = 0; i < rows_; ++i)
|
||||
for (uint64_t j = 0; j < cols_; ++j) {
|
||||
T a = (*this)(i,j);
|
||||
@@ -59,6 +59,22 @@ public:
|
||||
return true;
|
||||
}
|
||||
|
||||
void eye(uint64_t rows_cols){
|
||||
rows_ = cols_ = rows_cols;
|
||||
|
||||
data_.clear();
|
||||
data_.resize(rows_cols*rows_cols, T{0});
|
||||
for (uint64_t i = 0; i < rows_; ++i){
|
||||
data_[i * cols_ + i] = T{1};
|
||||
}
|
||||
}
|
||||
|
||||
void resize(uint64_t rows, uint64_t cols, const T& value = T(0)){
|
||||
rows_ = rows;
|
||||
cols_ = cols;
|
||||
data_.resize(rows_*cols_, value);
|
||||
|
||||
}
|
||||
|
||||
|
||||
//# MATRIX: row helpers (copy out) #
|
||||
|
||||
+27
-2
@@ -27,6 +27,7 @@ public:
|
||||
}
|
||||
|
||||
|
||||
|
||||
//##########################################################
|
||||
//# VECTOR: --- basic properties --- #
|
||||
//##########################################################
|
||||
@@ -399,13 +400,37 @@ public:
|
||||
void print() const{
|
||||
std::cout << *this << std::endl;
|
||||
}
|
||||
|
||||
|
||||
bool nearly_equal(const Vector<T>& a, T tol = static_cast<T>(1e-9)) const {
|
||||
if (v.size() != a.size()){
|
||||
return false;
|
||||
}
|
||||
for (uint64_t i = 0; i < v.size(); ++i){
|
||||
T val1 = v[i];
|
||||
T val2 = a[i];
|
||||
if (std::is_floating_point<T>::value) {
|
||||
if (std::fabs(val1 - val2) > tol) return false;
|
||||
} else {
|
||||
if (val1 != val2) return false;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
typedef Vector<int> Vi;
|
||||
typedef Vector<float> Vf;
|
||||
typedef Vector<double> Vd;
|
||||
|
||||
|
||||
/*
|
||||
if (std::is_floating_point<T>::value) {
|
||||
if (std::fabs(a - b) > tol) return false;
|
||||
} else {
|
||||
if (a != b) return false;
|
||||
}*/
|
||||
} // namespace utils
|
||||
|
||||
#endif // _vector_n_
|
||||
Reference in New Issue
Block a user