#include <iostream>
#include <stdlib.h>
#include <cstdlib>
#include <fstream>
#include <String>
#include <Math.h>
#include <sys/timeb.h>
using namespace std;

fstream file,file2,file3,file4,file5,file6,file7,file8,file9,file10;



//deklaracja stałych:
const double D = 1.0;
const double lambda_bez = 1.0; //dla metod bezposredniej
const double lam_pos = 1.0; //dla metod posrednich
const double h = 0.1; //krok dla x
const double dt = (lam_pos*h*h)/D; //krok dla t
const double PI=3.14159265359;
const double x_max=1; //maksymalny przedzial x
const double x_min=0; //minimalny przedzial x
const double t_max=0.5; //maksymalny przedzial x
const double t_min=0; //minimalny przedzial x


//wymiary macierzy do obliczeń numerycznych:
//const int N = (int)((2*x_max/h)+1);//liczba kolumn (os x)
//const int N = 40; //liczba kolumn(x)
//const int M = (int)((t_max/dt)+1);//liczba wierszy (os t)

const int M=(t_max-t_min)/dt+1;
const int N=(x_max-x_min)/h+1;

//**************************************************************************************************
//tworzy nową macierz MxN
double **macierz(int m, int n)
{
       double **a;
       a=new double *[m];
       for(int i=0; i<m; i++)
       a[i]=new double[n];
       return a;
}

//**************************************************************************************************
//usuwa macierz
void usun_macierz(double **a, int m)
{
     for(int i=m-1; i>=0; i--)
     delete []a[i];
     delete []a;
}

//**************************************************************************************************
//zapisywanie macierzy do pliku
void zapisz(double **macierz, fstream &plik)
{

     //plik.setf(ios::scientific, ios::floatfield);
     //plik.precision(4);


             for(int j = 0; j<M; j++){
				for(int i=0;i<N;i++){

					plik <<i*h<<" "<<j*dt<<" "<<macierz[j][i] << endl;
				}
			 }

     plik.close();
}

//**************************************************************************************************
void zapiszW(double *wek, fstream &plik)
{
             for(int j = 0; j<M; j++){
					plik << j*h<<" "<< wek[j]<< "\n";
				}
     plik.close();
}

/*********************************************************************************************/
/*metoda obliczająca czas wykonywania obliczen*/

//**************************************************************************************************
//Metoda rozwiązująca analitycznie równanie rózniczkowe
void analitycznie(double **A, double *x, double *t)
{

	//wypełnienie całej macierzy:
     /*for(int i =1; i<M; i++)
	 	{
             for(int j =1; j<N-1; j++)
					{
                     A[j][i] = (1.0 - exp(-1.0*PI*PI*D*j*dt))*sin(PI*i*h);
					}

        }
     */

     double y=0;
     for (int n = 0;n<N;n++)

	{
		for (int m = 0;m < M; m++)
		{
			A[m][n] = (1.0 - exp(-1.0*PI*PI*D*m*dt))*sin(PI*n*h);

		}
	}

     //zapisanie rozwiązania do pliku:
     zapisz(A, file);
}

//**************************************************************************************************
//klasyczna metoda bezposrednia

//**************************************************************************************************
void blad(double *BLAD, double **KMB, double **A, fstream &plik)
{
 	 double temp;
     for(int i =0; i<M; i++)
     {
	  		 double max = 0.0;
             for(int j =0; j<N-1; j++){
                     temp = fabs(KMB[i][j] - A[i][j]);
                     if(temp>max){
					 			  max = temp;
					 			  }
                     }
                     BLAD[i] = max;
     }

     zapiszW(BLAD, plik);

}

//**************************************************************************************************
double norma(double *dane){
	double w = dane[0];
	for( register int i = 1; i < N; i++ ){
		if( w < fabs(dane[i]) )
			w = dane[i] ;
	}
	return w;
}

//**************************************************************************************************
//funkcja dokonująca częściowego wyboru elementów podstawowych
void szukaj_zamien(double **A, int *indeksy, int start){
	int max = start;
	for(register  int i = start+1; i < M; i++ )
		if( fabs( A[indeksy[i]][start] ) > A[indeksy[max]][start] )
			max = i;
	//sprawdzanie czy udało się znaleść odpowiedni element
	if( A[indeksy[max]][start] == 0 ){
		printf("Bladna macierz A\n");
		exit(1);
	}
	indeksy[start] = max;
	indeksy[max] = start;
}

//**************************************************************************************************
//funkcja dokonująca dekompozycji LU
void rozklad_LU(double **LU, double **A, int *indeksy){
	double wsp;
	for(int i = 0; i < N; i++ ){
		for(int j = 0; j < N; j++ ){
			LU[indeksy[i]][j] = A[indeksy[i]][j];
		}
	}
	for( int i = 0; i < N; i++ ){
		for(int j = i+1; j < N; j++ ){
			if(LU[indeksy[i]][i] == 0 ) szukaj_zamien(LU,indeksy,i);
			wsp = LU[indeksy[j]][i]/LU[indeksy[i]][i];
			for(int k = i; k < M; k++ ){
				LU[indeksy[j]][k] = LU[indeksy[j]][k] - LU[indeksy[i]][k]*wsp;
				if( i < j ){
					LU[indeksy[j]][i] = wsp;
				}
			}
		}
	}
}

//**************************************************************************************************
//funkcja dokonująca obliczeń wektora y z równania Ly=b
void obliczanie_y(double **LU, double *Y, double *B, int *indeksy, const int n){
	double S;
	for(register  int i = 0; i < n; i++ ){
		S = 0;
		for(register  int j = 0; j < i; j++ ){
			S += Y[j]*LU[indeksy[i]][j];
		}
		Y[i] = B[indeksy[i]] - S;
	}
}

//**************************************************************************************************
//funkcja dokonująca obliczeń wektora x z równania Ux=y
void obliczanie_x(double **LU, double *X, double *Y, int *indeksy, const int n){
	double S;
	register int j;
	for( register int i = n-1; i >= 0; i-- ){
		S = 0;
		for( j = n-1; j > i; j-- )
			S += X[j]*LU[indeksy[i]][j];
		X[i] = (Y[i] - S)/LU[indeksy[i]][i];
	}
}

//**************************************************************************************************
//funkcja inicjalizująca tablicę indeksów potrzebna dla powyższych funkcji.
void init_kolumny(int *indeksy, const int n ){
	for( register int i = 0; i < n; i++ )
		indeksy[i] = i;
}

//**************************************************************************************************
void crank_nicolson_LU(double **A,double *x, double *t, double lam_pos){
	register int i,j,k;
	double *b,*z;

	double **LU = new double * [N];
	double **WSP = new double * [N];
	for( i = 0 ; i < N ; i++ ){
		LU[i] = new double [N];
		WSP[i] = new double [N];
	}
	double *Y = new double [N];
	int *indeksy = new int [N];

	b=new double[N];
	z=new double[N];

	for(j=0;j<N;j++)
		A[0][j]=0;

	for(i=0;i<M;i++){
		A[i][0]=0;
		A[i][N-1]=0;
	}

	for( i = 0; i < N; i++ )
		for(  j =0; j< N; j++ ){
			if( j == i-1 || j == i+1){
				WSP[i][j] = -lam_pos/2.0;
			}else if( j == i ){
				WSP[i][j] = 1 + lam_pos;
			}else{
				WSP[i][j] = 0;
			}
		}
	WSP[0][0] = 1;
	WSP[0][1] = 0;
	WSP[N - 1][N - 1] = 1;
	WSP[N - 1][N - 2] = 0;
		init_kolumny(indeksy,N);
		rozklad_LU(LU,WSP,indeksy);
		for( i=1;i<M;i++){
			b[0]=b[N-1]=0;
			for(j=1;j<N-1;j++)
				//b[j]=(lambda/2.0)*A[i][j-1]+(1-lambda)*A[i][j]+(lambda/2.0)*A[i][j+1];
				b[j]=((lam_pos/2.0)*A[i-1][j-1]+(1.0-lam_pos)*A[i-1][j]+(lam_pos/2.0)*A[i-1][j+1]);
			for(k=1;k<N-1;k++)
				b[k]=b[k]+(D*dt*PI*PI*sin(PI*k*h));


			obliczanie_y( LU,Y, b,indeksy,N);
			obliczanie_x( LU,z, Y,indeksy,N);

			for(j=1;j<N-1;j++)
				A[i][j]=z[j];
		}

		zapisz(A, file4);
}
double *reziduum(double **matrix, double *x, double *b) {
	double sum = 0;
	double *r = new double[N];

	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N; j++) {
			sum += matrix[i][j] * x[j];
		}
		r[i] = sum - b[i];
		sum = 0;
	}
	return r;
}
//**************************************************************************************************
//gauss Seidel
void GS(double **A, double *x, double *b){
	double tmp;
	double m=0;
	double *xp = new double[N];
	double *tmp_w = new double[N];
	int it=0;

	do{
		it++;
	for (int i = 0; i < N; i++) {
			for (int j = 0; j < N; j++) {
				if (i != j)
					m += A[i][j] * x[j];
			}
			xp[i] = x[i];
			x[i] = (1.0 / A[i][i]) * (b[i] - m);
			m = 0;
		for( int i = 0; i < N; i++ ){
			tmp_w[i] = fabs(x[i] - xp[i]);
			xp[i] = x[i];}

		}
	}while( norma(tmp_w) > 0.00000000001 );
}


void crank_nicolson_GS(double **A, double *x, double *t, double lam_pos){
	register int i,j,k;
	double *b,*z;

	double **LU = new double * [N];
	double **WSP = new double * [N];
	for( i = 0 ; i < N ; i++ ){
		LU[i] = new double [N];
		WSP[i] = new double [N];
	}

	b=new double[N];
	z=new double[N];


	for(j=0;j<N;j++){
		A[0][j]=0;
		z[j] = 1;
	}

	for(i=0;i<M;i++){
		A[i][0]=0;
		A[i][N-1]=0;
	}

	for( i = 0; i < N; i++ )
		for(  j =0; j< N; j++ ){
			if( j == i-1 || j == i+1){
				WSP[i][j] = -lam_pos/2.0;
			}else if( j == i ){
				WSP[i][j] = 1 + lam_pos;
			}else{
				WSP[i][j] = 0;
			}
		}
		WSP[0][0] = 1;
	WSP[0][1] = 0;
	WSP[N - 1][N - 1] = 1;
	WSP[N - 1][N - 2] = 0;
		for( i=1;i<M-1;i++){
			b[0]=b[N-1]=0;
			for(j=1;j<N-1;j++)
				//b[j]=(lambda/2.0)*A[i][j-1]+(1-lambda)*A[i][j]+(lambda/2.0)*A[i][j+1];
				b[j]=((lam_pos/2.0)*A[i-1][j-1]+(1.0-lam_pos)*A[i-1][j]+(lam_pos/2.0)*A[i-1][j+1]);

			for(k=1;k<N-1;k++)
				b[k]=b[k]+(D*dt*PI*PI*sin(PI*k*h));

			GS(WSP,z,b);

			for(j=1;j<N-1;j++)
				A[i][j]=z[j];
		}
		zapisz(A, file6);
}



//**************************************************************************************************
int main(){
	file.open("Analit.txt",ios::out);

	file6.open("CN-GS.txt",ios::out);
		file5.open("Bledy_CN-LU.txt",ios::out);
	file4.open("CN-LU.txt",ios::out);
	file7.open("Bledy_CN-GS.txt",ios::out);


	double *x = new double[N];
    double *t = new double[M];


    for(int i =1; i<N; i++)
            x[i] = x_min + (i*h);

    for(int i =1; i<M; i++)
            t[i] = t_min + (i*dt);

	double **A = macierz(M, N); //macierz do rozwiązania analitycznego
	double *BLAD = new double[M];
	double **CN_LU = macierz(M,N);
	double **CN_GS = macierz(M,N);
	analitycznie(A,x,t);
	crank_nicolson_GS(CN_GS,x,t,lambda_bez);
	blad(BLAD,CN_GS,A,file7);
	crank_nicolson_LU(CN_LU,x,t,lambda_bez);
	blad(BLAD,CN_LU,A,file5);

	return 0;
}