forked from TheAlgorithms/C-Plus-Plus
-
Notifications
You must be signed in to change notification settings - Fork 0
/
lu_decompose.cpp
126 lines (108 loc) · 3.55 KB
/
lu_decompose.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
/**
* \file
* \brief [LU decomposition](https://en.wikipedia.org/wiki/LU_decompositon) of a
* square matrix
* \author [Krishna Vedala](https://github.com/kvedala)
*/
#include <ctime>
#include <iomanip>
#include <iostream>
#include <vector>
#ifdef _OPENMP
#include <omp.h>
#endif
/** 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
*/
int lu_decomposition(const std::vector<std::vector<double>> &A,
std::vector<std::vector<double>> *L,
std::vector<std::vector<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;
}
/**
* operator to print a matrix
*/
template <typename T>
std::ostream &operator<<(std::ostream &out,
std::vector<std::vector<T>> const &v) {
const int width = 10;
const char separator = ' ';
for (size_t row = 0; row < v.size(); row++) {
for (size_t col = 0; col < v[row].size(); col++)
out << std::left << std::setw(width) << std::setfill(separator)
<< v[row][col];
out << std::endl;
}
return out;
}
/** Main function */
int main(int argc, char **argv) {
int mat_size = 3; // default matrix size
const int range = 50;
const int range2 = range >> 1;
if (argc == 2)
mat_size = atoi(argv[1]);
std::srand(std::time(NULL)); // random number initializer
/* Create a square matrix with random values */
std::vector<std::vector<double>> A(mat_size);
std::vector<std::vector<double>> L(mat_size); // output
std::vector<std::vector<double>> U(mat_size); // output
for (int i = 0; i < mat_size; i++) {
// calloc so that all valeus are '0' by default
A[i] = std::vector<double>(mat_size);
L[i] = std::vector<double>(mat_size);
U[i] = std::vector<double>(mat_size);
for (int j = 0; j < mat_size; j++)
/* create random values in the limits [-range2, range-1] */
A[i][j] = static_cast<double>(std::rand() % range - range2);
}
std::clock_t start_t = std::clock();
lu_decomposition(A, &L, &U);
std::clock_t end_t = std::clock();
std::cout << "Time taken: "
<< static_cast<double>(end_t - start_t) / CLOCKS_PER_SEC << "\n";
std::cout << "A = \n" << A << "\n";
std::cout << "L = \n" << L << "\n";
std::cout << "U = \n" << U << "\n";
return 0;
}