qr_eigen_values.cpp File Reference

Compute real eigen values and eigen vectors of a symmetric matrix using QR decomposition method. More...

#include <cassert>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include "./qr_decompose.h"

Include dependency graph for qr_eigen_values.cpp:

Namespaces
namespace	qr_algorithm
	Functions to compute QR decomposition of any rectangular matrix.

Macros
#define	LIMS   9

Functions
void	create_matrix (std::valarray< std::valarray< double > > *A)
void	mat_mul (const std::valarray< std::valarray< double > > &A, const std::valarray< std::valarray< double > > &B, std::valarray< std::valarray< double > > *OUT)
std::valarray< double >	qr_algorithm::eigen_values (std::valarray< std::valarray< double > > *A, bool print_intermediates=false)
void	test1 ()
void	test2 ()
int	main (int argc, char **argv)

Detailed Description

Compute real eigen values and eigen vectors of a symmetric matrix using QR decomposition method.

Author

Krishna Vedala

Macro Definition Documentation

LIMS

#define LIMS   9

limit of range of matrix values

Function Documentation

create_matrix()

void create_matrix

(

std::valarray< std::valarray< double > > *

A

)

create a symmetric square matrix of given size with random elements. A symmetric square matrix will always have real eigen values.

Parameters

[out]

A

matrix to create (must be pre-allocated in memory)

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

Here is the call graph for this function:

main()

int main	(	int	argc,
		char **	argv )

main function

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

Here is the call graph for this function:

mat_mul()

void mat_mul	(	const std::valarray< std::valarray< double > > &	A,
		const std::valarray< std::valarray< double > > &	B,
		std::valarray< std::valarray< double > > *	OUT )

Perform multiplication of two matrices.

• R2 must be equal to C1
• Resultant matrix size should be R1xC2

  Parameters

  [in]	A	first matrix to multiply
  [in]	B	second matrix to multiply
  [out]	OUT	output matrix (must be pre-allocated)

  Returns

  pointer to resultant matrix

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

test1()

void test1

(

)

test function to compute eigen values of a 2x2 matrix

\[\begin{bmatrix} 5 & 7\\ 7 & 11 \end{bmatrix}\]

which are approximately, {15.56158, 0.384227}

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

Here is the call graph for this function:

test2()

void test2

(

)

test function to compute eigen values of a 2x2 matrix

\[\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}\]

which are approximately, {9.27648, -9.26948, 2.0181, -1.03516, -5.98994}

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

Here is the call graph for this function: