Algorithms_in_C++
1.0.0
Set of algorithms implemented in C++.
Loading...
Searching...
No Matches
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>
19
using
matrix
=
std::vector<std::valarray<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>
29
int
lu_decomposition
(
const
matrix<T>
&A,
matrix<double>
*L,
matrix<double>
*U) {
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) {
91
matrix<double>
L(A.
size
(),
std::valarray<double>
(A.
size
()));
92
matrix<double>
U(A.
size
(),
std::valarray<double>
(A.
size
()));
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
}
std::cerr
determinant_lu
double determinant_lu(const matrix< T > &A)
Definition
lu_decomposition.h:90
lu_decomposition
int lu_decomposition(const matrix< T > &A, matrix< double > *L, matrix< double > *U)
Definition
lu_decomposition.h:29
mat_size
ll mat_size
Definition
matrix_exponentiation.cpp:45
std::vector::size
T size(T... args)
std::valarray
std::vector
numerical_methods
lu_decomposition.h
Generated by
1.12.0