#ifndef CTRIDIAG_H
#define CTRIDIAG_H
#include<complex>
#include"Cvec.h"
#include"Wavefunction.h"
typedef std::complex<double> Creal;

class Ctridiag {
public:
	Cvec sup;
	Cvec sub;
	Cvec diag;

	Ctridiag();
	// initialize Tridiagmatrix with n diagonal elements and (n-1) upper/lower diagonal elements
	Ctridiag(const int n);
	// initialize Tridiagmatrix with entries of "in"
	Ctridiag(const Ctridiag& in);
	// initialize Tridiagmatrix with n real diagonal elements (diag_in) and (n-1) upper(sup_in)/lower(sub_in) elements
	Ctridiag(const int n, double sub_in, double diag_in, double sup_in);
	// initialize Tridiagmatrix with n complex diagonal elements (diag_in) and (n-1) upper(sup_in)/lower(sub_in) elements
	Ctridiag(const int n, Creal& sub_in, Creal& diag_in, Creal& sup_in);
	~Ctridiag();

	// matrix-vector-multiplication (result: vector)
	Cvec operator*(const Cvec& invec);
	// matrix-vector-multiplication (result: wavefunctionvector)
	Wavefunction operator*(const Wavefunction& invec);
	// efficient matrix-vector-multiplication (result: wavefunctionvector)
	void matmul(Wavefunction& res, const Wavefunction& invec);
	// multiply every element of Ctridiag with (real) "in"
	Ctridiag operator*(const double& in);
	// multiply every element of Ctridiag with (complex) "in"
	Ctridiag operator*(const Creal& in);

	// set sub/sup to 0 and diag to diag_in 
	Ctridiag operator=(double& diag_in);
	// set sub/sup to 0-vectors with length (in.size-1) and diag to "in"
	Ctridiag operator=(Cvec& in); 

	// element-by-element addition of Ctridiag and "in"
	Ctridiag operator+(const Ctridiag& in);
	// element-by-element subtraction of "in" from Ctridiag
	Ctridiag operator-(const Ctridiag& in);

	// divide every element of Ctridiag by "in"
	Ctridiag operator/(double& in);

	// solve equation M*x = b with Tridiagmatrix M and known vector b
	Cvec solve(const Cvec& b) const; 
	// additional friend operator/  to solve x = b/M
	friend Cvec operator/(const Cvec& bin, const Ctridiag& Min) {
		return Min.solve(bin);
	}

	// solve equation M*x = b with Tridiagmatrix M and known wavefunction b
	Wavefunction solve(const Wavefunction& b) const;
	// additional friend operator/  to solve x = b/M
	friend Wavefunction operator/(const Wavefunction& bin, const Ctridiag& Min) {
		return Min.solve(bin);
	}
	// efficient solve operator
	void solve_ref(Wavefunction& x, const Wavefunction& b, Ctridiag& a);

	void init(int n);
	// print matrix (first column: sub, second: diag, third: sup)
	void print();
};

#endif
