Algorithms_in_C++  1.0.0
Set of algorithms implemented in C++.
lu_decomposition.h
Go to the documentation of this file.
1 /**
2  * @file lu_decomposition.h
3  * @author [Krishna Vedala](https://github.com/kvedala)
4  * @brief Functions associated with [LU
5  * Decomposition](https://en.wikipedia.org/wiki/LU_decomposition)
6  * of a square matrix.
7  */
8 #pragma once
9 
10 #include <iostream>
11 #include <valarray>
12 #include <vector>
13 #ifdef _OPENMP
14 #include <omp.h>
15 #endif
16 
17 /** Define matrix type as a `std::vector` of `std::valarray` */
18 template <typename T>
20 
21 /** Perform LU decomposition on matrix
22  * \param[in] A matrix to decompose
23  * \param[out] L output L matrix
24  * \param[out] U output U matrix
25  * \returns 0 if no errors
26  * \returns negative if error occurred
27  */
28 template <typename T>
30  int row, col, j;
31  int mat_size = A.size();
32 
33  if (mat_size != A[0].size()) {
34  // check matrix is a square matrix
35  std::cerr << "Not a square matrix!\n";
36  return -1;
37  }
38 
39  // regularize each row
40  for (row = 0; row < mat_size; row++) {
41  // Upper triangular matrix
42 #ifdef _OPENMP
43 #pragma omp for
44 #endif
45  for (col = row; col < mat_size; col++) {
46  // Summation of L[i,j] * U[j,k]
47  double lu_sum = 0.;
48  for (j = 0; j < row; j++) {
49  lu_sum += L[0][row][j] * U[0][j][col];
50  }
51 
52  // Evaluate U[i,k]
53  U[0][row][col] = A[row][col] - lu_sum;
54  }
55 
56  // Lower triangular matrix
57 #ifdef _OPENMP
58 #pragma omp for
59 #endif
60  for (col = row; col < mat_size; col++) {
61  if (row == col) {
62  L[0][row][col] = 1.;
63  continue;
64  }
65 
66  // Summation of L[i,j] * U[j,k]
67  double lu_sum = 0.;
68  for (j = 0; j < row; j++) {
69  lu_sum += L[0][col][j] * U[0][j][row];
70  }
71 
72  // Evaluate U[i,k]
73  L[0][col][row] = (A[col][row] - lu_sum) / U[0][row][row];
74  }
75  }
76 
77  return 0;
78 }
79 
80 /**
81  * Compute determinant of an NxN square matrix using LU decomposition.
82  * Using LU decomposition, the determinant is given by the product of diagonal
83  * elements of matrices L and U.
84  *
85  * @tparam T datatype of input matrix - int, unsigned int, double, etc
86  * @param A input square matrix
87  * @return determinant of matrix A
88  */
89 template <typename T>
90 double determinant_lu(const matrix<T> &A) {
93 
94  if (lu_decomposition(A, &L, &U) < 0)
95  return 0;
96 
97  double result = 1.f;
98  for (size_t i = 0; i < A.size(); i++) {
99  result *= L[i][i] * U[i][i];
100  }
101  return result;
102 }
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< double >
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
machine_learning::sum
T sum(const std::vector< std::valarray< T >> &A)
Definition: vector_ops.hpp:232
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