/*
  Iterativi.c
  Jacobi + GaussSeidel

Risolve sistemi lineari con i metodi iterativi classici.
Sono possibili ancora ottimizzazioni nei calcoli.

Riga di comando:
- "solve 10"  crea una matrice di Hilbert 10x10
- "solve 10 1 2 3" crea una tridiagonale 10x10 con elementi 1, 2, 3
  (quello centrale è l'elemento diagonale)
- "solve nomefile" legge la matrice dal file `nomefile' (prima la
  dimensione, poi tutti gli elementi per riga)
  Usare nomefile='-' per leggere da stdin

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "matutils.h"

#define EPS (5e-13f)
#define MAXITER 50

/*
   Metodo iterativo di Jacobi per risolvere Ax=b
   Restituisce 0 se non c'e' convergenza, altrimenti
   il numero di iterazioni effettuate.
   Controllo di convergenza in norma infinito:
          |x - xold|<EPS
*/

int JacobiSolve(int n, const matrice a, const double *b,
		double *x)
{
  int i,j, niter=0;
  double sum, *oldx, diff=1.0;
  oldx = calloc(n,sizeof *oldx);
  for (i=0; i<n; i++) { oldx[i]=x[i]=1.0; }

  while (diff>EPS && ++niter<MAXITER) {
    diff = 0.0;
    for (i=0; i<n; i++) {
      sum = b[i];
      for (j=0; j<n; j++)
	if (j!=i) sum -= a[i][j]*oldx[j];
      x[i] = sum / a[i][i];
      diff += fabs(x[i]-oldx[i]);
    }
    for (i=0; i<n; i++) oldx[i] = x[i];
    /*  printf("%d %g  ",niter,diff); StampaVet(n,x);  */
  }

  free(oldx);
  if (niter==MAXITER) return 0;
  return niter;
}

/*
   Metodo iterativo di Gauss Seidel per risolvere Ax=b
   Restituisce 0 se non c'e' convergenza, altrimenti
   il numero di iterazioni effettuate.
   Nota: gia' per n=4 e A=matrice di Hilbert e' rho(G)=0.999
   e praticamente non c'e' convergenza (anche se teoricamente si')
   * Inserire peso "w" per ottenere il metodo SOR
*/

int GaussSeidelSolve(int n, const matrice a, const double *b,
		      double *x)
{
  int i,j, niter=0;
  double sum, oldxi, diff=1.;
  for (i=0; i<n; i++) x[i]=1.;

  while (diff>EPS && ++niter<MAXITER) {
    diff = 0.0;

    for (i=0; i<n; i++) {
      sum = b[i];
      for (j=0; j<n; j++)
	if (j!=i) sum -= a[i][j]*x[j];
      oldxi = x[i];
      x[i] = sum / a[i][i];               /* <- Inserire peso "w" */
      diff += fabs(x[i]-oldxi);
    }
    /*  printf("%d %g  ",niter,diff); StampaVet(n,x);  */
  }

  if (niter==MAXITER) return 0;
  return niter;
}

int main(int argc, char *argv[])
{
  int i, n;
  matrice a, l;
  double *x1, *x2;

	/* Inizializzazione della matrice */

  if (argc!=2 && argc!=5) {
    fprintf(stderr,"Uso:\n"
	    "\t%s <nomefile>       Legge matrice dal file (usare '-' per stdin)\n"
	    "\t%s <N>              Matrice di Hilbert NxN\n"
	    "\t%s <N> <A> <B> <C>  Matrice tridiagonare NxN con elementi [A B C]\n",
	    argv[0],argv[0],argv[0]);
    return 1;
  }
  n = atoi(argv[1]);
  if (argc==2) {
    if (n==0)
      a = LeggiFile(argv[1],&n);
    else {
      a = alloca_matrice(n);
      hilbert(n,a);
    }
  } else {
    a = alloca_matrice(n);
    tridiagonale(n, atof(argv[2]), atof(argv[3]), atof(argv[4]), a);
  }
  

  l=alloca_matrice(n);

  printf("Matrice iniziale:\n"); 
  Stampa(n,a);

  x1 = calloc(n, sizeof *x1);
  for (i=0; i<n; i++) x1[i]=1.;
  x2 = calloc(n, sizeof *x2);

        /* Sistema lineare risolto con Jacobi */

  printf("Sol. di A.x = "); StampaVet(n,x1);
  i = JacobiSolve(n,a,x1,x2);
  if ( !i ) printf("Jacobi non converge\n");
  else { printf("Jacobi (%d it): ",i); StampaVet(n,x2); }

        /* Sistema lineare risolto con GaussSeidel */

  printf("Sol. di A.x = "); StampaVet(n,x1);
  i = GaussSeidelSolve(n,a,x1,x2);
  if ( !i ) printf("GaussSeidel non converge\n");
  else { printf("GaussSeidel (%d it): ",i); StampaVet(n,x2); }

        /* Libera memoria utilizzata */

  free(x2); free(x1);
  libera_matrice(n,a);
  libera_matrice(n,l);
  return 0;
}