/**
   * Az A*x=b egyenletrendszer megoldasahoz az LDLT felbontas.
   * Ehhez a felbontashoz az szükséges, hogy A n*n -es matrix
   * szimmetrikus és pozitív definit legyen.
   */
   /* Megoldás a könyvbeni pszeudo kóddal. */
   /*
   * Bemenet:
   * n
   * A
   * b
   */
typedef struct sor
{
       int oszlop;
       double ertek;
       struct sor *kov;
} SOR;

void LDLT(char in[], char out[]);
void beolvas2 ( char in[], SOR ** A, double* b );
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 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 );

void BESZUR_VEGE ( SOR ** f, int o, double e );
double ELEM ( SOR ** f, int j );
void TOROL_LISTA ( SOR ** f );

int szimmetrikus ( double *A, int n);
int pozitiv_definit ( double *A, int n);
int det ( double *A, int n);
void Ellenorzes(double *A, int n, double * b, double * x);
double Ellenorzes_legnagyobb_kulonbseg( double * x, int n);
double Euklideszi_norma_vektorra( double * vektor, int o );

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

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);
    
    LDLT(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 LDLT(char in[], char out[]){
    int n, i, j, k;
    n = i = j = k = 0;
    
    FILE *fin, *fout;
    
    fin = fopen(in, "r");
    fout = fopen(out, "w");
         
    fscanf(fin, "%d\n", &n);
    double *b = (double *) calloc( n, ( sizeof ( double ) ) );
    SOR *A[n];
    
    fclose(fin);
    beolvas2 ( in, A, b );

    	
    /*if ( szimmetrikus( (double *)A, n ) ){*/
         /* Az algoritmus:*/
	   	 double *r = (double *) calloc( n, ( sizeof ( double ) ) );
	   	 double *x = (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(&A[i], j); }
         /*R_i_first( A, 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(&A[i], k)*r[k]*ELEM(&A[j], k); }
             x[j] = x[j]-ELEM(&A[j], k)*x[k];	 
           
         }
           if ( fabs( r[j] ) <= pow ( 10,-15 ) ){
            fputs("hiba", fout);
            puts("hiba");
            fprintf(fout, "j=%lf",r[j]);
            fflush(fout);
            fclose(fout);
            free(b); free(r); free(x); free(d);
            for( i=0; i<n; i++ ) TOROL_LISTA(&A[i]);
            return;
            }
           else d[j]=1/r[j];
           for ( i=j+1; i<n; i++ ){
            /*Atemp[i][j] = r[i]*d[j];*/
            BESZUR ( A, i, j, r[i]*d[j] );
           }
		 }
	   	 for ( i=n-1; i>=0; i-- ){ 
         x[i]=d[i]*x[i];
         for ( k=i+1; k<n; k++ ){ x[i]=x[i]-ELEM(&A[k], i)*x[k]; }
		 }
		 for ( i=0; i<n; i++ ){ 
        if (i<n-1) fprintf(fout, "%.5lf ", x[i] );
        else fprintf(fout, "%.5lf", x[i] );
      }
    fputs("", fout);
    printf("A legnagyobb külünbség: ||kezdővektor-kapottmegoldás||= %.16lf\n", Ellenorzes_legnagyobb_kulonbseg(x, n) );
    fprintf(fout, "\n\nA legnagyobb külünbség: ||kezdővektor-kapottmegoldás||= %.16lf\n", Ellenorzes_legnagyobb_kulonbseg(x, n) );
    /*Ellenorzes( (double *)Atemp, n, b, x );*/
    free(r); free(x); free(d);
    /*}*/
    
    fflush(fout);
    fclose(fout);
    free(b);
    for( i=0; i<n; i++ ) TOROL_LISTA(&A[i]);
return;
}
/*____________________________________________________________________________*/
void beolvas2 ( char in[], SOR ** A, double* b ){
  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;
  }
  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);
  
  for ( i=0; i<n; i++ ){
    fscanf(fin, "%lf", &b[i] );
  }
}
/*____________________________________________________________________________*/
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 TOROL_LISTA ( SOR ** f ){
  SOR *temp;
  
  while ( *f != NULL)
  {  
     temp = (*f) -> kov;
     free (*f);
     *f = temp;
  }
  return;
}
/*____________________________________________________________________________*/
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 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("\n%lf += %lf*%lf ", temp[i], A[i*n+j], x[j]);*/
      }
      printf( "b[%d]=%lf\t\t%lf\t\t\n",i, b[i], temp[i] );
    }
}
/*____________________________________________________________________________*/
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;
}
