# Algorithms_in_C++: numerical_methods/lu_decomposition.h Source File

lu_decomposition.h
Go to the documentation of this file.
/**
* @file lu_decomposition.h
* @author [Krishna Vedala](https://github.com/kvedala)
* @brief Functions associated with [LU
* Decomposition](https://en.wikipedia.org/wiki/LU_decomposition)
* of a square matrix.
*/
#pragma once
#include <iostream>
#include <valarray>
#include <vector>
#ifdef _OPENMP
#include <omp.h>
#endif
/** Define matrix type as a `std::vector` of `std::valarray` */
template <typename T>
using matrix = std::vector<std::valarray<T>>;
/** Perform LU decomposition on matrix
* \param[in] A matrix to decompose
* \param[out] L output L matrix
* \param[out] U output U matrix
* \returns 0 if no errors
* \returns negative if error occurred
*/
template <typename T>
int lu_decomposition(const matrix<T> &A, matrix<double> *L, matrix<double> *U) {
int row, col, j;
int mat_size = A.size();
if (mat_size != A[0].size()) {
// check matrix is a square matrix
std::cerr << "Not a square matrix!\n";
return -1;
}
// regularize each row
for (row = 0; row < mat_size; row++) {
// Upper triangular matrix
#ifdef _OPENMP
#pragma omp for
#endif
for (col = row; col < mat_size; col++) {
// Summation of L[i,j] * U[j,k]
double lu_sum = 0.;
for (j = 0; j < row; j++) {
lu_sum += L[0][row][j] * U[0][j][col];
}
// Evaluate U[i,k]
U[0][row][col] = A[row][col] - lu_sum;
}
// Lower triangular matrix
#ifdef _OPENMP
#pragma omp for
#endif
for (col = row; col < mat_size; col++) {
if (row == col) {
L[0][row][col] = 1.;
continue;
}
// Summation of L[i,j] * U[j,k]
double lu_sum = 0.;
for (j = 0; j < row; j++) {
lu_sum += L[0][col][j] * U[0][j][row];
}
// Evaluate U[i,k]
L[0][col][row] = (A[col][row] - lu_sum) / U[0][row][row];
}
}
return 0;
}
/**
* Compute determinant of an NxN square matrix using LU decomposition.
* Using LU decomposition, the determinant is given by the product of diagonal
* elements of matrices L and U.
*
* @tparam T datatype of input matrix - int, unsigned int, double, etc
* @param A input square matrix
* @return determinant of matrix A
*/
template <typename T>
double determinant_lu(const matrix<T> &A) {
matrix<double> L(A.size(), std::valarray<double>(A.size()));
matrix<double> U(A.size(), std::valarray<double>(A.size()));
if (lu_decomposition(A, &L, &U) < 0)
return 0;
double result = 1.f;
for (size_t i = 0; i < A.size(); i++) {
result *= L[i][i] * U[i][i];
}
return result;
}
check_termination
bool check_termination(long double delta)
Definition: durand_kerner_roots.cpp:91
Sqrt
double Sqrt(double a)
Definition: sqrt_double.cpp:16
std::srand
T srand(T... args)
tests
void tests()
Definition: number_of_positive_divisors.cpp:70
lcm
unsigned int lcm(unsigned int x, unsigned int y)
Definition: least_common_multiple.cpp:43
tests
void tests()
Definition: least_common_multiple.cpp:50
prime_numbers
std::vector< int > prime_numbers
Definition: prime_factorization.cpp:16
main
int main()
Definition: primes_up_to_billion.cpp:26
MAX_ITERATIONS
#define MAX_ITERATIONS
Maximum number of iterations to check.
Definition: false_position.cpp:24
statistics::stats_computer2::operator>>
friend std::istream & operator>>(std::istream &input, stats_computer2 &stat)
Definition: realtime_stats.cpp:97
mat_size
ll mat_size
Definition: matrix_exponentiation.cpp:45
EPSILON
#define EPSILON
solution accuracy limit
Definition: golden_search_extrema.cpp:17
test3
void test3()
Test function to find maxima for the function in the interval Expected result: .
Definition: brent_method_extrema.cpp:188
std::string
STL class.
std::uniform_int_distribution
primes
std::vector< int > primes(int max)
Definition: prime_numbers.cpp:12
miller_rabin_primality_test
bool miller_rabin_primality_test(T num, T repeats)
Definition: miller_rabin.cpp:125
std::clock_t
std::pair
main
int main()
Definition: brent_method_extrema.cpp:204
std::cos
T cos(T... args)
EPSILON
#define EPSILON
system accuracy limit
Definition: brent_method_extrema.cpp:23
std::vector
STL class.
test2
void test2()
Definition: lu_decompose.cpp:66
std::vector::size
T size(T... args)
multiply
int multiply(int x, int res[], int res_size)
Definition: power_for_huge_numbers.cpp:25
ans
ll ans(ll n)
Definition: matrix_exponentiation.cpp:91
main
int main()
Definition: sieve_of_eratosthenes.cpp:65
sieve
std::vector< bool > sieve(uint32_t N)
Definition: sieve_of_eratosthenes.cpp:26
complex_str
const char * complex_str(const std::complex< double > &x)
Definition: durand_kerner_roots.cpp:76
test3
void test3()
Test function to find maxima for the function in the interval Expected result: .
Definition: golden_search_extrema.cpp:123
get_minima
double get_minima(const std::function< double(double)> &f, double lim_a, double lim_b)
Get the real root of a function in the given interval.
Definition: brent_method_extrema.cpp:35
std::setfill
T setfill(T... args)
std::function
sgn
int sgn(T val)
Definition: bisection_method.cpp:32
statistics::stats_computer1::std
double std() const
Definition: realtime_stats.cpp:48
factors
std::vector< std::pair< int, int > > factors
Definition: prime_factorization.cpp:19
main
int main()
Definition: prime_factorization.cpp:62
reverse_binary
std::vector< T > reverse_binary(T num)
Definition: miller_rabin.cpp:19
std::sqrt
T sqrt(T... args)
power
void power(int x, int n)
Definition: power_for_huge_numbers.cpp:56
test1
void test1()
Test function to find minima for the function in the interval Expected result = 2.
Definition: golden_search_extrema.cpp:78
is_prime
bool is_prime(T num)
Definition: check_prime.cpp:22
main
int main()
Definition: prime_numbers.cpp:26
std::mt19937
std::vector::push_back
T push_back(T... args)
main
int main()
Definition: modular_inverse_fermat_little_theorem.cpp:84
std::clock
T clock(T... args)
durand_kerner_algo
std::pair< uint32_t, double > durand_kerner_algo(const std::valarray< double > &coeffs, std::valarray< std::complex< double >> *roots, bool write_log=false)
Definition: durand_kerner_roots.cpp:109
test1
void test1()
Definition: lu_decompose.cpp:36
test
void test()
Definition: sum_of_digits.cpp:58
std::isnan
T isnan(T... args)
statistics
Statistical algorithms.
test2
void test2()
Definition: sum_of_digits.cpp:49
EPSILON
constexpr double EPSILON
system accuracy limit
Definition: newton_raphson_method.cpp:20
std::snprintf
T snprintf(T... args)
std::random_device
std::numeric_limits::infinity
T infinity(T... args)
std::cout
std::ofstream
STL class.
main
int main()
Definition: false_position.cpp:39
std::isinf
T isinf(T... args)
statistics::stats_computer2::variance
double variance() const
Definition: realtime_stats.cpp:89
gcd
unsigned int gcd(unsigned int x, unsigned int y)
Definition: least_common_multiple.cpp:16
main
int main()
Definition: sqrt_double.cpp:42
modular_exponentiation
T modular_exponentiation(T base, const std::vector< T > &rev_binary_exponent, T mod)
Definition: miller_rabin.cpp:43
poly_function
std::complex< double > poly_function(const std::valarray< double > &coeffs, std::complex< double > x)
Definition: durand_kerner_roots.cpp:53
main
int main(int argc, char **argv)
Definition: lu_decompose.cpp:84
std::valarray
STL class.
std::complex::real
T real(T... args)
test1
void test1()
Definition: durand_kerner_roots.cpp:207
main
int main()
Definition: gaussian_elimination.cpp:9
std::ofstream::open
T open(T... args)
ACCURACY
#define ACCURACY
Definition: durand_kerner_roots.cpp:45
statistics::stats_computer1::mean
double mean() const
Definition: realtime_stats.cpp:42
statistics::stats_computer2
Definition: realtime_stats.cpp:72
main
int main(int argc, char **argv)
Definition: realtime_stats.cpp:158
sgn
int sgn(T val)
Definition: false_position.cpp:34
print
void print(uint32_t N, const std::vector< bool > &is_prime)
Definition: sieve_of_eratosthenes.cpp:44
std::rand
T rand(T... args)
std::swap
T swap(T... args)
lu_decomposition.h
Functions associated with LU Decomposition of a square matrix.
test_function
void test_function(const float *test_data, const int number_of_samples)
Definition: realtime_stats.cpp:118
SieveOfEratosthenes
void SieveOfEratosthenes(int N)
Definition: prime_factorization.cpp:23
std::string::substr
T substr(T... args)
miller_test
bool miller_test(T d, T num)
Definition: miller_rabin.cpp:73
main
int main()
Definition: least_common_multiple.cpp:78
lu_decomposition
int lu_decomposition(const matrix< T > &A, matrix< double > *L, matrix< double > *U)
Definition: lu_decomposition.h:29
binExpo
int64_t binExpo(int64_t a, int64_t b, int64_t m)
Definition: modular_inverse_fermat_little_theorem.cpp:52
std::endl
T endl(T... args)
statistics::stats_computer1::new_val
void new_val(T x)
Definition: realtime_stats.cpp:32
test2
void test2()
Test function to find root for the function in the interval Expected result: .
Definition: brent_method_extrema.cpp:165
std::left
T left(T... args)
statistics::stats_computer2::new_val
void new_val(T x)
Definition: realtime_stats.cpp:77
statistics::stats_computer1::variance
double variance() const
Definition: realtime_stats.cpp:45
sum_of_digits
int sum_of_digits(int num)
Definition: sum_of_digits.cpp:23
main
int main()
Definition: number_of_positive_divisors.cpp:81
std::vector::cbegin
T cbegin(T... args)
main
int main()
Definition: miller_rabin.cpp:183
std::strtod
T strtod(T... args)
main
int main()
Definition: power_for_huge_numbers.cpp:82
prime_factorization
void prime_factorization(int num)
Definition: prime_factorization.cpp:40
main
int main()
Definition: sum_of_digits.cpp:68
determinant_lu
double determinant_lu(const matrix< T > &A)
Definition: lu_decomposition.h:90
std::fixed
T fixed(T... args)
statistics::stats_computer1
Definition: realtime_stats.cpp:27
test2
void test2()
Test function to find maxima for the function in the interval Expected result: .
Definition: golden_search_extrema.cpp:100
statistics::stats_computer2::std
double std() const
Definition: realtime_stats.cpp:92
main
int main()
Definition: golden_search_extrema.cpp:139
test1
void test1()
Definition: sum_of_digits.cpp:40
std::complex
STL class.
std::complex::imag
T imag(T... args)
eq
static double eq(double i)
Definition: false_position.cpp:28
test1
void test1()
Test function to find root for the function in the interval Expected result = 2.
Definition: brent_method_extrema.cpp:143
add
std::string add(std::string a, std::string b)
Definition: string_fibonacci.cpp:24
std::make_pair
T make_pair(T... args)
std::time
T time(T... args)
MAX_ITERATIONS
#define MAX_ITERATIONS
Maximum number of iterations to check.
Definition: bisection_method.cpp:22
std::vector::cend
T cend(T... args)
tests
void tests()
Definition: sieve_of_eratosthenes.cpp:56
std::setw
T setw(T... args)
std::max
T max(T... args)
main
int main()
Definition: string_fibonacci.cpp:81
number_of_positive_divisors
int number_of_positive_divisors(int n)
Definition: number_of_positive_divisors.cpp:33
eq
static double eq(double i)
Definition: bisection_method.cpp:26
get_minima
double get_minima(const std::function< double(double)> &f, double lim_a, double lim_b)
Get the minima of a function in the given interval. To get the maxima, simply negate the function....
Definition: golden_search_extrema.cpp:29
fib_Accurate
void fib_Accurate(uint64_t n)
Definition: string_fibonacci.cpp:68
Sieve
void Sieve(int64_t n)
Definition: primes_up_to_billion.cpp:13
std::cin
operator<<
std::ostream & operator<<(std::ostream &out, matrix< T > const &v)
Definition: lu_decompose.cpp:18
statistics::stats_computer2::mean
double mean() const
Definition: realtime_stats.cpp:86
isPrime
bool isPrime(int64_t m)
Definition: modular_inverse_fermat_little_theorem.cpp:68
std::ofstream::is_open
T is_open(T... args)
test2
void test2()
Definition: durand_kerner_roots.cpp:242
main
int main()
Definition: bisection_method.cpp:37
std::exit
T exit(T... args)
MAX
#define MAX
Definition: power_for_huge_numbers.cpp:10
statistics::stats_computer1::operator>>
friend std::istream & operator>>(std::istream &input, stats_computer1 &stat)
Definition: realtime_stats.cpp:53
tests
void tests()
Definition: miller_rabin.cpp:157
main
int main()
Definition: graph_coloring.cpp:96
prime
char prime[100000000]
Definition: primes_up_to_billion.cpp:10
std::pow
T pow(T... args)
isprime
bool isprime[1000006]
Definition: prime_factorization.cpp:13