Domanda Cosa non funziona nel programma? Compilazione pulita ma errore di calcolo

[RandomVariable]

Buongiorno,

Sto scrivendo un codice C++ il quale, attraverso l'algoritmo di Givens, diagonalizzi una matrice simmetrica (ovviamente quadrata) di ordine deciso dall'utente.

Non ho errori di compilazione, eppure quando il programma mostra la matrice trasposta qualcosa non va poiché restituisce una matrice di soli zeri. Ecco il codice.

Qualcuno ha qualche idea per favore?
Ecco il codice.


#include <random>
#include <iostream>
#include <sstream>
#include <fstream>
#include<vector>


// calculate the norm
long double calculate_norm( std::vector<std::vector<double>> matrixA, int order)
{
long double norm =0;
for (int i=0; i!=order; ++i)
for(int j=0; j!=order; ++j)
norm += matrixA[j]*matrixA[j];
return sqrt(norm);
}

//write random values
void random_val (int m, double extreme1, double extreme2, int i)
{
std::cout<<"put the extremes within which you want to generate the values"<<'\n';

std::cin>>extreme1;
std::cin>>extreme2;
double min= extreme1;
double max=extreme2;

std:fstream os("RANDOM");

while(true)
{

std::uniform_real_distribution<double> unif(min,max);
std::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 (std::vector<std::vector<double>> matrixA1, int m)
{
for(int i=0; i<m; ++i) //print matrix U transpose in output
{
for(int j=0; j<m; ++j)
{
std::cout<< matrixA1 [j];
std::cout<<"\t";
}
std::cout<<std::endl;
}
std::cout << "________________________________\n";
}


//building of U matrix
void buildU ( std::vector<std::vector<double>> matrixU, std::vector<std::vector<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] =1;
else
matrixU[j] =0;
std::cout<<std::endl;
}
}

double b;
b = ( matrixA - matrixA [j][j] )/(2*matrixA [j]);
//std::cout<<b<<'\n';



double s1, s2;
s1= b+sqrt(1+b*b);
s2= b-sqrt(1+b*b);

if(fabs(s1) < fabs(s2))
t= s1;
else
t = s2;

matrixU = matrixU [j][j] = 1/(sqrt(1+t*t)); //k row l column
matrixU [j] = matrixU [j] = t/(sqrt(1+t*t));
matrixU [j] *= -1.0;
std::cout<<"matrix U"<<'\n';
stamp(matrixU, order);
}

// U transpose matrix
void traspU (std::vector<std::vector<double>> matrixUt, std::vector<std::vector<double>> matrixU, int order, int j, int i)
{
for( j=0; j<order; ++j)
{
for(i=0; i<order; ++i)
{
matrixUt[j] = matrixU [j];
}
}
}

//matrix product
void prodd(std::vector<std::vector<double>> matrixA2, std::vector<std::vector<double>> matrixA1,
std::vector<std::vector<double>> matrixUt, std::vector<std::vector<double>> matrixU, std::vector<std::vector<double>> matrixA,
int order, int j, int i, int z)
{
for(i=0; i<order; ++i)
for(j=0; j<order; ++j)
{
matrixA1 [j] = 0;
for(z=0; z<order; ++z)
matrixA1 [j] += matrixUt[z] * matrixA[z][j];
}

std::cout<<"matrixA1"<<'\n';
stamp(matrixA1, order);

for(i=0; i<order; ++i)
for(j=0; j<order; ++j)
{
matrixA2 [j] = 0;
for(z=0; z<order; ++z)
matrixA2 [j] += matrixA1[z] * matrixU [z][j];
}

for(int j=0; j<order; ++j)
for(int i=0; i<order; ++i)
{
matrixA [j] = matrixA2 [j];
}
}



bool diagonalize(std::vector<std::vector<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[j]) >= norm*10e-9)
{
return(false);
}
}

}
return (true);
}



int main()
{
int order;
int i;
int j;
double t;
double extreme1;
double extreme2;
int z;
long double norm;
std::ifstream input_stream_da_disco; //declaration type+object

while (true)
{
//A
std::cout<<"This program diagonalises symmetric matrices by using the Givens method. Choose the order of the symmetric matrix by typing a positive integer "<<'\n';
std::cin>>order;
if (!order) break;

std::vector<std::vector<double>> matrixA(order, std::vector<double>(order));
std::vector<std::vector<double>> matrixU(order, std::vector<double>(order));
std::vector<std::vector<double>> matrixUt(order, std::vector<double>(order));
std::vector<std::vector<double>> matrixA1(order, std::vector<double>(order));
std::vector<std::vector<double>> matrixA2(order, std::vector<double>(order));
std::vector<std::vector<double>> matrixAf(order, std::vector<double>(order));

std::cout<<"Choose how you want to insert the values of the matrix: choose (1) for a keyboard input or (2) if you want to generate random values "<<'\n';
int k;
std::cin>>k;

if(k==1)
{
std::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)
{
std::cin>> matrixA [j];
matrixA[j] = matrixA[j];
}
}

stamp(matrixA, order);
}

else if(k==2)
{
//std::cout<<"random values";
random_val (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];
matrixA[j] = matrixA[j];
}
}

stamp(matrixA, order);
}

else break;



norm = calculate_norm( matrixA, order);
//std::cout<<"norm\n" <<norm;
std::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 [j]==0) continue;
else
buildU( matrixU, matrixA, order, i, j, t);


traspU (matrixUt, matrixU, order, i, j); // transpose of U
std::cout<<"the transpose of U is " <<'\n';
stamp(matrixUt, order);

prodd(matrixA2, matrixA1, matrixUt, matrixU, matrixA, order, j, i, z);

}
std::cout<<"the final matrix \n";
stamp(matrixA2, order);
if(diagonalize( matrixA2, order, i, j, norm)) break;
}
}
return 0;}

 

[--- Ra ---]

Ciao @RandomVariable, ti consiglio di indentare il codice e di usare l'apposita funzione del forum per postarlo (l'icona con le parentesi angolari), altrimenti diventa difficile/impossibile leggere e capire quello che hai scritto.

 

[DanyDollaro]

Ci sono molte cose che non vanno, per esempio nella funzione calculate_norm da quel che ho capito stai utilizzando l'agoritmo di Frobenius ma lo hai implementato male.

Prendiamo la funzione che hai scritto:

// calculate the norm
long double calculate_norm(std::vector<std::vector<double>> matrixA, int order)
{
    long double norm = 0;
    for (int i = 0; i != order; ++i)
    {
        for (int j = 0; j != order; ++j)
        {
            norm += matrixA[j] * matrixA[j];
        }
    }

    return sqrt(norm);
}

Come ben saprai la formula di Frobenius consiste nel calcolare la radice della somma dei quadrati di ogni scalare contenuto nella matrice, eppure vedendo il tuo codice non ti stai iterando sulla matrice bensì solo su una riga siccome stai utilizzando solo la varibaile j.

Sto utilizzando Visual Studio 2022 e non capisco come tu riesca a compilare questo codice siccome nella riga norm += matrixA[j] * matrixA[j]; quello che stai facendo è una moltiplicazione tra std::vector e non valori scalari, implementando quell'algoritmo correttamente in C++ avrai qualcosa come:

// calculate the norm
long double calculate_norm(std::vector<std::vector<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);
}

oppure in un modo più elegante:

// calculate the norm
long double calculate_norm(std::vector<std::vector<double>> matrix)
{
    long double norm = 0;

    for (auto& vec : matrix)
    {
        for (auto s : vec)
        {
            norm += s * s;
        }
    }

    return sqrt(norm);
}

Ho visto che hai commesso questo errore un pò ovunque nel tuo codice, ovvero quello di iterati su un vettore di vettori come se fosse un unico vettore.

Poi un altro errore che ho notato è che nella funzione random_val apri un file stream con std::ofstream os("RANDOM");, però non chiudi il file, per ragioni di performance le classi fstream scriveranno i dati sul disco solo dopo aver chiamato la funzione close() o altri metodi come flush(), finchè non lo fari i dati verrannò archiviati nella RAM ed il file sul disco rimarrà vuoto.

Nel caso di una tua prossima ripsosta ti prego di fornire il nome del compilatore che stai utilizzando, perchè mi pare strano che tu riesca a compilarlo.