mirror of
https://github.com/TheAlgorithms/C-Plus-Plus.git
synced 2026-02-11 14:36:25 +08:00
Merge pull request #17 from kvedala/numerical_methods
Numerical methods - LU and QR decomposition of matrices
This commit is contained in:
Krishna Vedala
@@ -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
|
||||
|
||||
126
numerical_methods/lu_decompose.cpp
Normal file
126
numerical_methods/lu_decompose.cpp
Normal file
@@ -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;
|
||||
}
|
||||
210
numerical_methods/qr_decompose.h
Normal file
210
numerical_methods/qr_decompose.h
Normal file
@@ -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_
|
||||
58
numerical_methods/qr_decomposition.cpp
Normal file
58
numerical_methods/qr_decomposition.cpp
Normal file
@@ -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;
|
||||
}
|
||||
269
numerical_methods/qr_eigen_values.cpp
Normal file
269
numerical_methods/qr_eigen_values.cpp
Normal file
@@ -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;
|
||||
}
|
||||
Reference in New Issue
Block a user