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

void Prekondicionalt_Konjugalt_gradiens(char in[], char out[]);
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 LDLT( double *A, double *b, int n, double *eredmeny, int hiba );
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);

#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 A[n][n], H[n][n], HT[n][n], Temp[n][n], b[n], x0[n], xk[n], h0[n], hk[n], g0[n], gk[n], mk, lk, e;*/    /* lk == lambda k */
    double (*A)[n] = NULL;
    A = calloc (n, sizeof (*A));
    
    double (*HHT)[n] = NULL;
    HHT = calloc (n, sizeof (*HHT));
    
    double (*Temp)[n] = NULL;
    Temp = calloc (n, sizeof (*Temp)); 
    
    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;
    
    /*----------------------------------------*/
    /* A matrix beolvasasa */
    for ( i=0; i<n; i++ ){
        for ( j=0; j<n; j++ ){
            fscanf(fin, "%lf", &A[i][j]);
        }
    }
    /* H mátrix beolvasása */
    for ( i=0; i<n; i++ ){
        for ( j=0; j<n; j++ ){
            fscanf(fin, "%lf", &HHT[i][j]);
        }
    }
    /* b vektor beolvasása */
    for ( i=0; i<n; i++ ){
    	fscanf(fin, "%lf", &b[i] );
    	/*if ( i < n-2 ) scanf(" ");
	    if ( i == n-2 ) scanf("\n");*/
    }
    /* x0 vektor beolvasása */
    for ( i=0; i<n; i++ ){
    	fscanf(fin, "%lf", &x0[i] );
    	xk[i] = x0[i];
    	/*if ( i < n-2 ) scanf(" ");
	    if ( i == n-2 ) scanf("\n");*/
    }
    fscanf(fin, "%d",&maxit);
    fscanf(fin, "%lf",&e);
    fclose(fin);
    
    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*/
	for ( i=0; i<n; i++ ){
    g0[i]=0;
    for ( j=0; j<n; j++ ){
        g0[i] += A[i][j] * x0[j];
    }
	}
	
	for ( i=0; i<n; i++ ){
    g0[i] = g0[i] - b[i];
    gk[i] = g0[i];
   }
	
	/* HHT e0 = g0 */ /* e0 az egy vektor, ek+1 az ek    ek az e0 */
	/* TODO e0 meghatározása.*/
  /* Temp = H*HT */
	for ( i=0; i<n; i++ ){
    for ( j=0; j<n; j++ ){
      Temp[i][j] = 0;
      for( k=0; k<n; k++ )
        Temp [i][j] += HHT[i][k] * HHT[j][k];     /* Itt Temp [i][j] += HHT[i][k] * HHT[k][j];  volt */
    }
  }
  
  /* e0 megh LDLT seg.vel */
  hiba = 0;
  LDLT( (double *)Temp, g0, n, e0, hiba );
  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);
      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) */
          for ( i=0; i<n; i++ ){
            temp[i]=0;
            for ( j=0; j<n; j++ ){
              temp[i] += A[i][j] * h0[j];
            }
          }
                    
          for ( i=0; i<n; i++ )
            nevezo += temp[i] * h0[i];
          for ( i=0; i<n; i++ )
            szamlalo += e0[i] * g0[i];
          lk = szamlalo / nevezo;
          
          /*printf("lk: %lf\n", lk);*/
          
          /*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*/
          for ( i=0; i<n; i++ ){
            temp[i]=0;
            for ( j=0; j<n; j++ ){
              temp[i] += A[i][j] * h0[j];
            }
          }
          for ( i=0; i<n; i++ )
            temp[i] = lk*temp[i];
          for ( i=0; i<n; i++ )
            gk[i] = g0[i] + temp[i];
           
          /* HHT ek+1 = gk+1 */
          /* TODO ek+1 meghatározása.*/
          for ( i=0; i<n; i++ ){
            for ( j=0; j<n; j++ ){
              Temp[i][j] = 0;
              for( k=0; k<n; k++ )
                Temp [i][j] += HHT[i][k] * HHT[j][k];     /* Itt Temp [i][j] += HHT[i][k] * HHT[k][j];  volt */
            }
          }
          
          /* ek+1 megh LDLT seg.vel */
          /*ek = LDLT( (double *)Temp, gk, n );*/
          free(ek);
          ek=(double *) calloc( n, ( sizeof ( double ) ) );
          hiba = 0;
          LDLT( (double *)Temp, gk, n, ek, hiba );
          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);*/
            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);
          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);
          return;
      }
      else      
        puts("HIBA! Rossz parancsot adott meg!");
                
    }
  }
        
    }
    /*puts("vmi");
    for ( i=0; i<n; i++ ){ printf( "%.5lf ", x0[i] ); }
    printf("\nLepesek szama: %d",m);
    printf("A legnagyobb külünbség: ||kezdővektor-kapottmegoldás||= %lf\n", Ellenorzes_legnagyobb_kulonbseg(x0, n) );*/
    /*Ellenorzes( (double *)A, n, b, x0 );*/
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 LDLT( double *A, double *b, int n, double *eredmeny, int hiba ){
int i, j, k;
    i = j = k = 0;
         
    /*double Atemp[n][n];*/
    double (*Atemp)[n] = NULL;
    Atemp = calloc (n, sizeof (*Atemp));
    
    /* A matrix beolvasasa */
    for ( i=0; i<n; i++ ){
        for ( j=0; j<n; j++ ){
            Atemp [i][j] = A[i*n+j];
        }
    }

    /*double x[n];*/
    double *x=(double *) calloc( n, ( sizeof ( double ) ) );
    if ( szimmetrikus( (double *)A, n ) ){    
         /* Az algoritmus:*/
	   	 /*double r[n], d[n];*/
	   	 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] = A[i*n+j]; }
          x[j]=b[j];
          for ( k=0; k<=j-1; k++ ){
 		  	 	 for ( i=j; i<n; i++ ){ r[i] = r[i]-A[i*n+k]*r[k]*A[j*n+k]; } /*!!!!*/
 		  	 	 x[j] = x[j]-A[j*n+k]*x[k];
         }
		     if ( fabs( r[j] ) <= pow ( 10,-15 ) ){
          /*puts("hiba"); printf("j=%lf\n",r[j]);*/
          hiba = 1;
            for ( i=0; i<n; i++ )
              eredmeny[i] = x[i];
            return;
		      }
		     else d[j]=1/r[j];
		     for ( i=j+1; i<n; i++ ){
		 	 	 A[i*n+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]-A[k*n+i]*x[k]; }/*k=i ről módosult*/
		 }
    }
    
    for ( i=0; i<n; i++ ){
      eredmeny[i] = x[i];
    }

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] );
    }
}
