#include <cmath>
using namespace std;
double** matrix(int, int, int, int);
double* Vector(int, int);
void CleanMatrix(double**, int, int);
void CleanVector(double*, int);
int Jacobi(double**, double**, double[], int);
const int n = 4;
int main() {
double** A;
double** X;
double* D;
int i, j, k;
int check;
double temp;
A = matrix(1, n, 1, n);
X = matrix(1, n, 1, n);
D = Vector(1, n);
/*
A is a real symmetric matrix (lower triangle)
X are the eigenvectors in columns
D is a vector of eigenvalues
n is the size of the problem
*/
A[1][1] = 1;
A[2][1] = 2; A[2][2] = 1;
A[3][1] = 3; A[3][2] = 2; A[3][3] = 1;
A[4][1] = 4; A[4][2] = 3; A[4][3] = 2; A[4][4] = 1;
for (i = 2; i <= n; i++) {
for (j = 1; j <= i - 1; j++) {
A[j][i] = A[i][j];
}
}
Jacobi(A, X, D, n);
cout << "Eigenvalues: ";
for (i = 1; i <= n; i++) {
cout << D[i] << ", ";
}
for (i = 2; i <= n; i++) {
for (j = 1; i <= i-1; i++) {
A[i][j] = A[j][i];
}
}
check = 0e0;
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
temp = -X[i][j]*D[j];
for (k = 1; k <= n; k++) {
temp += A[i][k]*X[k][j];
}
if (check < fabs(temp)) {
check = fabs(temp);
}
}
}
cout << endl;
cout << "Error = " << check << endl;
}
int Jacobi(double **A, double **X, double D[], int n) {
const double epsilon = 1e-34;
int maxIter = 100;
int i, j, k, iter;
double aii, aij, ajj, tang, sine, cosine, amax;
// Setting up eigenvalues in diagonals while zeroing off-diagonals
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
X[i][j] = 0e0;
}
X[i][i] = 1e0;
D[i] = A[i][i];
}
// Iterating
for (maxIter = 1; maxIter <= n; maxIter++) {
amax = 0e0;
for (i = 2; i <= n; i++) {
for (j = 1; j <= (i-1); j++) {
aii = D[i];
ajj = D[j];
aij = fabs(A[i][j]);
if (aij > amax) {
amax = aij;
}
// Performing rotation on the matrix
if (aij > epsilon) {
cosine = 0.5e0*(aii-ajj)/A[i][j];
tang = 1e0/(fabs(cosine)+sqrt((cosine*cosine)+1));
if (cosine < 0e0) {
tang = -tang;
}
cosine = 1e0/(sqrt((tang*tang)+1));
sine = cosine*tang;
for (k = 1; k <= (j-1); k++) {
tang = A[j][k]*cosine - A[i][k]*sine;
A[i][k] = A[i][k]*cosine + A[j][k]*sine;
A[j][k] = tang;
}
// Interchanging
for (k = (j+1); k <= (i-1); k++) {
tang = A[k][j]*cosine - A[i][k]*sine;
A[i][k] = A[i][k]*cosine + A[k][j]*sine;
A[k][j] = tang;
}
for (k = (i+1); k <= n; k++) {
tang = A[k][j]*cosine - A[k][i]*sine;
A[k][i] = A[k][i]*cosine + A[k][j]*sine;
A[k][j] = tang;
}
for (k = 1; k <= n; k++) {
tang = X[k][j]*cosine - A[k][i]*sine;
X[k][i] = X[k][i]*cosine + X[k][j]*sine;
X[k][j] = tang;
}
// Updating eigenvalues
tang = 2e0*sine*cosine*A[i][j];
D[i] = aii*cosine*cosine + ajj*sine*sine + tang;
D[j] = ajj*cosine*cosine + aii*sine*sine - tang;
A[i][j] = 0e0;
}
}
}
if (amax <= epsilon) {
break;
}
}
if (iter > maxIter) {
return 1;
}
return 0;
}
double **matrix(int imin, int imax, int jmin, int jmax) {
int i;
int i_n = imax-imin+1;
int j_n = jmax-jmin+1;
double **m;
m = (double**) malloc((size_t)(i_n*sizeof(double*)));
if (!m) {
printf("Allocation Error (Matrix)");
exit(1);
}
m -= imin;
m[imin] = (double*) malloc((size_t)(i_n*j_n*sizeof(double)));
if (!m[imin]) {
printf("Allocation Error (Matrix2)");
exit(2);
}
m[imin] -= jmin;
for (i = imin+1; i <= imax; i++) {
m[i] = m[i-1]+j_n;
}
return m;
}
double *Vector(int imin, int imax) {
double *v;
v = (double*) malloc((size_t)((imax-imin+1)*sizeof(double)));
if (!v) {
printf("Allocation Error (Vector)");
exit(1);
}
return v - imin;
}
void CleanMatrix(double **m, int imin, int jmin) {
free((void*)(m[imin]+jmin));
free((void*)(m+imin));
}
void CleanVector(double *v, int imin) {
free((void*)(v+imin));
}
