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

void LDLT(char in[], char out[]);
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 (*A)[n] = NULL;
    A = calloc (n, sizeof (*A));
    
    double (*Atemp)[n] = NULL;
    Atemp = calloc (n, sizeof (*Atemp));
    
    /*double b[n];*/
    double *b = (double *) calloc( n, ( sizeof ( double ) ) );
    
    /* A matrix beolvasasa */
    for ( i=0; i<n; i++ ){
        for ( j=0; j<n; j++ ){
            fscanf(fin, "%lf", &A[i][j]);
            Atemp[i][j] = A[i][j];
        }
    }
    /* b vektor beolvasása */
    for ( i=0; i<n; i++ ){
    	fscanf(fin, "%lf", &b[i] );
    }
    fclose(fin);
    
    
    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] = A[i][j]; }
         x[j]=b[j];
         for ( k=0; k<=j-1; k++ ){
             for ( i=j; i<n; i++ ){ r[i] = r[i]-A[i][k]*r[k]*A[j][k]; }
             x[j] = x[j]-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(A); free(Atemp); free(b); free(r); free(x); free(d);
            return;
            }
           else d[j]=1/r[j];
           for ( i=j+1; i<n; i++ ){
            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]-A[k][i]*x[k]; }
		 }
		 for ( i=0; i<n; i++ ){ 
        if (i<n-1){
          fprintf(fout, "%.5lf ", x[i] );
          printf("%.5lf ", x[i] );
        }
        else{
          fprintf(fout, "%.5lf", x[i] );
          printf("%.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(A); free(Atemp); free(b);
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;
}
