[C++, Algoritmi] Metodo di Gauss con pivoting parziale

[lucia88]

Ciao a tutti,
ho scritto il seguente codice in c per risolvere sistemi lineari con MEG con pivoting parziale.
Però quando lo vado ad eseguire,dopo aver letto da tastiera coefficienti e termini noti, il programma si blocca...
Precisamente il pc mi dice 'il programma ha smesso di funzionare...'

//MEG con PIVOTING PARZIALE:
//legge una matrice e il vettore dei termini noti, risolve il sistema col metodo di gauss con
//pivoting parziale, infine calcola il residuo
//Durante l'eliminazione calcola det(A)

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<iostream>
using namespace std;

void controlla_puntatore(float*);
float* alloca_vector(int);
float** alloca_matrix(int);

void copia_sistema(float**,float*,float**,float*,int);

//legge la dimensione del sistema
int leggi_dim();

//legge i coefficienti  e i termini noti del sistema
void input_read(float**,float*,int);

//METODO DELLE SOSTITUZIONI ALL'INDIETRO:risolve il sistema con la matrice 
//in forma triangolare superiore e il vettore dei termini noti trasformato
//scrive le soluzioni in un altro vettore
float* sost_indietro(float**,float*,int);

//MEG+PIVOTING PARZIALE:
//trasforma il sistema in forma triangolare superiore e calcola det(A)
double meg(float**,float*,int);

//ERRORE: si sostituisce alle incognite il valore trovato
//e si calcola la differenza con i termini noti
double err(float**,float*,float*,float*,int);

void stampa_sistema(float**,float*,int);
void stampa_soluz(float*,int);
void stampa_errore(double,float*,int);


int main()
{
	cout<<"Risoluzione sistema mediante metodo di eliminazione di Gauss"<<endl;
	
	int dim; //incognite sistema=n equazioni=dimensione
	dim=leggi_dim();
	
	float** coeff;
	coeff=alloca_matrix(dim);
	float* b;
	b=alloca_vector(dim);
	input_read(coeff,b,dim);
	
	//copio i coefficienti e i termini noti in un altra matrice e vettore
	float** copia_coeff;
	copia_coeff=alloca_matrix(dim);
	float* copia_b;
	copia_b=alloca_vector(dim);
	copia_sistema(coeff,b,copia_coeff,copia_b,dim);
	
	
	cout<<"Risoluzione sistema in corso..."<<endl;
	//trasformo il sist in triangolare superiore e calcolo il determinante
	double det;
	det=meg(coeff,b,dim);
	cout<<"Determinante matrice:"<<det<<endl;
	
	//calcolo soluzione con metodo delle sostituzioni all'indietro
	float* x;
	x=alloca_vector(dim);
	x=sost_indietro(coeff,b,dim);
	
	//stampo il sistema originale e la soluzione calcolata
	stampa_sistema(copia_coeff,copia_b,dim);
	stampa_soluz(x,dim);
	
	//ERRORE
	double err_max;
	float* error;
	error=alloca_vector(dim);
	err_max=err(copia_coeff,copia_b,x,error,dim);
	stampa_errore(err_max,error,dim);
	
	return 0;	
}

//DEFINIZIONE FUNZIONI

void controlla_puntatore(float* p)
{
	if(p!=NULL) return;
	else
	{
		cout<<"Memoria non allocata correttamente...";
		exit(1);
	}
}

float* alloca_vector(int n)
{
	float* v;
	v=new float[n];
	controlla_puntatore(v);
	return v;
}

float** alloca_matrix(int n)
{
	float** a;
	a=new float*[n];
	for(int i=0;i<n;i++)
	{
		a[i]=new float[n];
		controlla_puntatore(a[i]);
	}
	return a;
}

void copia_sistema(float** a,float* b,float** cop_a,float* cop_b,int n)
{
	for(int i=0;i<n;i++)
	{
		cop_b[i]=b[i];
		for(int j=0;j<n;j++) cop_a[i][j]=a[i][j];
	}
	return;
}

void stampa_sistema(float** A,float* b,int dim)
{
	cout<<endl;
	cout<<"Sistema originale:"<<endl;
	
	for(int i=0;i<dim;i++)
	{
		for(int j=0;j<dim;j++)
		{
			cout<<"("<<A[i][j]<<")x"<<j+1;
			if(j==dim-1) cout<<"=";
			else cout<<"+";
		}
		cout<<b[i]<<endl;
	}
	return;
}


void stampa_soluz(float* s,int n)
{
	cout<<endl<<"La soluzione del sistema e' "<<endl;
	for(int i=0;i<n;i++)
	{
		cout<<"x"<<i+1<<"="<<s[i]<<endl;
	}
	return;
}

void stampa_errore(double max,float* e,int n)
{
	cout<<endl<<"L'errore massimo commesso nella risoluzione del sistema e':"<<endl;
	cout<<max;
	
	cout<<endl<<"Vettore residuo=(";
	for(int i=0;i<n;i++)
	{
		cout<<e[i];
		if(i==n-1) cout<<")";
		cout<<",";
	}
	return;
}

int leggi_dim()
{
	int n;
	cout<<endl<<"Inserire il numero di equazioni e di incognite del sistema:";
	cin>>n;
	while(n<=0)
	{
		cout<<endl<<"Errore! Il numero inserito e' negativo...Riprovare."<<endl;
		cout<<"n=";
		cin>>n;
	}
	return n;
}

void input_read(float** c,float* b,int n)
{
	cout<<endl<<"Inserire il sistema da risolvere:"<<endl;
	for(int i=0;i<n;i++)
	{
		cout<<"EQUAZIONE"<<i+1<<":"<<endl;
		for(int j=0;j<n;j++)
		{
			cout<<"Coefficiente di x"<<j+1<<"=";
			cin>>c[i][j];
		}
		cout<<"Termine noto:";
		cin>>b[i];
	}
return;
}

//MEG con PIVOTING PARZIALE(trasformazione in triangolare superiore)
double meg(float** A,float* b,int n)
{
	double m,M,aux,det=1;
	int index,scont=0;
	
	//ciclo che fa le (n-1) iterazioni necessarie per eliminare
	//x_1,x_2,...,x_(n-1) dalle (di volta in volta) restanti equaz.
	//alla fine elimino (n-1) equ.
	for(int k=0;k<n-1;k++)
	{
		//ricerca elemento pivotale massimo:
		//scambio la riga che lo contiene con la riga k-esima
		M=fabs(A[k][k]);
		index=k;
		for(int i=k+1;i<n;k++){
			if(M<fabs(A[i][k])) 
			{
				M=fabs(A[i][k]);
				index=i;
			}
		}
		//se ho trovato un elemento pivotale maggiore scambio le righe tra loro
		if(index>k)
		{
			scont++; //contatore:indica scambi di riga effettuati
			
			//scambio di 2 vettori riga
			for(int j=k;j<n;j++)
			{
				aux=A[k][j];
				A[k][j]=A[index][j];
				A[index][j]=aux;
			}
			//scambio dei termini noti delle 2 equazioni
			aux=b[k];
			b[k]=b[index];
			b[index]=aux;
		}
		//quindi al k posto c'è l'eq.con il massimo k-esimo coeff
		//procedo con il MEG
		
		
		for(int i=k+1;i<n;i++)
		{
			m=A[i][k]/A[k][k];
			
			for(int j=k+1;j<n;j++)
			A[i][j]=A[i][j]-m*A[k][j];
			
			b[i]=b[i]-m*b[k];
		} 
	}
	
	
	//calcolo determinante della matrice
	for(int k=0;k<n;k++) det=det*A[k][k];
	//per ogni scambio di righe effettuato moltiplico il determinante per -1
	det=det*pow(-1,scont);
	cout<<endl<<"scambi di riga effettuati:"<<scont;
	return det;
}


float* sost_indietro(float** A,float* b,int n)
{
	float* x;
	x=alloca_vector(n);
	
	//somma parziale
	double s;
	
	//calcolo l'ultimo valore cioe x_n=b_n/a_n,n
	x[n-1]=b[n-1]/A[n-1][n-1];
	
	//calcolo delle incognite precedenti cioè x_1,x_2....
	for(int i=n-2;i>=0;i--)
	{
		s=0;
		for(int j=1+i;j<n;j++)
		s=s+A[i][j]*x[j]; // somma incognite trovate in precedenza moltipl. per loro coeff.
		
		// x_n=(b_n-s)/a_n,n
		x[i]=(b[i]-s)/A[i][i];
	}
	
	return x;
}


double err(float** A,float* b,float* x,float* e,int n)
{
	double max;
	
	for(int i=0;i<n;i++)
	{
		e[i]=0;
		for(int j=0;j<n;j++)
		{
			e[i]=e[i]+x[j]*A[i][j];
		}
		e[i]=fabs(b[i]-e[i]); 
		//ho calcolato l'errore di ciascuna riga come differenza tra il termine noto 
		//e la somma delle incognite per i rispettivi coefficienti 
	}
	
	//ricerca errore massimo
	max=fabs(e[0]);
	for(int i=1;i<n;i++)
	if(max<fabs(e[i])) max=fabs(e[i]);
	
	return max;
}

[apatriarca]

Alcune osservazioni generali mentre cerco il problema nel codice:
1. Usare puntatori a puntatori per rappresentare matrici multidimensionali non è molto efficiente e complica incredibilmente la gestione della memoria. Normalmente si preferisce allocare un singolo array monodimensionale e accedere ad esso con qualche calcolo. Nota che è così che gli array bidimensionali allocati staticamente appaiono in memoria.
2. È possibile dichiarare e inizializzare una variabile sulla stessa riga. È anzi normalmente preferibile (più leggibile).
3. Non verifichi che l'allocazione sia andata a buon fine nel caso dell'allocazione del vettore di puntatori. Ma in effetti in C++ viene generata una eccezione (se le eccezioni sono attive ovviamente) quando l'allocazione fallisce.
4. Il termine pow(-1,scont) è equivalente a ( (scont & 1) ? -1.0 : 1.0 ). Ma un modo migliore per implementare questo termine è quello di inizializzare un valore a 1 e poi moltiplicare per -1 ad ogni scambio.

Non ho provato a compilarlo, però mi sembra corretto. Hai provato ad usare un debugger per capire esattamente a che punto del programma va in crash?

[vict85]

Al di là di ciò che ti ha detto antonio aggiungo 2 cose:

1. L'errore numero (1) è qui

for(int i=k+1;i<n;k++){
         if(M<fabs(A[i][k])) 
         {
            M=fabs(A[i][k]);
            index=i;
         }
      }

Anche se è difficile da vedere l'errore consiste nell'aumentare k invece che i.

2. Una funzione dovrebbe fare una sola cosa. Quindi non dovresti calcolare il determinante all'interno della funzione che ti serve per risolvere sistemi. Al limite dovresti far ritornare \(\pm 1\).

3. È comodo che tu non debba inserire i coefficienti attraverso la console. Dovresti quindi creare una funzione che te ne fa una più o meno casuale oppure che ti legga la matrice da un file (è molto più comodo scrivere su un file di testo ). Io per i test ne ho creata una che me ne faceva una particolare (tra l'altro non negativa, simmetrica e a diagonale dominante ma ovviamente non è possibile usarne una così per controllare il pivoting).

[vict85]

Dopo quella correzione funzione tutto. Ho però un paio di commenti ancora.

1) Nel pivoting dovresti fare

if ( abs(A[i][k]) > M )

invece di

if(M < fabs(A[i][k]))

Oppure usare il <= perché dovresti evitare di cambiare riga se non ce n'è bisogno.

2) Stai usando una metodo di rappresentazione della memoria inefficiente. D'altra parte nell'unica cosa in cui questo metodo è tremendamente più efficiente tu eviti di utilizzarlo: nello scambio delle righe. Questo codice fa lo scambio delle righe:

if (index != k)
		{
			//scambio di 2 vettori riga
			double * t = A[k];
			A[k] = A[index];
			A[index] = t;
			
			//scambio dei termini noti delle 2 equazioni
			aux = b[k];
			b[k] = b[index];
			b[index] = aux;
		}