OpenSoure_Repos/TheAlgorithms-C-Plus-Plus
1
0
Fork 0
mirror of https://hub.njuu.cf/TheAlgorithms/C-Plus-Plus.git synced 2023-10-11 13:05:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

split lu_decomposition to a header file and templated the function

Browse Source
This commit is contained in:
Krishna Vedala 2020-06-25 13:38:11 -04:00
parent 8cd25f4f8a
commit 2c61414a83
No known key found for this signature in database
GPG Key ID: BA19ACF8FC8792F7
2 changed files with 73 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
numerical_methods/lu_decompose.cpp
View File

@ -7,73 +7,15 @@
#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;
}
#include "./lu_decomposition.h"
/**
* operator to print a matrix
*/
template <typename T>
std::ostream &operator<<(std::ostream &out,
std::vector<std::vector<T>> const &v) {
std::vector<std::valarray<T>> const &v) {
const int width = 10;
const char separator = ' ';
@ -99,14 +41,14 @@ int main(int argc, char **argv) {
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
std::vector<std::valarray<double>> A(mat_size,
std::valarray<double>(mat_size));
std::vector<std::valarray<double>> L(
mat_size, std::valarray<double>(mat_size)); // output
std::vector<std::valarray<double>> U(
mat_size, std::valarray<double>(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);

65
numerical_methods/lu_decomposition.h Normal file
View File

@ -0,0 +1,65 @@
#pragma once
#include <iostream>
#include <valarray>
#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
*/
template <typename T>
int lu_decomposition(const std::vector<std::valarray<T>> &A,
std::vector<std::valarray<double>> *L,
std::vector<std::valarray<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;
}