// QR ITERATION #############################################################
/*
QR iteration for finding eigenvalues
LANGUAGE: C++0x
PRECISION: double
AUTHOR: Hendrik Wernecke <wernecke@itp.uni-frankfurt.de>
VERSION DATE: October 19, 2016
LICENCE: free to share, distribute and modify
INCLUDES: standard libraries (and possible bugs)
PRELIMINARIES
works so far only with square matrices
only valid for matrices with real eigenvalues
DESCRIPTION
By means of Householder decomposition a m*m matrix A is decomposed
A = Q*R
where Q is a m*m orthogonal matrix and R is a m*m upper triangular matrix.
Using the iteration
A_n-1 = Q_n*R_n
and
A_n = R_n*Q_n
with initially
A0 = A
one can -- under certain circumstances -- converged in the limit
lim_{n\to\infty}A_n = A'
towards an upper triangular matrix A' which has the same eigenvalues
than the original matrix A.
*/
//###########################################################################
/* ### INCLUDES ### */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
/* ### FUNCTIONS forward definitions ### */
// definitions below MAIN
double** newMatrix(int m, int n);
void matrixCopy(double** mat0, double** mat1, int m, int n);
void UnitMatrix(double** mat, int m, int n);
void printMatrix(double** mat, int m, int n);
void printArray(double* arr, int m);
double norm(double* arr, int m);
void getColumn(double** mat, int k, double* arr, int m, int n);
void multiplyMatrices(double** mat1, double** mat2, double** erg, int m, int n);
void updateAMatrix(double** Ah, double** Qk, int k, int m, int n);
void transposeMatrix(double** mat, double** tmat, int m, int n);
bool qrHouseholderDecomp(double** A0, double** Q, double** R, int m, int n, double tol);
/* ### MAIN ### */
int main(int argLength, char* argValues[])
{
/* SETTINGS */
int maxitnum = 1000; // maximal number of iterations
double tol = 1e-4; // tolerance for QR decomposition
// and accuracy of eigenvalues
/*
// example 1) dimension 3*3
const int m = 3, n = 3; // dimensions of matrix (n!=m not yet implemented)
double A[m][n] = { // matrix to be decomposed
{12, -51, 4},
{ 6, 167, -68},
{-4, 24, -41}
};
*/
/*
// example 2) dimension 4*4
const int m = 4, n = 4; // dimensions of matrix (n!=m not yet implemented)
double A[m][n] = { // matrix to be decomposed
{ 4, 1, 3, -2},
{ 1, -2, 4, 1},
{ 3, 4, 1, 2},
{-2, 1, 2, 3}
};
*/
// example 3) symmetric matrix
int size = 8;
if (argLength==2)
size = atof(argValues[1]);
const int m=size, n=size;
double A[m][n];
for (int i=0; i<m; i++) {
for (int j=i; j<n; j++) {
A[i][j] = rand()%100;
A[j][i] = A[i][j];
}
}
/* COMPUTATIONS */
// map A to real 2D array
double** A0 = newMatrix(m, n);
for (int i=0; i<m; i++) {
for (int j=0; j<n; j++) {
A0[i][j] = A[i][j];
}
}
// allocate matrices for storage of results
double** Aprev = newMatrix(m, n);
matrixCopy(A0, Aprev, m, n);
double** Anext = newMatrix(m, n);
double** Q = newMatrix(m, n);
double** R = newMatrix(m, n);
// some auxiliar variables
int itnum = 0;
bool check = false;
// get started with the iteration
while ( itnum<maxitnum )
{
//printf("A_%i = \n", itnum); printMatrix(Aprev, m, n);
if (not qrHouseholderDecomp(Aprev, Q, R, m, n, tol)) {
printf("QR decomposition failed!\n");
return 1;
}
multiplyMatrices(R, Q, Anext, m, n);
// check if converged (comparing diagonal entries)
double diff = 0;
for (int i=0; i<m && i<n; i++) {
diff += fabs(Aprev[i][i]-Anext[i][i]);
}
//printf("diff = %f\n", diff);
if ( diff<=tol ) { check = true; break; }
matrixCopy(Anext, Aprev, m, n);
itnum++;
}
/* RESULTS */
// announce final result
if (check) {
printf("eigenvalues converged within %i steps of iteration\n", itnum);
printf("original matrix A_0 = \n"); printMatrix(A0, m, n);
printf("iterated matrix A_%i = \n", itnum); printMatrix(Anext, m, n);
printf("eigenvalues are: ");
for (int i=0; i<m; i++) { printf(" %f, ", Anext[i][i]); } printf("\n");
} else {
printf("iteration did not converge in %i steps\nresulting in\n", itnum);
printf("A_%i\n", itnum);
printMatrix(Anext, m, n);
}
/* CLEANING UP */
delete [] A0;
delete [] Aprev;
delete [] Anext;
delete [] Q;
delete [] R;
return 0;
}
// END of MAIN
/* ### FUNCTIONS DEFINITIONS ### */
bool qrHouseholderDecomp(double** A0, double** Q, double** R, const int m, const int n, double tol)
// performs a QR decomposition base on Householder reflection
// see <https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections>
// takes the matrix <A0> to be decomposed and the matrices <Q>, <R> (resulting composites)
// all matrices are <m>*<n> matrices
// for selfconsistent check A0=Q*R the tolerance <tol> is encountered
{
// auxiliar matrix H
double** H = newMatrix(m, n);
UnitMatrix(H, m, n);
// create local copy
double** Ak = newMatrix(m, n);
matrixCopy(A0, Ak, m, n);
for (int k = 0; k<m-1 && k<n ; k++) {
//printf("before k = %i, Ah = \n", k); printMatrix(Ak, m, n);
// get first column
double x[m];
getColumn(Ak, k, x, m, n);
//printf("x = \n"); printArray(x, m);
// compute alpha as norm of first column
double alpha = norm(x, m);
//printf("alpha = %f\n", alpha);
double e[m];
for (int i=0; i<m; i++) { e[i] = 0; } e[k] = 1; // k-th unit vector
double u[m];
for (int i=0; i<m; i++) { u[i] = x[i] - alpha * e[i]; }
double unorm = norm(u, m);
//printf("u = \n"); printArray(u, m);
double v[m];
for (int i=0; i<m; i++) { v[i] = u[i]/unorm; }
//printf("v = \n"); printArray(v, m);
// assemble l-th Householder matrix Qk = 1-2 v*v^T
double** Qk = newMatrix(m, n);
for (int i=0; i<m; i++) {
for (int j=0; j<n; j++) {
if (i==j) { // diagonal
Qk[i][j] = 1 - 2*v[i]*v[j];
} else { // non-diagonal
Qk[i][j] = - 2*v[i]*v[j];
}
}
}
//printf("Qk = \n"); printMatrix(Qk, m, n);
// updating Ah matrix
updateAMatrix(Ak, Qk, k, m, n);
//printf("A_k = \n"); printMatrix(Ak, m, n);
// store product of Q_k (H=Q_n*Q_n-1*...*Q_1)
double** Hh = newMatrix(m, n);
matrixCopy(H, Hh, m, n);
multiplyMatrices(Qk, H, Hh, m, n);
H = newMatrix(m, n);
matrixCopy(Hh, H, m, n);
//printf("H = \n"); printMatrix(H, m, n);
delete [] Hh;
delete [] Qk;
}
// Q = (Hn*...*H1)^T = H^T
UnitMatrix(Q, m, n);
transposeMatrix(H, Q, m, n);
//printf("Q = \n"); printMatrix(Q, m, n);
// R = Q^T*A = H*A
UnitMatrix(R, m, n);
multiplyMatrices(H, A0, R, m, n);
//printf("R = \n"); printMatrix(R, m, n);
// check that A = Q*R
double** Ah = newMatrix(m, n);
multiplyMatrices(Q, R, Ah, m, n);
//printf("A = \n"); printMatrix(Ah, m, n);
// check deviation
double dev = 0;
for (int i=0; i<m; i++) {
for (int j=0; j<n; j++) {
dev += fabs(Ah[i][j]-A0[i][j]);
}
}
if (dev<=tol) {return true;}
else {return false;}
}
double** newMatrix(int m, int n)
// returns pointer to new <m>*<n> matrix
{
// allocate new matrix
double** mat = (double**) malloc(m*sizeof(double));
for (int i=0; i<m; i++) { mat[i] = (double*) malloc(n*sizeof(double)); }
return mat;
}
void matrixCopy(double** mat0, double** mat1, int m, int n)
// copy content of matrix <mat0> to matrix <mat1>, both <m>*<n>
{
// copy entries
for (int i=0; i<m; i++) {
for (int j=0; j<n; j++) {
mat1[i][j] = mat0[i][j];
}
}
return;
}
void UnitMatrix(double** mat, int m, int n)
// fills the <m>*<n> matrix <mat> with a unity matrix
{
// write entries
for (int i=0; i<m; i++) {
for (int j=0; j<n; j++) {
if (i==j) { mat[i][j] = 1; }
else { mat[i][j] = 0; }
}
}
return;
}
void printMatrix(double** mat, int m, int n)
// prints the matrix <mat> with <m> rows and <n> columns
{
printf("\n");
for (int i=0; i<m; i++) {
for (int j=0; j<n; j++) {
printf(" %8.3f", mat[i][j]);
}
printf("\n");
}
printf("\n");
return;
}
void printArray(double* arr, int m)
// prints the array <arr> with length <m>
{
printf("\n");
for (int i=0; i<m; i++) {
printf(" %8.3f", arr[i]);
}
printf("\n");
return;
}
double norm(double* arr, int m)
// computes the norm of an array <arr> with length <m>
{
double h = 0;
for (int i=0; i<m; i++) { h += arr[i]*arr[i]; }
return sqrt(h);
}
void getColumn(double** mat, int k, double* arr, int m, int n)
// extracts the <k>-th column of a <m>*<n> matrix <mat> and saves it to array <arr>
{
if ( k<0 or k>=n ) { arr = NULL; return; }
for (int i=0; i<m; i++) { arr[i] = mat[i][k]; }
return;
}
void multiplyMatrices(double** mat1, double** mat2, double** erg, int m, int n)
// multiplies the matrices <mat1> and <mat2>
{
for (int i=0; i<m; i++) { for (int j=0; j<n; j++) { erg[i][j] = 0; } }
// compute product
for (int i=0; i<m; i++) {
for (int j=0; j<n; j++) {
for (int k=0; k<n; k++) {
erg[i][j] += mat1[i][k]*mat2[k][j];
}
}
}
return;
}
void updateAMatrix(double** Ah, double** Qk, int k, int m, int n)
// updateds the new matrix <Ah> for next Householder step
// takes the Householder matrix <Qk>, both <m>*<n> matrices
// and the Householder iteration step number <k>
{
// copy Qk*A
double** QA = newMatrix(m, n);
multiplyMatrices(Qk, Ah, QA, m, n);
// create dummy
double** H = newMatrix(m, n);
/* assemble new matrix */
// k+1 zero columns
for (int i=0; i<k+1; i++) { for (int j=0; j<n; j++) { H[i][j] = 0; } }
// k+1 zero rows
for (int i=0; i<m; i++) { for (int j=0; j<k+1; j++) { H[i][j] = 0; } }
// k+1 diagonal elements set to one
for (int i=0; i<k+1; i++) { H[i][i] = 1; }
// rest of the matrix filled with values from QA
for (int i=k+1; i<m; i++) { for (int j=k+1; j<n; j++) { H[i][j] = QA[i][j]; } }
// copy result to Ah
for (int i=0; i<m; i++) { for (int j=0; j<n; j++) { Ah[i][j] = H[i][j]; } }
// clean up
delete [] QA;
delete [] H;
return;
}
void transposeMatrix(double** mat, double** tmat, int m, int n)
// transposes the <m>*<n> matrix <mat>
// and stores the transpose in the <n>*<m> matrix <tmat>
{
// transpose
for (int j=0; j<n; j++) {
for (int i=0; i<m; i++) {
tmat[j][i] = mat[i][j];
}
}
return;
}
// ### END OF FILE ###
