So I tried to create a code which aims to diagonalise symmetric matrices via Givens algorithm. The program only works for symmetric matrices, for we know a symmetric matrix is always diagonalisable, and it takes values which will be set as columns. Due to the matrix being symmetrix, the user will type all the values for the first column, then $n-1$ values for the second column, $n-2£ values for the third column and so on.
For example if the matrix has order three, then, say, $1 3 4$ is the first column, then we may choose $9 2$ as the two bottom elements for the second column (the first one is $3$, directly understood by the program), and the final element say it's $10$, for the third column (first two will be $4$ and $2$).
This will serve as an example too, so let's diagonalise
1 3 4
3 9 2
4 2 4
I actually did not want to write here the whole code, but I really need some help.
The code is the following:
#include <random>
#include <iostream>
#include <sstream>
#include <fstream>
using namespace std;
// calculate the norm
long double calculate_norm( double **matrixA, int order)
{
long double norm =0;
for (int i=0; i!=order; ++i)
for(int j=0; j!=order; ++j)
norm += matrixA[i][j]*matrixA[i][j];
return sqrt(norm);
}
//write random values
void val_random (int m, double extreme1, double extreme2, int i)
{
cout<<"put the extremes within which you want to generate the values"<<'\n';
cin>>extreme1;
cin>>extreme2;
double min= extreme1;
double max=extreme2;
ofstream os("RANDOM");
while(true)
{
uniform_real_distribution<double> unif(min,max);
default_random_engine re;
for(i=0; i< (m*(m+1))/2 ; ++i)
{
double a_random_double = unif(re);
os<<a_random_double<<'\n';
}
break;
}
}
//write the output matrices
void stamp (double **matrixA1, int m)
{
for(int i=0; i<m; ++i) //print matrix U transpose in output
{
for(int j=0; j<m; ++j)
{
cout<< matrixA1 [i][j];
cout<<"\t";
}
cout<<endl;
}
cout << "________________________________\n";
}
//building of U matrix
void buildU ( double **matrixU, double **matrixA, int order, int i, int j, double t)
{
for(int j=0; j<order; ++j) // identity matrix
{
for(int i=0; i<order; ++i)
{
if (j==i)
matrixU[j] [i]=1;
else
matrixU[j] [i]=0;
cout<<endl;
}
}
double b;
b = ( matrixA [i][i] - matrixA [j][j] )/(2*matrixA [i][j]);
//cout<<b<<'\n';
double x= sin(t)/cos(t);
//cout<< x<<'\n';
double s1, s2, sol1, sol2, sol3, sol4;
s1= b+sqrt(1+b*b);
s2= b-sqrt(1+b*b);
if(fabs(s1) < fabs(s2))
t= s1;
else
t = s2;
matrixU [i][i] = matrixU [j][j] = 1/(sqrt(1+t*t)); //k row l column
matrixU [i][j] = matrixU [j][i] = t/(sqrt(1+t*t));
matrixU [j][i] *= -1.0;
cout<<"matrix U"<<'\n';
stamp(matrixU, order);
}
// U transpose matrix
void traspU (double **matrixUt, double **matrixU, int order, int j, int i)
{
for(int j=0; j<order; ++j)
{
for(int i=0; i<order; ++i)
{
matrixUt[j] [i]= matrixU[i] [j];
}
}
}
//matrix product
void prodd(double **matrixA2, double **matrixA1, double **matrixUt, double **matrixU, double **matrixA, int order, int j, int i, int z)
{
for(i=0; i<order; ++i)
for(j=0; j<order; ++j)
{
matrixA1 [i][j] = 0;
for(z=0; z<order; ++z)
matrixA1 [i][j] += matrixUt[i][z] * matrixA[z][j];
}
cout<<"matrixA1"<<'\n';
stamp(matrixA1, order);
for(i=0; i<order; ++i)
for(j=0; j<order; ++j)
{
matrixA2 [i][j] = 0;
for(z=0; z<order; ++z)
matrixA2 [i][j] += matrixA1[i][z] * matrixU [z][j];
}
for(int j=0; j<order; ++j)
for(int i=0; i<order; ++i)
{
matrixA [i][j] = matrixA2 [i][j];
}
}
bool diagonalize(double **matrixA2, int order, int i, int j, long double norm)
{
for(j=0; j<order; ++j)
{
for(i=j+1; i<order; ++i)
{
if ( fabs(matrixA2[i][j]) >= norm*10e-9)
{
cout<<""<<'\n';
return(false);
}
}
}
return (true);
}
int main()
{
int order;
int i;
int j;
double t;
int z;
double extreme1,extreme2;
double min;
double max;
long double norm;
ifstream input_stream_da_disco; //declaration type+object
while (true)
{
//A
cout<<"This program diagonalises symmetric matrices by using the Givens method. Choose the order of the symmetric matrix by typing a positive integer "<<'\n';
cin>>order;
if (!order) break;
double **matrixA = new double *[order];
for (i=0; i<order; ++i)
{
matrixA[i]=new double [order];
}
double **matrixU = new double *[order];
for (i=0; i<order; ++i)
{
matrixU[i]=new double [order];
}
double **matrixUt = new double *[order];
for (i=0; i<order; ++i)
{
matrixUt[i]=new double [order];
}
double **matrixA1 = new double *[order];
for (i=0; i<order; ++i)
{
matrixA1[i]=new double [order];
}
double **matrixA2 = new double *[order];
for (i=0; i<order; ++i)
{
matrixA2[i]=new double [order];
}
double **matrixAf = new double *[order];
for (i=0; i<order; ++i)
{
matrixAf[i]=new double [order];
}
cout<<"Choose how you want to insert the values of the matrix: choose (1) for a keyboard input or (2) if you want to insert values from a document "<<'\n';
int k;
cin>>k;
if(k==2)
{
cout<<"Choose (3) if you want to include a personal document or (4) to generate random values " <<'\n';
int o;
cin>>o;
if(o==3)
cout<<"Personal document"<<'\n';
else if(o==4)
{
//cout<<"random values";
val_random (order, extreme1, extreme2, i);
input_stream_da_disco . open("RANDOM");
for(j=0; j<order; ++j)//filling matrix A
{
for(i=j; i<order; ++i)
{
input_stream_da_disco>> matrixA [j][i];
matrixA[i][j] = matrixA[j][i];
}
}
stamp(matrixA, order);
}
if(!input_stream_da_disco)
{
cout<<"The document is not valid \n";
}
}
else if(k==1)
{
cout<<"Insert the values of the inferior triangular matrix, which will be inserted as columns "<<'\n';
for(j=0; j<order; ++j)//filling of matrix A
{
for(i=j; i<order; ++i)
{
cin>> matrixA [j][i];
matrixA[i][j] = matrixA[j][i];
}
}
stamp(matrixA, order);
}
else break;
norm = calculate_norm( matrixA, order);
//cout<<"norm\n" <<norm;
cout<<"Now A will be diagonalized through the following matrices U: "<<'\n';
while (true) //overwrite the identity matrix with sines and cosines
{
for(int i=0; i!=order-1; ++i) //cosine1
for(int j=i+1; j!=order; ++j) //cosine2
{
if (matrixA [i][j]==0) continue;
else
buildU( matrixU, matrixA, order, i, j, t);
traspU (matrixUt, matrixU, order, i, j); // transpose of U
cout<<"the transpose of U is " <<'\n';
stamp(matrixUt, order);
prodd(matrixA2, matrixA1, matrixUt, matrixU, matrixA, order, j, i, z);
}
cout<<"the final matrix \n";
stamp(matrixA2, order);
if(diagonalize( matrixA2, order, i, j, norm)) break;
}
}
return 0;}
I run with W. Mathematica the matrix above and I can say the output is quite correct (the eigenvalues are those ones). This worksin general.
My doubt and questions are basically if I can improve the program in some way. I am interested in the computation of matrices defined by the user, so I'm not really into "improvind the input document part", to say.
Is there something I can manage to get better precision? Or anyway: do you believe this program is sufficient, as a first try? Something strange, wrong? (it does compile without errors or warnings).
Thank you so much!