2
\$\begingroup\$

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!

\$\endgroup\$
2
  • \$\begingroup\$ I suggest you use some software to properly format the code in the future, the indentation is inconsistent and that makes reading the code more difficult. \$\endgroup\$
    – pacmaninbw
    Jan 12, 2023 at 0:51
  • \$\begingroup\$ Indents really matter, they affect how we reason about code, even if the compiler ignores them. I found this code quite illegible, until I ran it through astyle. \$\endgroup\$
    – J_H
    Jan 12, 2023 at 2:15

1 Answer 1

2
\$\begingroup\$

Use the Proper Compiler Flags

There are several issues here that could be caught using the proper compiler flags.

Using C++ from Visual Studio 2022

1>C:\CodeReview\Givens\Givens\main.cpp(97,38): warning C4101: 'sol4': unreferenced local variable
1>C:\CodeReview\Givens\Givens\main.cpp(97,26): warning C4101: 'sol2': unreferenced local variable
1>C:\CodeReview\Givens\Givens\main.cpp(97,20): warning C4101: 'sol1': unreferenced local variable
1>C:\CodeReview\Givens\Givens\main.cpp(97,32): warning C4101: 'sol3': unreferenced local variable
1>C:\CodeReview\Givens\Givens\main.cpp(186,12): warning C4101: 'max': unreferenced local variable
1>C:\CodeReview\Givens\Givens\main.cpp(185,12): warning C4101: 'min': unreferenced local variable
1>C:\CodeReview\Givens\Givens\main.cpp(250): error C4700: uninitialized local variable 'extreme2' used
1>C:\CodeReview\Givens\Givens\main.cpp(250): error C4700: uninitialized local variable 'extreme1' used
1>C:\CodeReview\Givens\Givens\main.cpp(299): error C4700: uninitialized local variable 't' used
1>C:\CodeReview\Givens\Givens\main.cpp(306): error C4700: uninitialized local variable 'z' used

Using gcc on Ubuntu with the following compiler flags: -Wall -Wextra -pedantic -Werror

main.cpp: In function ‘void buildU(double**, double**, int, int, int, double)’:
main.cpp:93:12: error: unused variable ‘x’ [-Werror=unused-variable]
   93 |     double x = sin(t) / cos(t);
      |            ^
main.cpp:97:20: error: unused variable ‘sol1’ [-Werror=unused-variable]
   97 |     double s1, s2, sol1, sol2, sol3, sol4;
      |                    ^~~~
main.cpp:97:26: error: unused variable ‘sol2’ [-Werror=unused-variable]
   97 |     double s1, s2, sol1, sol2, sol3, sol4;
      |                          ^~~~
main.cpp:97:32: error: unused variable ‘sol3’ [-Werror=unused-variable]
   97 |     double s1, s2, sol1, sol2, sol3, sol4;
      |                                ^~~~
main.cpp:97:38: error: unused variable ‘sol4’ [-Werror=unused-variable]
   97 |     double s1, s2, sol1, sol2, sol3, sol4;
      |                                      ^~~~
main.cpp: In function ‘void traspU(double**, double**, int, int, int)’:
main.cpp:115:65: error: unused parameter ‘j’ [-Werror=unused-parameter]
  115 | void traspU(double** matrixUt, double** matrixU, int order, int j, int i)
      |                                                             ~~~~^
main.cpp:115:72: error: unused parameter ‘i’ [-Werror=unused-parameter]
  115 | void traspU(double** matrixUt, double** matrixU, int order, int j, int i)
      |                                                                    ~~~~^
main.cpp: In function ‘int main()’:
main.cpp:185:12: error: unused variable ‘min’ [-Werror=unused-variable]
  185 |     double min;
      |            ^~~
main.cpp:186:12: error: unused variable ‘max’ [-Werror=unused-variable]
  186 |     double max;
      |            ^~~
main.cpp:250:27: error: ‘extreme1’ may be used uninitialized [-Werror=maybe-uninitialized]
  250 |                 val_random(order, extreme1, extreme2, i);
      |                 ~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
main.cpp:250:27: error: ‘extreme2’ may be used uninitialized [-Werror=maybe-uninitialized]
main.cpp:299:31: error: ‘t’ may be used uninitialized [-Werror=maybe-uninitialized]
  299 |                         buildU(matrixU, matrixA, order, i, j, t);
      |                         ~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
main.cpp:306:26: error: ‘z’ may be used uninitialized [-Werror=maybe-uninitialized]
  306 |                     prodd(matrixA2, matrixA1, matrixUt, matrixU, matrixA, order, j, i, z);
      |                     ~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Please note that in this section of code j is uninitialized in the call to diagonalize() because the call is outside the for loop where j is defined. That line corresponds to the warning/error message for line 311 above.

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

Avoid using namespace std;

If you are coding professionally you probably should get out of the habit of using the using namespace std; statement. The code will more clearly define where cout and other identifiers are coming from (std::cin, std::cout). As you start using namespaces in your code it is better to identify where each function comes from because there may be function name collisions from different namespaces. The identifiercout you may override within your own classes, and you may override the operator << in your own classes as well. This stack overflow question discusses this in more detail.

Declare the Variables as Needed

In the original version of C back in the 1970s and 1980s variables had to be declared at the top of the function. That is no longer the case, and a recommended programming practice to declare the variable as needed. In C++ the language doesn't provide a default initialization of the variable so variables should be initialized as part of the declaration. For readability and maintainability each variable should be declared and initialized on its own line.

DRY Code

There is a programming principle called the Don't Repeat Yourself Principle sometimes referred to as DRY code. If you find yourself repeating the same code multiple times it is better to encapsulate it in a function. If it is possible to loop through the code that can reduce repetition as well.

Here is some of the code repetition in main() that can be fixed with a function:

        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];
        }

In the previous code a function could be created that allocates the memory for the matrix.

Complexity

The main() function is 140 lines of code and comments. A general good practice in programming is to keep every function to one screen in the IDE or editor so that anyone that needs to maintain the code can easily understand the function. One screen in an IDE or editor is generally less than 60 lines, so the main() function is 1.5 screens too large.

The function main() is too complex (does too much). As programs grow in size the use of main() should be limited to calling functions that parse the command line, calling functions that set up for processing, calling functions that execute the desired function of the program, and calling functions to clean up after the main portion of the program.

There is also a programming principle called the Single Responsibility Principle that applies here. The Single Responsibility Principle states:

that every module, class, or function should have responsibility over a single part of the functionality provided by the software, and that responsibility should be entirely encapsulated by that module, class or function.

In addition to the function suggested in the DRY Code section above create a function for the outer while (true) loop, then break that function into smaller functions with only one task per function.

Memory Leaks

There are 6 matrices with memory allocated, none of these matrices are every de-allocated. If this was part of a larger longer running program there would definitely be memory leaks and the process could possibly run out of memory at some point.

\$\endgroup\$
5
  • \$\begingroup\$ Thank you for the information. I am trying to fix what I can. I already fixed the std part. Now I'm tryng to do something about the memory leaks. I tried to control them with address Sanitizer but when launghing the command it does not return any message or error nor anything. Is there something I can do about? If not I'm proceding with "delete" (de allocating) \$\endgroup\$
    – Henry
    Jan 18, 2023 at 11:58
  • \$\begingroup\$ P.s. I'm trying to compile with visual studio too, but I do not obtain any of your error/warning messages, using gcc... \$\endgroup\$
    – Henry
    Jan 18, 2023 at 12:11
  • \$\begingroup\$ I use Visual Studio Professional, I use warning level 3 (W3). You can change the warning level in the Project Properties. Project Properties -> C/C++ -> Warning Level. \$\endgroup\$
    – pacmaninbw
    Jan 18, 2023 at 14:28
  • \$\begingroup\$ I have visual studio Code. Is it the same, ish? I cannot find "project properties"! \$\endgroup\$
    – Henry
    Jan 18, 2023 at 14:38
  • \$\begingroup\$ @Henry Not the same, your version is free and works on other operating systems, my version cost $800 USD a year, I haven't tried it on Linux (I use VI on Linux). If you are a student you might be able to get my version for much less. \$\endgroup\$
    – pacmaninbw
    Jan 18, 2023 at 14:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.