/*
   Calcola la matrice di iterazione del metodo
   di Gauss Seidel e di Jacobi, e ne calcola
   il raggio spettrale col metodo delle potenze.

Parametro sulla linea di comando:
nomefile     -> legge la matrice dal file ('-' per tastiera)
n            ->  matrice di Hilbert n per n
n a b c      ->  matrice tridiagonale n per n con elementi [a b c]
*/

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

#define PREC (5e-16f)

/*
   Matrice J
   La diagonale e' nulla, per il resto
   basta dividere ogni riga di A per l'elemento
   diagonale e cambiare segno.
*/
void MatriceJacobi(int n,const matrice a, matrice m)
{
  int i, j;
  for (i=0; i<n; i++)
    for (j=0; j<n; j++)
      if (i==j) m[i][j]=0.;
      else m[i][j] = -a[i][j]/a[i][i];
}

/*
   Matrice G
   Calcolata con una "forward substitution" dove la matrice
   da invertire e' la parte triangolare inferiore di A
   e il "termine noto" e' la parte strettamente superiore.
   Nota: prima riga di G == prima riga di J, inoltre
   la prima colonna di G e' nulla.
*/
void MatriceGaussSeidel(int n,const matrice a, matrice g)
{
  int i, j, k;
  double sum;

  for (i=0; i<n; i++)
    for (k=0; k<n; k++) {
      if (i<k) sum=a[i][k]; else sum=0.;
      for (j=0; j<i; j++) sum += a[i][j]*g[j][k];
      g[i][k] = -sum/a[i][i];
    }
}

#define MAXITER 200
#define LOOSEPREC (PREC*8192)

/*
   Metodo delle potenze.
   Spesso non converge (ad es. perche' ci sono piu' autovalori
   con lo stesso modulo) e restituisce -1.
   ATTENZIONE: -1 non puņ essere confuso con un vero raggio
   spettrale, essendo un numero negativo!
*/
double RaggioSpettrale(int n,const matrice a)
{
  int i, j, niter=0;
  double sum, norma, *x, *newx, beta=0, oldbeta=1;
  x = calloc(n,sizeof *x);
  newx = calloc(n,sizeof *newx);
  for (i=0; i<n; i++) x[i]=1;

  while ( ++niter<MAXITER && fabs(beta-oldbeta)>LOOSEPREC ) {
    oldbeta=beta;
    norma = -1;
    for (i=0; i<n; i++) {
      sum = 0.;
      for (j=0; j<n; j++) sum += a[i][j]*x[j];
      newx[i] = sum;
      if ( fabs(sum)>norma ) norma=fabs(sum);
    }
    beta = newx[0]/x[0];
    for (i=0; i<n; i++) x[i]=newx[i]/norma;
  }
  if (niter==MAXITER) return -1;
  return fabs(beta);
}

int main(int argc, char *argv[])
{
  int n=0;
  matrice a, l;

	/* 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);

  MatriceJacobi(n,a,l);
  printf("Matrice J (raggio spettrale %f):\n",RaggioSpettrale(n,l));
  Stampa(n,l);

  MatriceGaussSeidel(n,a,l);
  printf("Matrice G (raggio spettrale %f):\n",RaggioSpettrale(n,l));
  Stampa(n,l);

        /* Libera memoria utilizzata */

  libera_matrice(n,a);
  libera_matrice(n,l);

  return 0;
}