// ............................................................. seidelrelax.cc
//                                                         versione 08.APR.2003
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "corso.h"

// ================================================================ seidelrelax

void SeidelRelax(double **a,double *b,double *x,int n,double accuracy,
		 double relax)
{
  int i,ii,j;
  double delta;
  double *dx,*xx;

  // .......................................................... memoria per xx
  dx=new double [n];
  xx=new double [n];
  // ................................................ inizializza la soluzione
  for (i=0;i<n;i++) x[i]=0;
  // ..................................................... soluzione iterativa
  for (ii=0;ii<10000;ii++)
    { // .......................................................... deviazione
      delta=0;
      // ........................................................ nuovi valori
      for (i=0;i<n;i++)
	{ // .................................................... termine noto
	  xx[i]=b[i];
	  // ........................................... termini non diagonali
	  for (j=0;j<n;j++) if (j!=i) xx[i]-=a[i][j]*x[j];
	  // ............................................... termine diagonale
	  xx[i]/=a[i][i];
	  dx[i]=xx[i]-x[i];
	  // ...................................................... deviazione
	  delta+=pow((xx[i]-x[i]),2.0);
	  x[i]+=relax*dx[i];
	}
      // ......................................... mostra i risultati parziali
      printf("\n");
      for (i=0;i<n;i++) printf("x[%d]=%.8le\n",i,x[i]);
      // ............................................................... stop?
      if (sqrt(delta/(double) n)<accuracy) break;
    }
  printf("\n%d iterazioni\n\n",ii+1);
  // .............................................................. pulisci xx
  delete [] xx;
  delete [] dx;
}