/*Programma che legge una dimensione n minore
di 15, gli elementi di una matrice A, del vettore b
e della tolleranza 0<t<=10^-4 e approssima la soluzione
del sistema lineare con il metodo di Jacobi con arresto a
posteriori e un massimo di 100 iterazioni*/
#include <stdio.h>
#include <math.h>
#define N 15
#define MAX_ITER 100
#define tol 0.0001
int leggi_dato (int&);
int controllo_dato (double);
int controllo_tolleranza (double);
void leggi_matrice (double [N][N], int);
void leggi_vettore (double[], int);
double leggi_tolleranza (double&);
double modulo (double);
double costante_lambda (double [N][N], int);
int algoritmo_Jacobi (int, double [N][N], double [], double [], double, double);
double norma_max (int, double [], double []);
int main(){
int n= leggi_dato (n);
if (controllo_dato (n)==0) {printf ("Numero errato\n"); return 0;}
double A[N][N];
leggi_matrice (A, n);
double b[n];
leggi_vettore (b, n);
printf ("Matrice A=\n");
for (int i=0; i<n; i++){for (int j=0; j<n; j++) printf ("%.2lf\t", A[i][j]); printf ("\n");}
double t=leggi_tolleranza (t);
if (controllo_tolleranza (t)==0) {printf ("Numero errato\n"); return 0;}
double l=costante_lambda (A, n);
if (l>=1-tol) {printf ("Diagonale non strettamente dominante"); return 0;}
double eps=l/(1-l);
double x[N], xk[N];
for (int i=0; i<n; i++) x[i]=1;
int iterazioni= algoritmo_Jacobi (n, A, b, x, eps, t);
printf ("\nNumero di iterazioni compiute=%d", iterazioni);
return 0;}
int leggi_dato (int&a)
{printf ("Inserire numero minore di 16="); scanf("%d", &a);
return a;}
int controllo_dato (double a)
{if (a>15||a<1) return 0; else return 1;}
int controllo_tolleranza (double a)
{if (a>0.0001||a<=0) return 0; else return 1;}
void leggi_matrice (double A[N][N], int n)
{for (int i=0; i<n; i++) {for (int j=0; j<n; j++) {printf ("Matrice %d %d=", i+1, j+1); scanf("%lf", &A[i][j]);}}
return;
}
void leggi_vettore (double a[], int n)
{for (int i=0; i<n; i++) {printf ("Inserire %d componente del vettore b=", i+1);
scanf ("%lf", &a[i]);} return;
}
double leggi_tolleranza (double&a)
{printf ("Inserire tolleranza 0<t<=0.0001 t="); scanf("%lf", &a);
return a;}
double modulo (double a)
{if (a<0) return -a;
else return a; }
double costante_lambda (double a[N][N], int n)
{double l;
double lambda[N];
l=lambda[0];
for (int i=0; i<n; i++) {for (int j=0; j<n; j++)
{if (j!=i) lambda[i]+=modulo (a[i][j]);
lambda[i]=lambda[i]/modulo(a[i][i]);}
if (lambda[i]>l) l=lambda[i];}
return l;}
int algoritmo_Jacobi (int n, double a[N][N], double b[], double x[], double eps, double t)
{int cont=0;
double xk[N];
while ((cont<MAX_ITER)&&(norma_max(n,x, xk)*eps>=t))
{ double s=0; for (int i=0; i<n; i++){for (int j=0; j<n; j++)
{if (j!=i) {s+=a[i][j]*xk[j];} xk[i]=(b[i]-s)/a[i][i];}
x[i]=xk[i];}
cont++;}
for (int i=0; i<n; i++) printf ("x=%.2lf\n", xk[i]);
return cont;
}
double norma_max (int n, double w[], double v[])
{double abs[N], diff [N];
for (int i=0; i<n; i++) {diff [i]= w[i]-v[i];}
double normainf= abs[0];
for (int i=0; i<n; i++) {abs[i]= modulo(diff[i]);
if (abs[i]>normainf) normainf=abs[i];
}
return normainf;}