#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define MAXDIM 4
typedef double matrice [MAXDIM][MAXDIM];
typedef double vettore[MAXDIM];
void stampavettore(vettore,int);
void riduzionegauss(matrice,vettore,int);
int ricercapivot(matrice,int,int);
double modulo(double);
void scambio(double&,double&);
void scambiorighe(matrice,vettore,int,int,int);
void soluzione_sistematriangolaresuperiore(matrice,vettore,vettore,int);
int main()
{int n;
do{printf("Inserire la dimensione del sistema lineare Ax=b:\n");scanf("%d",&n);} while(n!=3);
matrice A;
vettore b;
printf("Matrice A={{3,3,3},{1,3,1},{2,-1,3}}\n");
//costruzione matrice A:
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
printf("Inserire l'elemento a[%d][%d] della matrice:\n",i+1,j+1); scanf("%lf",&A[i][j]);
}
}
//costruzione vettore b(le cui componenti sono le somme degli elementi di ogni riga di A):
for(int i=0;i<n;i++) {
for(int j=0;j<n;j++)
{b[i]+=A[i][j];}
}
printf("Vettore b dei termini noti:\n");
for(int i=0;i<n;i++)
printf("%lf\t",b[i]);
vettore x;
riduzionegauss(A,b,n);
printf("\nMatrice dopo riduzione Gauss:\n");
for(int i=0;i<n;i++){printf("\n");
for(int j=0;j<n;j++){
printf("%lf\t",A[i][j]);
}
}
printf("\nVettore dopo riduzione Gauss:\n");
for(int i=0;i<n;i++)
printf("%lf\t",b[i]);
return 0;
}
//Metodo della riduzione di Gauss(con riordinamento pivotale)applicato al sistema lineare Ax=b:
void riduzionegauss(matrice a,vettore v,int d)
{
double m;//coefficienti moltiplicatori di ogni riga della matrice
int k_pivot; //indice della riga in cui si trova il pivot
for(int k=0;k<d-1;k++) /*ordine della riduzione delle righe:in particolare,partendo dalla riga successiva alla prima(che è la riga 0)si eseguono operazioni
sulle rimanenti righe,per d-1 volte(e arrivando quindi all'ordine (d-2) compreso)*/
{k_pivot=ricercapivot(a,d,k);
if (k_pivot!=k)
scambiorighe(a,v,d,k,k_pivot); //scambio tra loro delle righe k-esima e k_pivot-esima(cioè quella in cui c'è il pivot)
for(int i=k+1;i<d;i++) //si parte considerando la "sottomatrice" della matrice di partenza(per far annullare gli elementi sotto la k-esima colonna)
{m=a[i][k]/a[k][k]; //calcolo del moltiplicatore di ogni k-esima riga(della matrice di partenza e di ogni sottomatrice che si creerà dopo)
for(int j=k+1;j<d;j++)
{
a[i][j]=a[i][j]-m*a[k][j];/*calcolo dell'elemento a[i][j] della matrice ridotta, al passo k-esimo:gli elementi sotto la prima colonna di ogni
sottomatrice dovranno annullarsi*/
}
v[i]=v[i]-m*v[k]; //calcolo dell'elemento v[i] della colonna dei termini noti ridotta, al passo k-esimo
}
}
return;
}
//Calcolo dei moduli degli elementi della matrice(per poi trovare il pivot,cioè l'elemento con maggior valore assoluto sulla colonna):
double modulo(double x)
{if(x<0) return -x;
else return x;
}
//Ricerca dell'indice della riga in cui è presente il pivot:
int ricercapivot(matrice a,int d,int k)
{ double mk,max;
int indice_pivot=k;
for(int i=k+1;i<d;i++){
max=modulo(a[k][k]);
mk=modulo(a[i][k]);
if(mk>max)
{max=mk;
indice_pivot=i;
}
}
return indice_pivot;
}
//Funzione scambio(per poi scambiare tra loro le componenti delle due righe k-esima e k_pivot-esima della matrice):
void scambio(double& x,double& y){
double aux=x;x=y;y=aux;
return;
}
//Funzione che scambia tra loro le righe k-esima e k_pivot-esima della matrice:
void scambiorighe(matrice a,vettore v,int d,int k,int k_pivot)
{
for(int j=k;j<d;j++) //scorrendo sulla colonna k-esima,dalla componente k-esima fino all'ultima,si scambiano gli elementi della riga k e della riga k_pivot
scambio(a[k][j],a[k_pivot][j]);
scambio(v[k],v[k_pivot]);
return;
}
//Funzione che risolve il sistema triangolare superiore (ottenuto dopo la riduzione di Ax=b con l'algoritmo di Gauss) tramite l'algoritmo backward:
void soluzione_sistematriangolaresuperiore(matrice a,vettore v,vettore x,int d)
{for(int i=d-1;i>=0;i--) //partendo dall'ultima riga del sistema,si trova la componente x[i] del vettore soluzione x[d] del sistema:si fa poi l'algoritmo backward
{x[i]=v[i];
for(int j=i+1;j<d;j++)
x[i]-=a[i][j]*x[j];
x[i]/=a[i][i];
}
return;
}