Algorithms_in_C++  1.0.0
Set of algorithms implemented in C++.
lu_decomposition.h File Reference

Functions associated with LU Decomposition of a square matrix. More...

#include <iostream>
#include <valarray>
#include <vector>
Include dependency graph for lu_decomposition.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Typedefs

template<typename T >
using matrix = std::vector< std::valarray< T > >
 

Functions

template<typename T >
int lu_decomposition (const matrix< T > &A, matrix< double > *L, matrix< double > *U)
 
template<typename T >
double determinant_lu (const matrix< T > &A)
 

Detailed Description

Functions associated with LU Decomposition of a square matrix.

Author
Krishna Vedala

Typedef Documentation

◆ matrix

template<typename T >
using matrix = std::vector<std::valarray<T> >

Define matrix type as a std::vector of std::valarray

Function Documentation

◆ determinant_lu()

template<typename T >
double determinant_lu ( const matrix< T > &  A)

Compute determinant of an NxN square matrix using LU decomposition. Using LU decomposition, the determinant is given by the product of diagonal elements of matrices L and U.

Template Parameters
Tdatatype of input matrix - int, unsigned int, double, etc
Parameters
Ainput square matrix
Returns
determinant of matrix A
90  {
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 }
Here is the call graph for this function:

◆ lu_decomposition()

template<typename T >
int lu_decomposition ( const matrix< T > &  A,
matrix< double > *  L,
matrix< double > *  U 
)

Perform LU decomposition on matrix

Parameters
[in]Amatrix to decompose
[out]Loutput L matrix
[out]Uoutput U matrix
Returns
0 if no errors
negative if error occurred
29  {
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 }
Here is the call graph for this function:
mat_size
ll mat_size
Definition: matrix_exponentiation.cpp:45
std::vector
STL class.
std::vector::size
T size(T... args)
std::cerr
std::valarray< double >
lu_decomposition
int lu_decomposition(const matrix< T > &A, matrix< double > *L, matrix< double > *U)
Definition: lu_decomposition.h:29