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

T	datatype of input matrix - int, unsigned int, double, etc

Parameters

A	input square matrix

Returns

determinant of matrix A

                                          {
    matrix<double> L(A.size(), std::valarray<double>(A.size()));
    matrix<double> U(A.size(), std::valarray<double>(A.size()));
 
    if (lu_decomposition(A, &L, &U) < 0)
        return 0;
 
    double result = 1.f;
    for (size_t i = 0; i < A.size(); i++) {
        result *= L[i][i] * U[i][i];
    }
    return result;
}

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]	A	matrix to decompose
[out]	L	output L matrix
[out]	U	output U matrix

Returns

0 if no errors

negative if error occurred

                                                                               {
    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;
}

Here is the call graph for this function: