qr_eigen_values.c File Reference

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

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "qr_decompose.h"

Include dependency graph for qr_eigen_values.c:

Macros
#define	LIMS   9
	limit of range of matrix values
#define	EPSILON   1e-10
	accuracy tolerance limit

Functions
void	create_matrix (double **A, int N)
	create a square matrix of given size with random elements
double **	mat_mul (double **A, double **B, double **OUT, int R1, int C1, int R2, int C2)
	Perform multiplication of two matrices.
double	eigen_values (double **A, double *eigen_vals, int mat_size, char debug_print)
	Compute eigen values using iterative shifted QR decomposition algorithm as follows:
void	test1 ()
	test function to compute eigen values of a 2x2 matrix
void	test2 ()
	test function to compute eigen values of a 2x2 matrix
int	main (int argc, char **argv)
	main function

Detailed Description

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

Author

Krishna Vedala

Function Documentation

create_matrix()

void create_matrix	(	double **	A,
		int	N
	)

create a square matrix of given size with random elements

Parameters

[out]	A	matrix to create (must be pre-allocated in memory)
[in]	N	matrix size

{
    int i, j, tmp, lim2 = LIMS >> 1;
 
#ifdef _OPENMP
#pragma omp for
#endif
    for (i = 0; i < N; i++)
    {
        A[i][i] = (rand() % LIMS) - lim2;
        for (j = i + 1; j < N; j++)
        {
            tmp = (rand() % LIMS) - lim2;
            A[i][j] = tmp;
            A[j][i] = tmp;
        }
    }
}

eigen_values()

double eigen_values	(	double **	A,
		double *	eigen_vals,
		int	mat_size,
		char	debug_print
	)

Compute eigen values using iterative shifted QR decomposition algorithm as follows:

1. Use last diagonal element of A as eigen value approximation \(c\)
2. Shift diagonals of matrix \(A' = A - cI\)
3. Decompose matrix \(A'=QR\)
4. Compute next approximation \(A'_1 = RQ \)
5. Shift diagonals back \(A_1 = A'_1 + cI\)
6. Termination condition check: last element below diagonal is almost 0
   1. If not 0, go back to step 1 with the new approximation \(A_1\)
   2. If 0, continue to step 7
7. Save last known \(c\) as the eigen value.
8. Are all eigen values found?
   1. If not, remove last row and column of \(A_1\) and go back to step 1.
   2. If yes, stop.

Note

The matrix \(A\) gets modified

Parameters

[in,out]	A	matrix to compute eigen values for
[out]	eigen_vals	resultant vector containing computed eigen values
[in]	mat_size	matrix size
[in]	debug_print	1 to print intermediate Q & R matrices, 0 for not to

Returns

time for computation in seconds

{
    if (!eigen_vals)
    {
        perror("Output eigen value vector cannot be NULL!");
        return -1;
    }
    double **R = (double **)malloc(sizeof(double *) * mat_size);
    double **Q = (double **)malloc(sizeof(double *) * mat_size);
    if (!Q || !R)
    {
        perror("Unable to allocate memory for Q & R!");
        if (Q)
        {
            free(Q);
        }
        if (R)
        {
            free(R);
        }
        return -1;
    }
 
    /* allocate dynamic memory for matrices */
    for (int i = 0; i < mat_size; i++)
    {
        R[i] = (double *)malloc(sizeof(double) * mat_size);
        Q[i] = (double *)malloc(sizeof(double) * mat_size);
        if (!Q[i] || !R[i])
        {
            perror("Unable to allocate memory for Q & R.");
            for (; i >= 0; i--)
            {
                free(R[i]);
                free(Q[i]);
            }
            free(Q);
            free(R);
            return -1;
        }
    }
 
    if (debug_print)
    {
        print_matrix(A, mat_size, mat_size);
    }
 
    int rows = mat_size, columns = mat_size;
    int counter = 0, num_eigs = rows - 1;
    double last_eig = 0;
 
    clock_t t1 = clock();
    while (num_eigs > 0) /* continue till all eigen values are found */
    {
        /* iterate with QR decomposition */
        while (fabs(A[num_eigs][num_eigs - 1]) > EPSILON)
        {
            last_eig = A[num_eigs][num_eigs];
            for (int i = 0; i < rows; i++) A[i][i] -= last_eig; /* A - cI */
            qr_decompose(A, Q, R, rows, columns);
 
            if (debug_print)
            {
                print_matrix(A, rows, columns);
                print_matrix(Q, rows, columns);
                print_matrix(R, columns, columns);
                printf("-------------------- %d ---------------------\n",
                       ++counter);
            }
 
            mat_mul(R, Q, A, columns, columns, rows, columns);
            for (int i = 0; i < rows; i++) A[i][i] += last_eig; /* A + cI */
        }
 
        /* store the converged eigen value */
        eigen_vals[num_eigs] = last_eig;
 
        if (debug_print)
        {
            printf("========================\n");
            printf("Eigen value: % g,\n", last_eig);
            printf("========================\n");
        }
 
        num_eigs--;
        rows--;
        columns--;
    }
    eigen_vals[0] = A[0][0];
    double dtime = (double)(clock() - t1) / CLOCKS_PER_SEC;
 
    if (debug_print)
    {
        print_matrix(R, mat_size, mat_size);
        print_matrix(Q, mat_size, mat_size);
    }
 
    /* cleanup dynamic memory */
    for (int i = 0; i < mat_size; i++)
    {
        free(R[i]);
        free(Q[i]);
    }
    free(R);
    free(Q);
 
    return dtime;
}

Here is the call graph for this function:

main()

int main	(	int	argc,
		char **	argv
	)

main function

{
    srand(time(NULL));
 
    int mat_size = 5;
    if (argc == 2)
    {
        mat_size = atoi(argv[1]);
    }
    else
    {  // if invalid input argument is given run tests
        test1();
        test2();
        printf("Usage: ./qr_eigen_values [mat_size]\n");
        return 0;
    }
 
    if (mat_size < 2)
    {
        fprintf(stderr, "Matrix size should be > 2\n");
        return -1;
    }
 
    int i;
 
    double **A = (double **)malloc(sizeof(double *) * mat_size);
    /* number of eigen values = matrix size */
    double *eigen_vals = (double *)malloc(sizeof(double) * mat_size);
    if (!eigen_vals)
    {
        perror("Unable to allocate memory for eigen values!");
        free(A);
        return -1;
    }
    for (i = 0; i < mat_size; i++)
    {
        A[i] = (double *)malloc(sizeof(double) * mat_size);
        eigen_vals[i] = 0.f;
    }
 
    /* create a random matrix */
    create_matrix(A, mat_size);
 
    print_matrix(A, mat_size, mat_size);
 
    double dtime = eigen_values(A, eigen_vals, mat_size, 0);
    printf("Eigen vals: ");
    for (i = 0; i < mat_size; i++) printf("% 9.4g\t", eigen_vals[i]);
    printf("\nTime taken to compute: % .4g sec\n", dtime);
 
    for (int i = 0; i < mat_size; i++) free(A[i]);
    free(A);
    free(eigen_vals);
    return 0;
}

Here is the call graph for this function:

mat_mul()

double ** mat_mul	(	double **	A,
		double **	B,
		double **	OUT,
		int	R1,
		int	C1,
		int	R2,
		int	C2
	)

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)
  [in]	R1	number of rows of first matrix
  [in]	C1	number of columns of first matrix
  [in]	R2	number of rows of second matrix
  [in]	C2	number of columns of second matrix

  Returns

  pointer to resultant matrix

{
    if (C1 != R2)
    {
        perror("Matrix dimensions mismatch!");
        return OUT;
    }
 
    int i;
#ifdef _OPENMP
#pragma omp for
#endif
    for (i = 0; i < R1; i++)
    {
        for (int j = 0; j < C2; j++)
        {
            OUT[i][j] = 0.f;
            for (int k = 0; k < C1; k++) OUT[i][j] += A[i][k] * B[k][j];
        }
    }
    return OUT;
}

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}

{
    int mat_size = 2;
    double X[][2] = {{5, 7}, {7, 11}};
    double y[] = {15.56158, 0.384227};  // corresponding y-values
    double eig_vals[2] = {0, 0};
 
    // The following steps are to convert a "double[][]" to "double **"
    double **A = (double **)malloc(mat_size * sizeof(double *));
    for (int i = 0; i < mat_size; i++) A[i] = X[i];
 
    printf("------- Test 1 -------\n");
 
    double dtime = eigen_values(A, eig_vals, mat_size, 0);
 
    for (int i = 0; i < mat_size; i++)
    {
        printf("%d/5 Checking for %.3g --> ", i + 1, y[i]);
        char result = 0;
        for (int j = 0; j < mat_size && !result; j++)
        {
            if (fabs(y[i] - eig_vals[j]) < 0.1)
            {
                result = 1;
                printf("(%.3g) ", eig_vals[j]);
            }
        }
 
        // ensure that i^th expected eigen value was computed
        assert(result != 0);
        printf("found\n");
    }
    printf("Test 1 Passed in %.3g sec\n\n", dtime);
    free(A);
}

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}

{
    int mat_size = 5;
    double X[][5] = {{-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
    double eig_vals[5];
 
    // The following steps are to convert a "double[][]" to "double **"
    double **A = (double **)malloc(mat_size * sizeof(double *));
    for (int i = 0; i < mat_size; i++) A[i] = X[i];
 
    printf("------- Test 2 -------\n");
 
    double dtime = eigen_values(A, eig_vals, mat_size, 0);
 
    for (int i = 0; i < mat_size; i++)
    {
        printf("%d/5 Checking for %.3g --> ", i + 1, y[i]);
        char result = 0;
        for (int j = 0; j < mat_size && !result; j++)
        {
            if (fabs(y[i] - eig_vals[j]) < 0.1)
            {
                result = 1;
                printf("(%.3g) ", eig_vals[j]);
            }
        }
 
        // ensure that i^th expected eigen value was computed
        assert(result != 0);
        printf("found\n");
    }
    printf("Test 2 Passed in %.3g sec\n\n", dtime);
    free(A);
}

Here is the call graph for this function: