/**
   * Az Ax=b egyenletrendszer megoldasahoz a konjugált gradiens módszer.
   * x* az Ax*=b megoldása.
   */
   
   /*
   * Bemenet:
   * n
   * A
   * b
   * x0
   * maxit
   * e
   */
   
   /* Ahol
   * n       : nxn -es mátrix.
   * A       : 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 Konjugalt_gradiens(char in[], char out[]);

void beolvas ( char in[], SOR ** A, double* b, double *x0, double *xk, int* maxit, double * e );
void Beszuro_rendezes(double *A, int *oszlop, int hossz);
void rendez( SOR ** f, int n);
void BESZUR_ELEJE ( SOR ** f, int o, double e );
void beolvas2 ( char in[], SOR ** A, 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 *g0, double *gk, double *d0, double *dk );
void G_k( SOR **sorok, double *g0, double lk, double *d0, int n, double *gk );
double LAMDA_k( SOR ** sorok, double *d0, double norma2, int n );
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 Ellenorzes(double *A, int n, double * b, double * x);
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <sys\timeb.h>
/*#include "leak_detector_c.h"*/

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);
    
    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 Konjugalt_gradiens(char in[], char out[]){
    int n, i, j, k, maxit;
    n = i = j = k = maxit = 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 *b = (double *) calloc( n, ( sizeof ( double ) ) );
    double *x0 = (double *) calloc( n, ( sizeof ( double ) ) );
    double *xk = (double *) calloc( n, ( sizeof ( double ) ) );
    double *g0 = (double *) calloc( n, ( sizeof ( double ) ) );
    double *gk = (double *) calloc( n, ( sizeof ( double ) ) );
    double *d0 = (double *) calloc( n, ( sizeof ( double ) ) );
    double *dk = (double *) calloc( n, ( sizeof ( double ) ) );
    
    double mk, lk, e;
    SOR *A_sorai[n];

    fclose(fin);
    beolvas2 ( in, A_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 */
	
	/*printf("%d\n", n);
	puts("A mátrix:");
	for(i=0;i<n;i++){
    SOR *x = A_sorai[i];
      for( ; x!=NULL; x=x->kov){
        printf("%d %d: %lf\n", i,x->oszlop,x->ertek);
    }
  }
  
  puts("b:");
  for(i=0;i<n;i++)
    printf("%d: %lf\n", i, b[i]);
  puts("x0:");
  for(i=0;i<n;i++)
    printf("%d: %lf\n", i, x0[i]);
  printf("maxit: %d\n", maxit);
  printf("e: %lf\n", e);*/
  
  
	G_0( A_sorai, x0, b, n, g0, gk, d0, dk );
	/* g0=A*x0-b, d0=g0 */
	/*for ( i=0; i<n; i++ ){
    g0[i]=0;
    for ( j=0; j<n; j++ ){
        g0[i] += ELEM(&A_sorai[i],j) * x0[j];
    }
	}
	for ( i=0; i<n; i++ ){
    g0[i] = g0[i] - b[i];
    gk[i] = g0[i];
    d0[i] = g0[i]*(-1);
    dk[i] = d0[i];
   }*/
	/*puts("g0:");
  for(i=0;i<n;i++)
    printf("%d: %lf\n", i, g0[i]);
  return;*/
    
	double *temp = (double *) calloc( n, ( sizeof ( double ) ) );
	for ( k=0; k<maxit; k++ ){
        /* [?||gk||<=e?[stop, eredmény: xk] */
        double norma = Euklideszi_norma_vektorra ( g0, n );
        norma = sqrt (norma);
        if ( fabs( norma ) >= e ){
          double nevezo=0.0;
          /*double temp[n];*/
          
          /* lk=||gk||2/(Adk,dk) */
          norma = Euklideszi_norma_vektorra ( g0, n );
          norma = sqrt (norma);
          norma = norma * norma;
          lk = LAMDA_k( A_sorai, d0, norma, n );
          /*for ( i=0; i<n; i++ ){
            temp[i]=0;
            for ( j=0; j<n; j++ ){
              temp[i] += ELEM(&A_sorai[i],j) * d0[j];
            }
          }
          for ( i=0; i<n; i++ )
            nevezo += temp[i] * d0[i];
          lk = norma / nevezo;*/
          
          /*printf("lk: %lf\n", lk);*/
          
          /*xk+1=xk+lkdk*/ /*xk+1 az xk    xk az x0*/
          for ( i=0; i<n; i++ )
            temp[i] = lk*d0[i];
          for ( i=0; i<n; i++ )
            xk[i] = x0[i] + temp[i];
          
          /*gk+1=gk+lkAdk*/
          G_k( A_sorai, g0, lk, d0, n, gk );
          /*for ( i=0; i<n; i++ ){
            temp[i]=0;
            for ( j=0; j<n; j++ ){
              temp[i] += ELEM(&A_sorai[i],j) * d0[j];
            }
          }
          for ( i=0; i<n; i++ )
            temp[i] = lk*temp[i];
          for ( i=0; i<n; i++ )
            gk[i] = g0[i] + temp[i];*/
          
          /*puts("xk+1:");
          for(i=0;i<n;i++)
            printf("%d: %lf\n", i, xk[i]);
          puts("gk+1:");
          for(i=0;i<n;i++)
            printf("%d: %lf\n", i, gk[i]);*/
          
          /*mk=||gk+1||2/||gk||2*/    /* gk+1 az gk   gk az g0*/
          nevezo = norma = Euklideszi_norma_vektorra ( g0, n );
          nevezo = sqrt(nevezo);
          nevezo = nevezo * nevezo;
          norma = Euklideszi_norma_vektorra ( gk, n );
          norma = sqrt (norma);
          norma = norma * norma;
          mk = norma / nevezo;
          /*dk+1=-gk+1 + mkdk*/   /* dk+1 az dk   dk az d0*/
          for ( i=0; i<n; i++ )
            temp[i] = mk * d0[i];
          for ( i=0; i<n; i++ )
            dk[i] = (gk[i] * (-1)) + temp[i];          
        }
        else {
          fflush(fout);
          Eredmeny_kiiras( out, n, x0, k );
          /*Ellenorzes( (double *)A, n, b, x0 );*/
          fflush(fout);
          fclose(fout);
          /*free(A);*/   free(b);   free(x0);   free(xk);   free(g0);   free(gk);   free(d0);   free(dk);
          for( i=0; i<n; i++ ) TOROL_LISTA(&A_sorai[i]);
          return;
        }
        for ( i=0; i<n; i++ ){
          x0[i] = xk[i];
          g0[i] = gk[i];
          d0[i] = dk[i];
        }
        
        printf("%d.\n", k);
  
        if( k+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;
      else      
        puts("HIBA! Rossz parancsot adott meg!");
                
    }
  }
        
    }
    free(temp);
    
    fflush(fout);
    Eredmeny_kiiras( out, n, x0, k );
    /*Ellenorzes( (double *)A, n, b, x0 );*/
    fflush(fout);
    fclose(fout);
    /*free(A);*/   free(b);   free(x0);   free(xk);   free(g0);   free(gk);   free(d0);   free(dk);
    for( i=0; i<n; i++ ) TOROL_LISTA(&A_sorai[i]);
return ;
}
/*____________________________________________________________________________*/
void beolvas2 ( char in[], SOR ** A, 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;    
  }
  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++){
    SOR *x=A[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=A[i], sorhossz=0; x!=NULL; sorhossz++, x=x->kov ){
        rendezett_ertek[sorhossz]=x->ertek;
        rendezett_oszlop[sorhossz]=x->oszlop;
      }
      printf("%d.sor:", i);
      Beszuro_rendezes(rendezett_ertek, rendezett_oszlop, sorhossz);
      for( x=A[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);
    }    
  }*/
  
    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 ){
/*puts("BESZUR_ELEJE");*/
  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;
      }
      /*printf("%d.sor:", i);*/
      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){
/*puts("\nBeszuro_rendezes");*/
  if(hossz<2) return; /* nem kell rendezni */
  int i,j;
  
  /*puts("rendezetlenul:");
  for( j=0; j<hossz; j++ )
    printf("%d.:  %d # %lf\n", j, oszlop[j], A[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;
      /*int k;
      for( k=0; k<hossz; k++ )
        printf("%d.:  %d # %lf\n", k, oszlop[k], A[k]);*/
    }
    A[i+1]=kulcs_e;
    oszlop[i+1]=kulcs_o;
  }
  /*puts("rendezve:");
  for( j=0; j<hossz; j++ )
    printf("%d.:  %d # %lf\n", j, oszlop[j], A[j]);
  puts("Beszuro_rendezes vege");*/
}
/*____________________________________________________________________________*/
void beolvas ( char in[], SOR ** A, 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++ ){
    	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_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 ;
}
/*____________________________________________________________________________*/
void TOROL_LISTA ( SOR ** f ){
  SOR *temp;
  
  while ( *f != NULL)
  {  
     temp = (*f) -> kov;
     free (*f);
     *f = temp;
  }
  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 *g0, double *gk, double *d0, double *dk ){
  int i, j;    
  /*g0=A*x0-b*/
  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;
    }
    g0[i] = gk[i] = sum-b[i];
  }
  for ( i=0; i<n; i++ ){
    /*g0[i] = g0[i] - b[i];
    gk[i] = g0[i];*/
    d0[i] = g0[i]*(-1);
    dk[i] = d0[i];
   }
  return;
}
/*____________________________________________________________________________*/
void G_k( SOR **sorok, double *g0, double lk, double *d0, int n, double *gk ){
  int i, j;
  /*gk+1=gk+lkAdk*/
  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 && d0[j]!=0.0 ){
        sum += x->ertek * d0[j];
        x=x->kov;
      }
      else if ( x->oszlop==j )
        x=x->kov;
    }
    gk[i] = g0[i]+lk*sum;
  }
}
/*____________________________________________________________________________*/
double LAMDA_k( SOR ** sorok, double *d0, double norma2, int n ){
  int i, j;
   /* lk=||gk||2/(Adk,dk) */
   double lk=0.0;
  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 && d0[j]!=0.0 ){
        sum += x->ertek * d0[j];
        x=x->kov;
      }
      else if ( x->oszlop==j )
        x=x->kov;
    }
    lk += sum*d0[i];
  }
  return norma2 / lk;
}
/*____________________________________________________________________________*/
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 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("%lf ", A[i*n+j]);*/
      }
      printf( "b[%d]=%lf\t\t%lf\t\t\n",i, b[i], temp[i] );
    }
}
