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` */
18template <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 */
28template <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 */
89template <typename T>
90double 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}
double determinant_lu(const matrix< T > &A)
Definition lu_decomposition.h:90
int lu_decomposition(const matrix< T > &A, matrix< double > *L, matrix< double > *U)
Definition lu_decomposition.h:29
ll mat_size
Definition matrix_exponentiation.cpp:45
T size(T... args)