/**
   * Az Ax=b egyenletrendszer megoldasahoz a prekondicionalt_konjugalt_gradiens módszer.
   * x* az Ax*=b megoldása.
   */
   
   /*
   * Bemenet:
   * n
   * A
   * H
   * b
   * x0
   * maxit
   * e
   */
   
   /* Ahol
   * n       : nxn -es mátrix.
   * A       : mátrix
   * H       : prekondícionált mátrix
   * b       : Ax=b egyenletrendszer jobb oldala.
   * x0      : kezdővektor
   * maxit   : maximális iteráció száma.
   * e       : küszöbszám. xk és xk+1 közötti különbség ha ettől kisebb, akkor vége.
   */
   
typedef struct sor
{
       int oszlop;
       double ertek;
       struct sor *kov;
} SOR;

void Prekondicionalt_Konjugalt_gradiens(char in[], char out[]);

void beolvas2 ( char in[], SOR ** A, SOR **H, double* b, double *x0, double *xk, int* maxit, double * e );
void rendez( SOR ** f, int n);
void BESZUR_ELEJE ( SOR ** f, int o, double e );
void Beszuro_rendezes(double *A, int *oszlop, int hossz);
void beolvas ( char in[], SOR ** A, SOR ** H, /*SOR **Temp,*/ double* b, double *x0, double *xk, int* maxit, double * e );
void BESZUR_VEGE ( SOR ** f, int o, double e );
double ELEM ( SOR ** f, int j );
void G_0( SOR **sorok, double *x0, double *b, int n, double *eredmeny1, double *eredmeny2);
void HxHT ( SOR ** f, SOR ** eredmeny, int n );
double LAMBDA_k( SOR **sorok, double *ek, double *gk, double *hk, int n);
void G_k( SOR **sorok, double *gk, double lambda, double *hk, int n, double *eredmeny);
void TOROL_LISTA ( SOR ** f );

void Eredmeny_kiiras( char out[], int n, double * vektor, int it );
double Ellenorzes_legnagyobb_kulonbseg( double * x, int n);
double Euklideszi_norma_vektorra( double * vektor, int o );
void lin_e_r( SOR **H, double *b, int n );
void LDLT( SOR **sorok/*double *A*/, double *b, int n, double *eredmeny, int hiba );
void R_i_first( SOR **sorok, int j, int n, double *eredmeny);
void R_i( SOR **sorok, int j, int k, int n, double *eredmeny);
void X_j( SOR **sorok,  int i, int n, double *eredmeny);
void BESZUR ( SOR ** f, int s, int o, double e );
int szimmetrikus ( double *A, int n);
void Ellenorzes(double *A, int n, double *b, double * x);

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <sys\timeb.h>
#include <string.h>
#include <ctype.h>
/*#include "leak_detector_c.h"*/

int hiba=0;

struct timeb t_start, t_current;
     int t_diff_in_sec, t_diff_in_ms;

int main (int argc, char* argv[]){
    
     if (argc<3){
      puts("Hiba! Kevés argumentumot adott meg.");
      puts("A program helyes futtatása:\nprogramnév bemeneti_fájl_neve_kiterjesztéssel kimeneti_fájl_neve_kiterjesztéssel\n");
      return 1;
     }
     else if (argc>3){
       puts("Hiba! Túl sok argumentumot adott meg.");
       puts("A program helyes futtatása:\nprogramnév bemeneti_fájl_neve_kiterjesztéssel kimeneti_fájl_neve_kiterjesztéssel\n");
       return 1;
     }

     /*struct timeb t_start, t_current;
     int t_diff_in_sec, t_diff_in_ms;*/
     ftime(&t_start);
    
    Prekondicionalt_Konjugalt_gradiens(argv[1], argv[2]);
    
      ftime(&t_current);
      t_diff_in_sec = (int)  (t_current.time - t_start.time);
      t_diff_in_ms = (int ) (t_current.millitm - t_start.millitm);
    
    FILE *fout;
    fout = fopen(argv[2], "a");    
    printf("\nEltelt idő: %d ms",t_diff_in_sec*1000+t_diff_in_ms);
    fprintf(fout, "\nEltelt idő: %d ms",t_diff_in_sec*1000+t_diff_in_ms);
    fclose(fout);
    
    /*atexit(report_mem_leak);*/
    return 0;
}
/*____________________________________________________________________________*/
void Prekondicionalt_Konjugalt_gradiens(char in[], char out[]){
    int n, i, j, k, m, maxit;
    n = i = j = k = m = maxit = 0;
    FILE *fin, *fout;
    
    fin = fopen(in, "r");
    fout = fopen(out, "w");
    
    fscanf(fin, "%d\n",&n);
    
    double *b=(double *) calloc( n, ( sizeof ( double ) ) );
    double *x0=(double *) calloc( n, ( sizeof ( double ) ) );
    double *xk=(double *) calloc( n, ( sizeof ( double ) ) );
    double *h0=(double *) calloc( n, ( sizeof ( double ) ) );
    double *hk=(double *) calloc( n, ( sizeof ( double ) ) );
    double *g0=(double *) calloc( n, ( sizeof ( double ) ) );
    double *gk=(double *) calloc( n, ( sizeof ( double ) ) );
    
    double *e0, *ek;
    e0=(double *) calloc( n, ( sizeof ( double ) ) );
    ek=(double *) calloc( n, ( sizeof ( double ) ) );
    
    double mk, lk, e;
    SOR *A_sorai[n], *H_sorai[n]/*, *HT_sorai[n]*/; /* Temp ben lesz H*HT letárolva */
    
    /*----------------------------------------*/
    fclose(fin);
    beolvas2 ( in, A_sorai, H_sorai, /*Temp_sorai,*/ b, x0, xk, &maxit, &e );
    ftime(&t_current);
      t_diff_in_sec = (int)  (t_current.time - t_start.time);
      t_diff_in_ms = (int ) (t_current.millitm - t_start.millitm);
    printf("\nbeolvasási idő: %d ms",t_diff_in_sec*1000+t_diff_in_ms);
    
	/*----------------------------------------*/
	/* Az algoritmus */
	
	/* g0=A*x0-b h0=-e0*/
	G_0( A_sorai, x0, b, n, g0, gk);
  
  /* e0 megh LDLT seg.vel */
  hiba = 0;
  for ( i=0; i<n; i++ )
    e0[i] = g0[i];
  lin_e_r( H_sorai, e0, n );

  if ( hiba == 1 ){
      puts("hiba"); fputs("hiba", fout);
      fflush(fout);
      fclose(fout);
      /*free(A);*/ /*free(HHT); free(Temp);*/ free(b); free(x0); free(xk); free(h0); free(hk); free(g0); free(gk); free(e0); free(ek);
      for( i=0; i<n; i++ ) TOROL_LISTA(&A_sorai[i]);
      for( i=0; i<n; i++ ) TOROL_LISTA(&H_sorai[i]);
      /*for( i=0; i<n; i++ ) TOROL_LISTA(&Temp_sorai[i]);*/
      return;
    }
  
	for ( i=0; i<n; i++ ){
    h0[i] = e0[i]*(-1);
  }
  
	for ( m=0; m<=maxit; m++ ){
        /* [?||gk||<=e?[stop, eredmény: xk] */
        double norma = Euklideszi_norma_vektorra ( g0, n );
        norma = sqrt (norma);
        if ( fabs( norma ) >= e ){
          double nevezo=0.0, szamlalo=0.0;
          /*double temp[n];*/
          double *temp=(double *) calloc( n, ( sizeof ( double ) ) );
          
          /* lk=(ek,gk)/(Ahk,hk) */
          lk = LAMBDA_k( A_sorai, e0, g0, h0, n );
                    
          /*xk+1=xk+lkhk*/ /*xk+1 az xk    xk az x0*/
          for ( i=0; i<n; i++ )
            temp[i] = lk*h0[i];
          for ( i=0; i<n; i++ )
            xk[i] = x0[i] + temp[i];
          
          /*gk+1=gk+lkAhk*/
          G_k(A_sorai, g0, lk, h0, n, gk);
                     
          /* HHT ek+1 = gk+1 */
          /* TODO ek+1 meghatározása.*/
          /*HxHT(H_sorai, Temp_sorai, n);*/
          
          /* ek+1 megh LDLT seg.vel */
          /*ek = LDLT( (double *)Temp, gk, n );*/
          hiba = 0;
          /*puts("ldlt eleje");*/
          for ( i=0; i<n; i++ )
            ek[i] = gk[i];
          lin_e_r( H_sorai, ek, n );
                   
          if ( hiba == 1 ){
            puts("hiba"); fputs("hiba", fout);
            fflush(fout);
            fclose(fout);
            /*free(A);*/ /*free(HHT); free(Temp);*/ free(b); free(x0); free(xk); free(h0); free(hk); free(g0); free(gk); free(e0); free(ek);
            /*free(temp);*/
            for( i=0; i<n; i++ ) TOROL_LISTA(&A_sorai[i]);
            for( i=0; i<n; i++ ) TOROL_LISTA(&H_sorai[i]);
            /*for( i=0; i<n; i++ ) TOROL_LISTA(&Temp_sorai[i]);*/
            return;
          }
           
          /* mk = (ek+1,gk+1)/(ek,gk) */    /* gk+1 az gk   gk az g0*/
          szamlalo = nevezo = 0.0;
          for ( i=0; i<n; i++ )
            szamlalo += ek[i] * gk[i];
          for ( i=0; i<n; i++ )
            nevezo += e0[i] * g0[i];
          mk = szamlalo / nevezo;
          /*hk+1 = -ek+1 + mkhk*/   /* dk+1 az dk   dk az d0*/
          for ( i=0; i<n; i++ )
            temp[i] = mk * h0[i];
          for ( i=0; i<n; i++ )
            hk[i] = (ek[i] * (-1)) + temp[i];
            
          for ( i=0; i<n; i++ ){
            x0[i] = xk[i];
            g0[i] = gk[i];
            h0[i] = hk[i];
            e0[i] = ek[i];
          }
          /*free(ek);*/
          free(temp);
        }
        else {
          fflush(fout);
          Eredmeny_kiiras( out, n, x0, m );
          /*Ellenorzes( (double *)A, n, b, x0 );*/
          fflush(fout);
          fclose(fout);
          /*free(A);*/ /*free(HHT); free(Temp);*/ free(b); free(x0); free(xk); free(h0); free(hk); free(g0); free(gk); free(e0); free(ek);
          for( i=0; i<n; i++ ) TOROL_LISTA(&A_sorai[i]);
          for( i=0; i<n; i++ ) TOROL_LISTA(&H_sorai[i]);
          /*for( i=0; i<n; i++ ) TOROL_LISTA(&Temp_sorai[i]);*/
          return;
        }
        
        printf("%d.\n", m);
  
      if( m+1 ==maxit ){
    puts("A következő iterációs lépésben véget ér a program!\n");
    printf("A legnagyobb külünbség: kezdővektor-kapottmegoldás: \n");
    puts("Szeretné folytatni az iterációt, annak érdkében, hogy jobb eredményt kapjon?\n");
    int jo=0;
    char valasz[10];
    
    while( jo == 0 ){
      puts("Parancsok: [i|n], [igen|nem], [y|n], [yes|no]");
      scanf("%s", valasz);
      if( strncmp(valasz, "i", strlen("i"))==0 || strncmp(valasz, "igen", strlen("igen"))==0 || 
          strncmp(valasz, "y", strlen("y"))==0 || strncmp(valasz, "yes", strlen("yes"))==0){
        int ujmaxit=0;
        printf("uj maxit: "); scanf("%d", &ujmaxit);
        while( ujmaxit<=maxit ){
          puts("HIBA! A megadott új maximális iterációszám kisebb, mint a jelenlegi!\n");
          printf("Adjon %d -tól/-től nagyobb számot!\n", maxit);
          printf("uj maxit: "); scanf("%d", &ujmaxit);
        }
        maxit=ujmaxit;
        jo=1;
      }
      else if( strncmp(valasz, "n", strlen("n"))==0 || strncmp(valasz, "nem", strlen("nem"))==0 || 
                strncmp(valasz, "no", strlen("no"))==0){
                jo=1;
        fflush(fout);
          Eredmeny_kiiras( out, n, x0, m );
          /*Ellenorzes( (double *)A, n, b, x0 );*/
          fflush(fout);
          fclose(fout);
          /*free(A);*/ /*free(HHT); free(Temp);*/ free(b); free(x0); free(xk); free(h0); free(hk); free(g0); free(gk); free(e0); free(ek);
          for( i=0; i<n; i++ ) TOROL_LISTA(&A_sorai[i]);
          for( i=0; i<n; i++ ) TOROL_LISTA(&H_sorai[i]);
          /*for( i=0; i<n; i++ ) TOROL_LISTA(&Temp_sorai[i]);*/
          return;
      }
      else      
        puts("HIBA! Rossz parancsot adott meg!");
                
    }
  }
    }
return ;
}
/*____________________________________________________________________________*/
void beolvas2 ( char in[], SOR ** A, SOR **H, double* b, double *x0, double *xk, int* maxit, double * e ){
  int n, i, j, nnz, k;
  
  FILE *fin;    
  fin = fopen(in, "r");
  
  fscanf(fin, "%d\n",&n);
  fscanf(fin, "%d\n",&nnz);
   
  for ( i=0; i<n; i++ ){
    A[i]=NULL;
    H[i]=NULL;
  }
  for ( k=0; k<nnz; k++ ){
    double temp;
    fscanf(fin, "%d %d %lf", &i, &j, &temp);
    BESZUR_ELEJE(&A[i], j, temp);
    if ( i!=j )
      BESZUR_ELEJE(&A[j], i, temp);
  }  
  
  rendez(A, n);
  
  fscanf(fin, "%d\n",&nnz);
  for ( k=0; k<nnz; k++ ){
    double temp;
    fscanf(fin, "%d %d %lf", &i, &j, &temp);
    BESZUR_ELEJE(&H[i], j, temp);
    if ( i!=j )
      BESZUR_ELEJE(&H[j], i, temp);
  }
  rendez(H, n);
  
    for ( i=0; i<n; i++ ){
    	fscanf(fin, "%lf", &b[i] );
    }
    for ( i=0; i<n; i++ ){
    	fscanf(fin, "%lf", &x0[i] );
    	xk[i] = x0[i];
    }
    fscanf(fin, "%d", maxit);
    fscanf(fin, "%lf", e);
}
/*____________________________________________________________________________*/
void BESZUR_ELEJE ( SOR ** f, int o, double e ){
  SOR *uj;
  uj = (SOR *) malloc (sizeof (SOR));
  uj->oszlop=o;
  uj->ertek=e;
    
  uj->kov=*f;
  *f=uj;
  return ;
}
/*______________________________________________________________________________________*/
void rendez( SOR ** f, int n){
  int i;
  
  for( i=0; i<n; i++){
    SOR *x=f[i];
    int sorhossz;
    for( sorhossz=0; x!=NULL; sorhossz++ ) x=x->kov;
    if ( sorhossz > 1 && sorhossz>0){
      int *rendezett_oszlop = (int *) calloc( sorhossz, ( sizeof ( int ) ) );
      double *rendezett_ertek = (double *) calloc( sorhossz, ( sizeof ( double ) ) );
      for( x=f[i], sorhossz=0; x!=NULL; sorhossz++, x=x->kov ){
        rendezett_ertek[sorhossz]=x->ertek;
        rendezett_oszlop[sorhossz]=x->oszlop;
      }
      Beszuro_rendezes(rendezett_ertek, rendezett_oszlop, sorhossz);
      for( x=f[i], sorhossz=0; x!=NULL; sorhossz++, x=x->kov ){
        x->ertek=rendezett_ertek[sorhossz];
        x->oszlop=rendezett_oszlop[sorhossz];
      }
      free(rendezett_ertek);
      free(rendezett_oszlop);
    }    
  }
}
/*______________________________________________________________________________________*/
void Beszuro_rendezes(double *A, int *oszlop, int hossz){
  if(hossz<2) return; /* nem kell rendezni */
  int i,j;
        
  double kulcs_e;
  int kulcs_o;
  for( j=1; j<hossz; j++ ){
    kulcs_e=A[j];
    kulcs_o=oszlop[j];
    i=j-1;
    /*puts("while");*/
    while ( i>=0 && oszlop[i]>kulcs_o ){      
      A[i+1]=A[i];
      oszlop[i+1]=oszlop[i];
      i=i-1;
    }
    A[i+1]=kulcs_e;
    oszlop[i+1]=kulcs_o;
  }
}
/*____________________________________________________________________________*/
void beolvas ( char in[], SOR ** A, SOR ** H,/* SOR **Temp,*/ double* b, double *x0, double *xk, int* maxit, double * e ){
  int n, i, j;
  
  FILE *fin;    
  fin = fopen(in, "r");
  
  fscanf(fin, "%d\n",&n);
  
  for ( i=0; i<n; i++ ){
      A[i] = NULL;
      SOR *tail = NULL;
      for ( j=0; j<n; j++ ){
        double temp;
        fscanf(fin, "%lf", &temp);
        if ( fabs(temp) > pow(10,-15)){
          SOR *uj;
          uj = (SOR *) malloc (sizeof (SOR));
          uj->oszlop=j;
          uj->ertek=temp;            
          uj->kov=NULL;
          
          if ( tail == NULL ){
            A[i] = uj;
            tail = A[i];
          }
          else{
            tail->kov=uj;
            tail = uj;
          }
        }
      }
    }
    
    for ( i=0; i<n; i++ ){
      H[i] = NULL;
      /*Temp[i] = NULL;*/
      SOR *tail = NULL;
      for ( j=0; j<n; j++ ){
        double temp;
        fscanf(fin, "%lf", &temp);
        if ( fabs(temp) > pow(10,-15)){
          SOR *uj;
          uj = (SOR *) malloc (sizeof (SOR));
          uj->oszlop=j;
          uj->ertek=temp;            
          uj->kov=NULL;
          
          if ( tail == NULL ){
            H[i] = uj;
            tail = H[i];
          }
          else{
            tail->kov=uj;
            tail = uj;
          }
        }
      }
    }
    
    for ( i=0; i<n; i++ ){
    	fscanf(fin, "%lf", &b[i] );
    }
    for ( i=0; i<n; i++ ){
    	fscanf(fin, "%lf", &x0[i] );
    	xk[i] = x0[i];    }
    fscanf(fin, "%d", maxit);
    fscanf(fin, "%lf", e);
 
 fclose(fin);
}
/*____________________________________________________________________________*/
void BESZUR_VEGE ( SOR ** f, int o, double e ){
  SOR *uj,*x;
  uj = (SOR *) malloc (sizeof (SOR));
  uj->oszlop=o;
  uj->ertek=e;
    
  uj->kov=NULL;
  if ( *f==NULL ){ *f=uj; }
  else
  {
     x=*f;
     while ( x->kov!=NULL )
           x=x->kov;
     x->kov=uj;
  }
  return ;
}
/*____________________________________________________________________________*/
double ELEM ( SOR ** f, int j ){
  if( f==NULL ) return 0.0;
  
  SOR *x;
  x = *f;
  if( x->oszlop > j ) return 0.0;
  while ( x != NULL ){
    if ( x->oszlop == j )
      return x->ertek;
    else if( x->oszlop >j )
      return 0.0;
    x=x->kov;
  }
  return 0.0;
}
/*____________________________________________________________________________*/
void G_0( SOR **sorok, double *x0, double *b, int n, double *eredmeny1, double *eredmeny2){
  int i, j;
  
  for(i=0; i<n; i++){
    SOR *x=*(sorok+i);
    double sum=0.0;
    for ( j=0; j<n && x!=NULL; j++ ){
      if ( x->oszlop==j /*&& x->ertek!=0.0*/ && x0[j]!=0.0 ){
        sum += x->ertek * x0[j];
        x=x->kov;
      }
      else if ( x->oszlop==j )
        x=x->kov;
    }
    eredmeny1[i] = eredmeny2[i] = sum-b[i];
  }
  return;
}
/*____________________________________________________________________________*/
void HxHT ( SOR ** H, SOR ** eredmeny, int n ){
  int i, j;
  
  for ( i=0; i<n; i++ ){
    *(eredmeny+i)=NULL;
    for ( j=0; j<n; j++ ){
       SOR *y=*(H+j);
       SOR *x=*(H+i);
      double sum = 0;
      for( ; y!=NULL && x!=NULL; ){
          /*printf("%lf += %lf*%lf\n", sum, x->ertek, y->ertek);*/
          if ( x->oszlop == y->oszlop ){
            sum += x->ertek * y->ertek;
          }
          
          if ( x->oszlop > y->oszlop )
            y=y->kov;
          else if ( x->oszlop < y->oszlop )
            x=x->kov;
          else{
            x=x->kov;
            y=y->kov;
          }
      }
      /*puts("---");*/
      if ( fabs(sum) > pow(10,-15) )
        BESZUR_VEGE( (eredmeny+i), j, sum );
    }    
  }
  return;
}
/*____________________________________________________________________________*/
double LAMBDA_k( SOR **sorok, double *ek, double *gk, double *hk, int n){
  double szamlalo=0.0, lambda=0.0;
  int i, j;
  
  for ( i=0; i<n; i++ )
    szamlalo += ek[i] * gk[i];
  for(i=0; i<n; i++){
    SOR *x=*(sorok+i);
    double sum=0.0;
    for ( j=0; j<n && x!=NULL; j++ ){
      if ( x->oszlop==j /*&& x->ertek!=0.0*/ && hk[j]!=0.0 ){
        sum += x->ertek * hk[j];
        x=x->kov;
      }
      else if ( x->oszlop==j )
        x=x->kov;
    }
    lambda += sum * hk[i];
  }
  return szamlalo/lambda;
}
/*____________________________________________________________________________*/
void G_k( SOR **sorok, double *gk, double lambda, double *hk, int n, double *eredmeny){
  int i, j;
  
  for(i=0; i<n; i++){
    SOR *x=*(sorok+i);
    double sum=0.0;
    for ( j=0; j<n && x!=NULL; j++ ){
      if ( x->oszlop==j /*&& x->ertek!=0.0*/ && hk[j]!=0.0 ){
        sum += x->ertek * hk[j];
        x=x->kov;
      }
      else if ( x->oszlop==j )
        x=x->kov;
    }
    eredmeny[i] = sum*lambda+gk[i];
  }
  return;
}
/*____________________________________________________________________________*/
void TOROL_LISTA ( SOR ** f ){
  SOR *temp;
  
  while ( *f != NULL)
  {  
     temp = (*f) -> kov;
     free (*f);
     *f = temp;
  }
  return;
}
/*____________________________________________________________________________*/
void Eredmeny_kiiras( char out[], int n, double * vektor, int it ){
  
  FILE *fout;

  fout = fopen(out, "a");
  
  int i;
  for ( i=0; i<n; i++ ){
    if (i<n-1)
      fprintf(fout, "%.5lf ", vektor[i] );
     else
      fprintf(fout, "%.5lf", vektor[i] );
  }
  fputs("", fout);
  printf("A legnagyobb külünbség: ||kezdővektor-kapottmegoldás||= %.16lf\n", Ellenorzes_legnagyobb_kulonbseg(vektor, n) );
  fprintf(fout, "\n\nA legnagyobb külünbség: ||kezdővektor-kapottmegoldás||= %.16lf\n", Ellenorzes_legnagyobb_kulonbseg(vektor, n) );
  printf("Lépések száma: %d", it);
  fprintf(fout, "Lépések száma: %d", it);
  fflush(fout);
  fclose(fout);
return;
}
/*____________________________________________________________________________*/
double Ellenorzes_legnagyobb_kulonbseg( double * x, int n){
  double *temp=(double *) calloc( n, ( sizeof ( double ) ) );
  int i=0;
  for( ; i<n; i++ )
    temp[i] = 1-x[i];
  
  return sqrt( Euklideszi_norma_vektorra(temp, n) );
  
}
/*____________________________________________________________________________*/
double Euklideszi_norma_vektorra( double * vektor, int o ){
       double temp=0.0;
       int szamlal;
       for ( szamlal=0; szamlal<o; szamlal++ ){
           temp = temp + (vektor[szamlal]*vektor[szamlal]);
       }
       return temp;
}
/*____________________________________________________________________________*/
void lin_e_r( SOR **H, /*SOR **HT,*/ double *b, int n ){
int i,j;
  /* Hy=b */
  double * hii =(double *) calloc( n, ( sizeof ( double ) ) );
  for( i=0; i<n; i++ ){
    SOR *x=*(H+i);
    for( j=0; j<=i-1 && x!=NULL; j++ ){
      if ( x->oszlop==j && x->ertek!=0 && b[j]!=0.0 ){
        b[i] = b[i] - x->ertek * b[j];
        x=x->kov;
      }
      else if ( x->oszlop==j ){        
        x=x->kov;
      }
    }
    if( i==j ) hii[i] = x->ertek;
    if ( fabs( hii[i] ) <= pow ( 10,-15 ) ){
      hiba = 1;
      free(hii);
      return;
    }
    b[i] = b[i]/hii[i];
  }
  
  /* HTx=y */  
  for(j=n-1;j>=0;j--){
    SOR *x=*(H+j);
    if ( fabs( hii[j] ) <= pow ( 10,-15 ) ){
      hiba = 1;
      free(hii);
      return;
    }
    /*b[j]=b[j]/ELEM(&H[j],j);*/
    b[j]=b[j]/hii[j];
   for(i=0;i<j&&x!=NULL;i++){
    if ( x->oszlop==i && x->ertek!=0 && b[i]!=0.0 ){
      b[i]=b[i]-x->ertek*b[j];
      x=x->kov;
    }
    else if ( x->oszlop==i ){
      x=x->kov;
    }
   }
  }
  free(hii);
  return;
}

/*____________________________________________________________________________*/
void LDLT( SOR **sorok/*double *A*/, double *b, int n, double *eredmeny, int hiba ){
int i, j, k;
    i = j = k = 0;

    double *x=(double *) calloc( n, ( sizeof ( double ) ) );
    /*if ( szimmetrikus( (double *)Atemp, n ) ){*/
         /* Az algoritmus:*/
	   	 double *r=(double *) calloc( n, ( sizeof ( double ) ) );
	   	 double *d=(double *) calloc( n, ( sizeof ( double ) ) );
	   	 
	   	 for ( j=0; j<n; j++ ){
          /*for ( i=j; i<n; i++ )
            r[i] = ELEM(&sorok[i], j);*/ /*Atemp[i][j]*/              
          
          R_i_first( sorok, j, n, r );
          
          x[j]=b[j];
          for ( k=0; k<=j-1; k++ ){
 		  	 	 for ( i=j; i<n; i++ ) r[i] = r[i]-ELEM(&sorok[i], k)*r[k]*ELEM(&sorok[j], k);
 		  	 	 /*R_i( sorok, j, k, n, r);*/
 		  	 	 x[j] = x[j]-ELEM(&sorok[j], k)*x[k];
         }
         
		     if ( fabs( r[j] ) <= pow ( 10,-15 ) ){
          hiba = 1;
            for ( i=0; i<n; i++ )
              eredmeny[i] = x[i];
            free(x); free(r); free(d);
            return;
		      }
		     else d[j]=1/r[j];
		     for ( i=j+1; i<n; i++ ){
		 	 	 /*Atemp[i][j] = r[i]*d[j];*/
		 	 	 BESZUR ( sorok, i, j, r[i]*d[j] );
		  	 }
		 }
	   	 for ( i=n-1; i>=0; i-- ){ /*i<1 ről módosult*/
         x[i]=d[i]*x[i];
         for ( k=i+1; k<n; k++ ){ x[i]=x[i]-ELEM(&sorok[k], i)*x[k]; }
         /*X_j( sorok,  i, n, x);*/
		 }
    /*}*/
    
    for ( i=0; i<n; i++ ){
      eredmeny[i] = x[i];
    }
    
free(x); free(r); free(d);
return;
}
/*____________________________________________________________________________*/
void R_i_first( SOR **sorok, int j, int n, double *eredmeny){
  int k; 
  SOR *x=*(sorok+j);
  
  for ( k=0; k!=j && k<n && x!=NULL; k++ ){
    if ( k== x->oszlop )
      x=x->kov;
  }
  for ( ; k<n && x!=NULL; k++ ){
    if ( x->oszlop==k ){
      eredmeny[k] = x->ertek;
      x=x->kov;
    }
    else
      eredmeny[k] = 0.0;
  }
  return;
}
/*____________________________________________________________________________*/
void R_i( SOR **sorok, int j, int k, int n, double *eredmeny){
  
  SOR *x=*(sorok+k); /* ITT MEG LETT FORDÍTVA lik volt és én  a szimmetria miatt lki -t használom, mert így 1 soron tudok végig menni */
  double ljk = ELEM(&sorok[j], k);
  int v;

  for ( v=0; v!=j && v<n && x!=NULL; v++ ){
    if ( v== x->oszlop )
      x=x->kov;
  }
  
  for( v=0; v<n && x!=NULL; v++ ){
    if ( v==x->oszlop ){
      eredmeny[v] -= x->ertek * eredmeny[k] * ljk;
      x=x->kov;
    }
  }
  return;
}
/*____________________________________________________________________________*/
void X_j( SOR **sorok,  int i, int n, double *eredmeny){
  /*for ( k=i+1; k<n; k++ ){ x[i]=x[i]-ELEM(&sorok[k], i)*x[k]; }*/
  
  SOR *x=*(sorok+i);  /* ITT MEG LETT FORDÍTVA lki volt és én  a szimmetria miatt lik -t használom, mert így 1 soron tudok végig menni */
  int v;
  
  for ( v=0; v!=i+1 && v<n && x!=NULL; v++ ){
    if ( v== x->oszlop )
      x=x->kov;
  }
  
  for( v=0; v<n && x!=NULL; v++ ){
    if ( v==x->oszlop ){
      eredmeny[v] -= x->ertek * eredmeny[v];
      x=x->kov;
    }
  }
  return;
}
/*____________________________________________________________________________*/
void BESZUR ( SOR ** f, int s, int o, double e ){/*ha létezik, akkor fellírja*/
/*puts("BESZUR FGV");*/
  SOR *x=*(f+s);
  /*puts("1. while elott");*/
  for( ;x!=NULL; x=x->kov){
    /*puts("whileban");*/
    if (x->oszlop > o){
      /*puts("if");*/
      break;
    }
    else if (x->oszlop == o){
      /*puts("else-> már benne");*/
      x->ertek = e;
      return;
    }    
  }
  
  /*puts("nincs benne-> lefoglalok");*/
  
  SOR *uj, *elozo;
  uj = (SOR *) malloc (sizeof (SOR));
  uj->oszlop=o;
  uj->ertek=e;
  
  x=*(f+s);
  elozo=NULL;
  /*puts("2.while elott");*/
  while ( x!=NULL && x->oszlop<o )     /* < : 1, 2, 3, ...      > : ..., 3, 2, 1 */
  {
       elozo=x;
       x=x->kov;
  }
  /*puts("2.while után");*/
  uj->kov=x;
  if (elozo==NULL) *(f+s)=uj;
  else elozo->kov=uj;
  
  return ;
}
/*____________________________________________________________________________*/
int szimmetrikus ( double *A, int n){
    int i,j;
    i = j = 0;
    
    for ( i=0; i<n; i++ ){
        for ( j=0; j<n; j++ ){
          if ( A[i*n+j] != A [i*n+j] )
            return 0;
        }
    }
    return 1;
}

/*____________________________________________________________________________*/
void Ellenorzes(double *A, int n, double * b, double * x){
  int i, j;
  double temp[n];
  puts("\n\nEllenorzes: Ax=b");    
  puts("b:\t\t\tAx:");
    for ( i=0; i<n; i++ ){
      temp[i]=0;
      for ( j=0; j<n; j++ ){
        temp[i] += A[i*n+j]*x[j];
      }
      printf( "b[%d]=%lf\t\t%lf\t\t\n",i, b[i], temp[i] );
    }
}
