diff --git a/DIRECTORY.md b/DIRECTORY.md
index 31093341a..45143d3a1 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -129,8 +129,12 @@
   * [Durand Kerner Roots](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/durand_kerner_roots.cpp)
   * [False Position](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/false_position.cpp)
   * [Gaussian Elimination](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/gaussian_elimination.cpp)
+  * [Lu Decompose](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/lu_decompose.cpp)
   * [Newton Raphson Method](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/newton_raphson_method.cpp)
   * [Ordinary Least Squares Regressor](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/ordinary_least_squares_regressor.cpp)
+  * [Qr Decompose](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_decompose.h)
+  * [Qr Decomposition](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_decomposition.cpp)
+  * [Qr Eigen Values](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/qr_eigen_values.cpp)
   * [Successive Approximation](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/numerical_methods/successive_approximation.cpp)
 
 ## Operations On Datastructures
diff --git a/numerical_methods/lu_decompose.cpp b/numerical_methods/lu_decompose.cpp
new file mode 100644
index 000000000..a0a2d00ab
--- /dev/null
+++ b/numerical_methods/lu_decompose.cpp
@@ -0,0 +1,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;
+}
diff --git a/numerical_methods/qr_decompose.h b/numerical_methods/qr_decompose.h
new file mode 100644
index 000000000..c9c369144
--- /dev/null
+++ b/numerical_methods/qr_decompose.h
@@ -0,0 +1,210 @@
+/**
+ * @file
+ * \brief Library functions to compute [QR
+ * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given
+ * matrix.
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+
+#ifndef NUMERICAL_METHODS_QR_DECOMPOSE_H_
+#define NUMERICAL_METHODS_QR_DECOMPOSE_H_
+
+#include <cmath>
+#include <cstdlib>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <numeric>
+#include <valarray>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+/** \namespace qr_algorithm
+ * \brief Functions to compute [QR
+ * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of any
+ * rectangular matrix
+ */
+namespace qr_algorithm {
+/**
+ * operator to print a matrix
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &out,
+                         std::valarray<std::valarray<T>> const &v) {
+    const int width = 12;
+    const char separator = ' ';
+
+    out.precision(4);
+    for (size_t row = 0; row < v.size(); row++) {
+        for (size_t col = 0; col < v[row].size(); col++)
+            out << std::right << std::setw(width) << std::setfill(separator)
+                << v[row][col];
+        out << std::endl;
+    }
+
+    return out;
+}
+
+/**
+ * operator to print a vector
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &out, std::valarray<T> const &v) {
+    const int width = 10;
+    const char separator = ' ';
+
+    out.precision(4);
+    for (size_t row = 0; row < v.size(); row++) {
+        out << std::right << std::setw(width) << std::setfill(separator)
+            << v[row];
+    }
+
+    return out;
+}
+
+/**
+ * Compute dot product of two vectors of equal lengths
+ *
+ * If \f$\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\f$ and
+ * \f$\vec{b}=\left[b_0,b_1,b_1,...,b_L\right]\f$ then
+ * \f$\vec{a}\cdot\vec{b}=\displaystyle\sum_{i=0}^L a_i\times b_i\f$
+ *
+ * \returns \f$\vec{a}\cdot\vec{b}\f$
+ */
+template <typename T>
+inline double vector_dot(const std::valarray<T> &a, const std::valarray<T> &b) {
+    return (a * b).sum();
+    // could also use following
+    // return std::inner_product(std::begin(a), std::end(a), std::begin(b),
+    // 0.f);
+}
+
+/**
+ * Compute magnitude of vector.
+ *
+ * If \f$\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\f$ then
+ * \f$\left|\vec{a}\right|=\sqrt{\displaystyle\sum_{i=0}^L a_i^2}\f$
+ *
+ * \returns \f$\left|\vec{a}\right|\f$
+ */
+template <typename T>
+inline double vector_mag(const std::valarray<T> &a) {
+    double dot = vector_dot(a, a);
+    return std::sqrt(dot);
+}
+
+/**
+ * Compute projection of vector \f$\vec{a}\f$ on \f$\vec{b}\f$ defined as
+ * \f[\text{proj}_\vec{b}\vec{a}=\frac{\vec{a}\cdot\vec{b}}{\left|\vec{b}\right|^2}\vec{b}\f]
+ *
+ * \returns NULL if error, otherwise pointer to output
+ */
+template <typename T>
+std::valarray<T> vector_proj(const std::valarray<T> &a,
+                             const std::valarray<T> &b) {
+    double num = vector_dot(a, b);
+    double deno = vector_dot(b, b);
+
+    /*! check for division by zero using machine epsilon */
+    if (deno <= std::numeric_limits<double>::epsilon()) {
+        std::cerr << "[" << __func__ << "] Possible division by zero\n";
+        return a;  // return vector a back
+    }
+
+    double scalar = num / deno;
+
+    return b * scalar;
+}
+
+/**
+ * Decompose matrix \f$A\f$ using [Gram-Schmidt
+ *process](https://en.wikipedia.org/wiki/QR_decomposition).
+ *
+ * \f{eqnarray*}{
+ * \text{given that}\quad A &=&
+ *\left[\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_{N-1},\right]\\
+ * \text{where}\quad\mathbf{a}_i &=&
+ * \left[a_{0i},a_{1i},a_{2i},\ldots,a_{(M-1)i}\right]^T\quad\ldots\mbox{(column
+ * vectors)}\\
+ * \text{then}\quad\mathbf{u}_i &=& \mathbf{a}_i
+ *-\sum_{j=0}^{i-1}\text{proj}_{\mathbf{u}_j}\mathbf{a}_i\\
+ * \mathbf{e}_i &=&\frac{\mathbf{u}_i}{\left|\mathbf{u}_i\right|}\\
+ * Q &=& \begin{bmatrix}\mathbf{e}_0 & \mathbf{e}_1 & \mathbf{e}_2 & \dots &
+ * \mathbf{e}_{N-1}\end{bmatrix}\\
+ * R &=& \begin{bmatrix}\langle\mathbf{e}_0\,,\mathbf{a}_0\rangle &
+ * \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle &
+ * \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots \\
+ *                  0 & \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle &
+ * \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\
+ *                  0 & 0 & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle &
+ * \dots\\ \vdots & \vdots & \vdots & \ddots
+ *      \end{bmatrix}\\
+ * \f}
+ */
+template <typename T>
+void qr_decompose(
+    const std::valarray<std::valarray<T>> &A, /**< input matrix to decompose */
+    std::valarray<std::valarray<T>> *Q,       /**< output decomposed matrix */
+    std::valarray<std::valarray<T>> *R        /**< output decomposed matrix */
+) {
+    std::size_t ROWS = A.size();        // number of rows of A
+    std::size_t COLUMNS = A[0].size();  // number of columns of A
+    std::valarray<T> col_vector(ROWS);
+    std::valarray<T> col_vector2(ROWS);
+    std::valarray<T> tmp_vector(ROWS);
+
+    for (int i = 0; i < COLUMNS; i++) {
+        /* for each column => R is a square matrix of NxN */
+        int j;
+        R[0][i] = 0.; /* make R upper triangular */
+
+        /* get corresponding Q vector */
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (j = 0; j < ROWS; j++) {
+            tmp_vector[j] = A[j][i]; /* accumulator for uk */
+            col_vector[j] = A[j][i];
+        }
+        for (j = 0; j < i; j++) {
+            for (int k = 0; k < ROWS; k++) {
+                col_vector2[k] = Q[0][k][j];
+            }
+            col_vector2 = vector_proj(col_vector, col_vector2);
+            tmp_vector -= col_vector2;
+        }
+
+        double mag = vector_mag(tmp_vector);
+
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (j = 0; j < ROWS; j++) Q[0][j][i] = tmp_vector[j] / mag;
+
+            /* compute upper triangular values of R */
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (int kk = 0; kk < ROWS; kk++) {
+            col_vector[kk] = Q[0][kk][i];
+        }
+
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (int k = i; k < COLUMNS; k++) {
+            for (int kk = 0; kk < ROWS; kk++) {
+                col_vector2[kk] = A[kk][k];
+            }
+            R[0][i][k] = (col_vector * col_vector2).sum();
+        }
+    }
+}
+}  // namespace qr_algorithm
+
+#endif  // NUMERICAL_METHODS_QR_DECOMPOSE_H_
diff --git a/numerical_methods/qr_decomposition.cpp b/numerical_methods/qr_decomposition.cpp
new file mode 100644
index 000000000..237a5c946
--- /dev/null
+++ b/numerical_methods/qr_decomposition.cpp
@@ -0,0 +1,58 @@
+/**
+ * @file
+ * \brief Program to compute the [QR
+ * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given
+ * matrix.
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+
+#include <array>
+#include <cmath>
+#include <cstdlib>
+#include <ctime>
+#include <iostream>
+
+#include "./qr_decompose.h"
+
+using qr_algorithm::qr_decompose;
+using qr_algorithm::operator<<;
+
+/**
+ * main function
+ */
+int main(void) {
+    unsigned int ROWS, COLUMNS;
+
+    std::cout << "Enter the number of rows and columns: ";
+    std::cin >> ROWS >> COLUMNS;
+
+    std::cout << "Enter matrix elements row-wise:\n";
+
+    std::valarray<std::valarray<double>> A(ROWS);
+    std::valarray<std::valarray<double>> Q(ROWS);
+    std::valarray<std::valarray<double>> R(COLUMNS);
+    for (int i = 0; i < std::max(ROWS, COLUMNS); i++) {
+        if (i < ROWS) {
+            A[i] = std::valarray<double>(COLUMNS);
+            Q[i] = std::valarray<double>(COLUMNS);
+        }
+        if (i < COLUMNS) {
+            R[i] = std::valarray<double>(COLUMNS);
+        }
+    }
+
+    for (int i = 0; i < ROWS; i++)
+        for (int j = 0; j < COLUMNS; j++) std::cin >> A[i][j];
+
+    std::cout << A << "\n";
+
+    clock_t t1 = clock();
+    qr_decompose(A, &Q, &R);
+    double dtime = static_cast<double>(clock() - t1) / CLOCKS_PER_SEC;
+
+    std::cout << Q << "\n";
+    std::cout << R << "\n";
+    std::cout << "Time taken to compute: " << dtime << " sec\n ";
+
+    return 0;
+}
diff --git a/numerical_methods/qr_eigen_values.cpp b/numerical_methods/qr_eigen_values.cpp
new file mode 100644
index 000000000..991453775
--- /dev/null
+++ b/numerical_methods/qr_eigen_values.cpp
@@ -0,0 +1,269 @@
+/**
+ * @file
+ * \brief Compute real eigen values and eigen vectors of a symmetric matrix
+ * using [QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition)
+ * method.
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+#include <ctime>
+#include <iostream>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#include "./qr_decompose.h"
+using qr_algorithm::operator<<;
+
+#define LIMS 9 /**< limit of range of matrix values */
+
+/**
+ * create a symmetric square matrix of given size with random elements. A
+ * symmetric square matrix will *always* have real eigen values.
+ *
+ * \param[out] A matrix to create (must be pre-allocated in memory)
+ */
+void create_matrix(std::valarray<std::valarray<double>> *A) {
+    int i, j, tmp, lim2 = LIMS >> 1;
+    int N = A->size();
+
+#ifdef _OPENMP
+#pragma omp for
+#endif
+    for (i = 0; i < N; i++) {
+        A[0][i][i] = (std::rand() % LIMS) - lim2;
+        for (j = i + 1; j < N; j++) {
+            tmp = (std::rand() % LIMS) - lim2;
+            A[0][i][j] = tmp;  // summetrically distribute random values
+            A[0][j][i] = tmp;
+        }
+    }
+}
+
+/**
+ * Perform multiplication of two matrices.
+ * * R2 must be equal to C1
+ * * Resultant matrix size should be R1xC2
+ * \param[in] A first matrix to multiply
+ * \param[in] B second matrix to multiply
+ * \param[out] OUT output matrix (must be pre-allocated)
+ * \returns pointer to resultant matrix
+ */
+void mat_mul(const std::valarray<std::valarray<double>> &A,
+             const std::valarray<std::valarray<double>> &B,
+             std::valarray<std::valarray<double>> *OUT) {
+    int R1 = A.size();
+    int C1 = A[0].size();
+    int R2 = B.size();
+    int C2 = B[0].size();
+    if (C1 != R2) {
+        perror("Matrix dimensions mismatch!");
+        return;
+    }
+
+    for (int i = 0; i < R1; i++) {
+        for (int j = 0; j < C2; j++) {
+            OUT[0][i][j] = 0.f;
+            for (int k = 0; k < C1; k++) {
+                OUT[0][i][j] += A[i][k] * B[k][j];
+            }
+        }
+    }
+}
+
+namespace qr_algorithm {
+/** Compute eigen values
+ * \param[in,out] A matric to compute eigen values for \note This matrix gets
+ * modified
+ * \param[in] print_intermediates (optional) whether to print intermediate A, Q
+ * and R matrices (default = `false`)
+ */
+std::valarray<double> eigen_values(std::valarray<std::valarray<double>> *A,
+                                   bool print_intermediates = false) {
+    int rows = A->size();
+    int columns = rows;
+    int counter = 0, num_eigs = rows - 1;
+    double last_eig = 0;
+
+    std::valarray<std::valarray<double>> Q(rows);
+    std::valarray<std::valarray<double>> R(columns);
+
+    /* number of eigen values = matrix size */
+    std::valarray<double> eigen_vals(rows);
+    for (int i = 0; i < rows; i++) {
+        Q[i] = std::valarray<double>(columns);
+        R[i] = std::valarray<double>(columns);
+    }
+
+    /* continue till all eigen values are found */
+    while (num_eigs > 0) {
+        /* iterate with QR decomposition */
+        while (std::abs(A[0][num_eigs][num_eigs - 1]) >
+               std::numeric_limits<double>::epsilon()) {
+            // initial approximation = last diagonal element
+            last_eig = A[0][num_eigs][num_eigs];
+            for (int i = 0; i < rows; i++) {
+                A[0][i][i] -= last_eig; /* A - cI */
+            }
+
+            qr_decompose(*A, &Q, &R);
+
+            if (print_intermediates) {
+                std::cout << *A << "\n";
+                std::cout << Q << "\n";
+                std::cout << R << "\n";
+                printf("-------------------- %d ---------------------\n",
+                       ++counter);
+            }
+
+            // new approximation A' = R * Q
+            mat_mul(R, Q, A);
+
+            for (int i = 0; i < rows; i++) {
+                A[0][i][i] += last_eig; /* A + cI */
+            }
+        }
+
+        /* store the converged eigen value */
+        eigen_vals[num_eigs] = last_eig;
+        // A[0][num_eigs][num_eigs];
+        if (print_intermediates) {
+            std::cout << "========================\n";
+            std::cout << "Eigen value: " << last_eig << ",\n";
+            std::cout << "========================\n";
+        }
+
+        num_eigs--;
+        rows--;
+        columns--;
+    }
+    eigen_vals[0] = A[0][0][0];
+
+    if (print_intermediates) {
+        std::cout << Q << "\n";
+        std::cout << R << "\n";
+    }
+
+    return eigen_vals;
+}
+
+}  // namespace qr_algorithm
+
+/**
+ * test function to compute eigen values of a 2x2 matrix
+ * \f[\begin{bmatrix}
+ * 5 & 7\\
+ * 7 & 11
+ * \end{bmatrix}\f]
+ * which are approximately, {15.56158, 0.384227}
+ */
+void test1() {
+    std::valarray<std::valarray<double>> X = {{5, 7}, {7, 11}};
+    double y[] = {15.56158, 0.384227};  // corresponding y-values
+
+    std::cout << "------- Test 1 -------" << std::endl;
+    std::valarray<double> eig_vals = qr_algorithm::eigen_values(&X);
+
+    for (int i = 0; i < 2; i++) {
+        std::cout << i + 1 << "/2 Checking for " << y[i] << " --> ";
+        bool result = false;
+        for (int j = 0; j < 2 && !result; j++) {
+            if (std::abs(y[i] - eig_vals[j]) < 0.1) {
+                result = true;
+                std::cout << "(" << eig_vals[j] << ") ";
+            }
+        }
+        assert(result);  // ensure that i^th expected eigen value was computed
+        std::cout << "found\n";
+    }
+    std::cout << "Test 1 Passed\n\n";
+}
+
+/**
+ * test function to compute eigen values of a 2x2 matrix
+ * \f[\begin{bmatrix}
+ * -4& 4& 2& 0& -3\\
+ * 4& -4& 4& -3& -1\\
+ * 2& 4& 4& 3& -3\\
+ * 0& -3& 3& -1&-1\\
+ * -3& -1& -3& -3& 0
+ * \end{bmatrix}\f]
+ * which are approximately, {9.27648, -9.26948, 2.0181, -1.03516, -5.98994}
+ */
+void test2() {
+    std::valarray<std::valarray<double>> X = {{-4, 4, 2, 0, -3},
+                                              {4, -4, 4, -3, -1},
+                                              {2, 4, 4, 3, -3},
+                                              {0, -3, 3, -1, -3},
+                                              {-3, -1, -3, -3, 0}};
+    double y[] = {9.27648, -9.26948, 2.0181, -1.03516,
+                  -5.98994};  // corresponding y-values
+
+    std::cout << "------- Test 2 -------" << std::endl;
+    std::valarray<double> eig_vals = qr_algorithm::eigen_values(&X);
+
+    std::cout << X << "\n"
+              << "Eigen values: " << eig_vals << "\n";
+
+    for (int i = 0; i < 5; i++) {
+        std::cout << i + 1 << "/5 Checking for " << y[i] << " --> ";
+        bool result = false;
+        for (int j = 0; j < 5 && !result; j++) {
+            if (std::abs(y[i] - eig_vals[j]) < 0.1) {
+                result = true;
+                std::cout << "(" << eig_vals[j] << ") ";
+            }
+        }
+        assert(result);  // ensure that i^th expected eigen value was computed
+        std::cout << "found\n";
+    }
+    std::cout << "Test 2 Passed\n\n";
+}
+
+/**
+ * main function
+ */
+int main(int argc, char **argv) {
+    int mat_size = 5;
+    if (argc == 2) {
+        mat_size = atoi(argv[1]);
+    } else {  // if invalid input argument is given run tests
+        test1();
+        test2();
+        std::cout << "Usage: ./qr_eigen_values [mat_size]\n";
+        return 0;
+    }
+
+    if (mat_size < 2) {
+        fprintf(stderr, "Matrix size should be > 2\n");
+        return -1;
+    }
+
+    // initialize random number generator
+    std::srand(std::time(nullptr));
+
+    int i, rows = mat_size, columns = mat_size;
+
+    std::valarray<std::valarray<double>> A(rows);
+
+    for (int i = 0; i < rows; i++) {
+        A[i] = std::valarray<double>(columns);
+    }
+
+    /* create a random matrix */
+    create_matrix(&A);
+
+    std::cout << A << "\n";
+
+    clock_t t1 = clock();
+    std::valarray<double> eigen_vals = qr_algorithm::eigen_values(&A);
+    double dtime = static_cast<double>(clock() - t1) / CLOCKS_PER_SEC;
+
+    std::cout << "Eigen vals: ";
+    for (i = 0; i < mat_size; i++) std::cout << eigen_vals[i] << "\t";
+    std::cout << "\nTime taken to compute: " << dtime << " sec\n";
+
+    return 0;
+}