#include<iostream>
#include<vector>
#include"Tridiagonal.h"

using namespace std;



Tridiagonal::Tridiagonal()
{
	subDiagonal.push_back(0.0);
	diagonal.push_back(0.0);
	superDiagonal.push_back(0.0);
};

Tridiagonal::Tridiagonal(const Tridiagonal& B){

	subDiagonal=B.getSubDiagonal();
	diagonal=B.getDiagonal();
	superDiagonal=B.getSuperDiagonal();
};

Tridiagonal::Tridiagonal(vector<double> aVal,vector<double> bVal,vector<double> cVal)
{
	subDiagonal=aVal;
	diagonal=bVal;
	superDiagonal=cVal;
	LUDecomposition();

};

void Tridiagonal::setSubDiagonal(vector<double> &a){

	subDiagonal=a;
};

void Tridiagonal::setDiagonal(vector<double> &b){

	diagonal=b;
};

void Tridiagonal::setSuperDiagonal(vector<double> &c){

	superDiagonal=c;
};

vector<double> Tridiagonal::getSubDiagonal()const{

		return subDiagonal;

	};

vector<double> Tridiagonal::getDiagonal()const{

	return diagonal;

};

vector<double> Tridiagonal::getSuperDiagonal()const{

	return superDiagonal;

};


vector<double> Tridiagonal::getlValues()const{

	return lValues;

};

vector<double> Tridiagonal::getuValues()const{

	return uValues;
};
//Computing the LU decomposition of the tridiagonal matrix

void Tridiagonal::LUDecomposition(){

	double u=diagonal[0];
	uValues.push_back(u);

	double l=subDiagonal[0]/u;
	lValues.push_back(l);

	for(size_t k=1;k<diagonal.size();k++){

		u=diagonal[k]-lValues[k-1]*superDiagonal[k-1];
		uValues.push_back(u);
		if(k<diagonal.size()-1){
			l=subDiagonal[k]/u;
			lValues.push_back(l);
		}
	}
};
	
		


vector<double> Tridiagonal::solver(vector<double> rhs) {

		//First loop, solving Lv=rhs

		vector<double> v;

		v.push_back(rhs[0]);

		for(size_t k=1;k<rhs.size();k++){
			v.push_back(rhs[k]-lValues[k-1]*v[k-1]);
		}
		
		//Second loop, solving Uw=v

		vector<double> w;
		size_t N=v.size();
		w[N-1]=v[N-1]/uValues[N-1];

		for(size_t k=N-2;k>=0;k--){
			w[k]=(v[k]-superDiagonal[k]*w[k+1])/uValues[k];
		}

		return w;

	}