Algorithms_in_C++  1.0.0
Set of algorithms implemented in C++.
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

 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]Amatrix to create (must be pre-allocated in memory)
28  {
29  int i, j, tmp, lim2 = LIMS >> 1;
30  int N = A->size();
31 
32 #ifdef _OPENMP
33 #pragma omp for
34 #endif
35  for (i = 0; i < N; i++) {
36  A[0][i][i] = (std::rand() % LIMS) - lim2;
37  for (j = i + 1; j < N; j++) {
38  tmp = (std::rand() % LIMS) - lim2;
39  A[0][i][j] = tmp; // summetrically distribute random values
40  A[0][j][i] = tmp;
41  }
42  }
43 }
Here is the call graph for this function:

◆ main()

int main ( int  argc,
char **  argv 
)

main function

243  {
244  int mat_size = 5;
245  if (argc == 2) {
246  mat_size = atoi(argv[1]);
247  } else { // if invalid input argument is given run tests
248  test1();
249  test2();
250  std::cout << "Usage: ./qr_eigen_values [mat_size]\n";
251  return 0;
252  }
253 
254  if (mat_size < 2) {
255  fprintf(stderr, "Matrix size should be > 2\n");
256  return -1;
257  }
258 
259  // initialize random number generator
260  std::srand(std::time(nullptr));
261 
262  int i, rows = mat_size, columns = mat_size;
263 
265 
266  for (int i = 0; i < rows; i++) {
267  A[i] = std::valarray<double>(columns);
268  }
269 
270  /* create a random matrix */
271  create_matrix(&A);
272 
273  std::cout << A << "\n";
274 
275  clock_t t1 = clock();
277  double dtime = static_cast<double>(clock() - t1) / CLOCKS_PER_SEC;
278 
279  std::cout << "Eigen vals: ";
280  for (i = 0; i < mat_size; i++) std::cout << eigen_vals[i] << "\t";
281  std::cout << "\nTime taken to compute: " << dtime << " sec\n";
282 
283  return 0;
284 }
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]Afirst matrix to multiply
    [in]Bsecond matrix to multiply
    [out]OUToutput matrix (must be pre-allocated)
    Returns
    pointer to resultant matrix
56  {
57  int R1 = A.size();
58  int C1 = A[0].size();
59  int R2 = B.size();
60  int C2 = B[0].size();
61  if (C1 != R2) {
62  perror("Matrix dimensions mismatch!");
63  return;
64  }
65 
66  for (int i = 0; i < R1; i++) {
67  for (int j = 0; j < C2; j++) {
68  OUT[0][i][j] = 0.f;
69  for (int k = 0; k < C1; k++) {
70  OUT[0][i][j] += A[i][k] * B[k][j];
71  }
72  }
73  }
74 }

◆ 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}

177  {
178  std::valarray<std::valarray<double>> X = {{5, 7}, {7, 11}};
179  double y[] = {15.56158, 0.384227}; // corresponding y-values
180 
181  std::cout << "------- Test 1 -------" << std::endl;
183 
184  for (int i = 0; i < 2; i++) {
185  std::cout << i + 1 << "/2 Checking for " << y[i] << " --> ";
186  bool result = false;
187  for (int j = 0; j < 2 && !result; j++) {
188  if (std::abs(y[i] - eig_vals[j]) < 0.1) {
189  result = true;
190  std::cout << "(" << eig_vals[j] << ") ";
191  }
192  }
193  assert(result); // ensure that i^th expected eigen value was computed
194  std::cout << "found\n";
195  }
196  std::cout << "Test 1 Passed\n\n";
197 }
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}

210  {
211  std::valarray<std::valarray<double>> X = {{-4, 4, 2, 0, -3},
212  {4, -4, 4, -3, -1},
213  {2, 4, 4, 3, -3},
214  {0, -3, 3, -1, -3},
215  {-3, -1, -3, -3, 0}};
216  double y[] = {9.27648, -9.26948, 2.0181, -1.03516,
217  -5.98994}; // corresponding y-values
218 
219  std::cout << "------- Test 2 -------" << std::endl;
221 
222  std::cout << X << "\n"
223  << "Eigen values: " << eig_vals << "\n";
224 
225  for (int i = 0; i < 5; i++) {
226  std::cout << i + 1 << "/5 Checking for " << y[i] << " --> ";
227  bool result = false;
228  for (int j = 0; j < 5 && !result; j++) {
229  if (std::abs(y[i] - eig_vals[j]) < 0.1) {
230  result = true;
231  std::cout << "(" << eig_vals[j] << ") ";
232  }
233  }
234  assert(result); // ensure that i^th expected eigen value was computed
235  std::cout << "found\n";
236  }
237  std::cout << "Test 2 Passed\n\n";
238 }
Here is the call graph for this function:
std::srand
T srand(T... args)
mat_size
ll mat_size
Definition: matrix_exponentiation.cpp:45
create_matrix
void create_matrix(std::valarray< std::valarray< double >> *A)
Definition: qr_eigen_values.cpp:28
qr_algorithm::eigen_values
std::valarray< double > eigen_values(std::valarray< std::valarray< double >> *A, bool print_intermediates=false)
Definition: qr_eigen_values.cpp:98
std::clock
T clock(T... args)
std::fprintf
T fprintf(T... args)
std::cout
std::atoi
T atoi(T... args)
std::perror
T perror(T... args)
std::valarray
STL class.
test2
void test2()
Definition: qr_eigen_values.cpp:210
std::rand
T rand(T... args)
test1
void test1()
Definition: qr_eigen_values.cpp:177
std::endl
T endl(T... args)
std::time
T time(T... args)
LIMS
#define LIMS
Definition: qr_eigen_values.cpp:20