#ifndef _INCLUDED_L3_MAT5_H
#define _INCLUDED_L3_MAT5_H

// Copyright (C) Krzysztof Bosak, 1999-10-28...2000-06-03
// All rights reserved.
// kbosak@box43.pl
// http://www.kbosak.prv.pl

// *Mat5Diag with LU and Solve*

#include "l3_array.h"

template<class real>
class Mat5Diag
{	// Pentadiagonal matrix.
	basearray<real> _data;
	int _order;
	bool _decomposed;
public:
	inline explicit Mat5Diag(int order)
		: _data(order*5+10),
		_order(order),
		_decomposed(false)
	{
		assert(order>=4);
	}
	inline void Decompose()
	{
		assert(_decomposed==false);
		_decomposed=true;
		real * const E=&_data[0];
		real * const D=&_data[1];
		real * const A=&_data[2];
		real * const B=&_data[3];
		const real * const C=&_data[4];
		const int sm1=_order-1;
		int i, i5;
		for(i=1, i5=5; i<sm1; i++, i5+=5)
		{
			const int i5m5=i5-5;
			real l=D[i5]/=A[i5m5];
			A[i5]-=B[i5m5]*l;
			B[i5]-=C[i5m5]*l;
			const int i5p5=i5+5;
			// was l=E[i5p5]/=A[i5m5];!!!
			E[i5p5]/=A[i5m5];
			l=E[i5p5];
			D[i5p5]-=B[i5m5]*l;
			A[i5p5]-=C[i5m5]*l;
		}
		i5=sm1*5;
		const real l=D[i5]/=A[i5-5];
		A[i5]-=B[i5-5]*l;
	}
	inline basearray<real> Solve(
		basearray<real>& result,
		const basearray<real>& rhs
		)
	{	// Assumes that matrix is already LU decomposed!
		// result and rhs can be the same vector.
		assert(_decomposed==true);
		assert(_order==result.size());
		assert(_order==rhs.size());
		static basearray<real> y(_order);
		y.setsize(_order);
		const real * const E=&_data[0];
		const real * const D=&_data[1];
		y[0]=rhs[0];
		y[1]=rhs[1]-D[5]*rhs[0];
		int i, i5;
		for(i=2, i5=10; i<_order; i++, i5+=5)
		{
			y[i]=rhs[i]-E[i5]*y[i-2]-D[i5]*y[i-1];
		}
		const real * const A=&_data[2];
		const real * const B=&_data[3];
		const real * const C=&_data[4];

		const int sm1=_order-1;
		result[sm1]=y[sm1]/A[sm1*5];
		const int sm2=_order-2;
		result[sm2]=(y[sm2]-B[sm2*5]*result[sm1])/A[sm2*5];
		for(i=sm2-1, i5=i*5; i>=0; i--, i5-=5)
		{
			result[i]=(y[i]-B[i5]*result[i+1]-C[i5]*result[i+2])/A[i5];
		}
		return result;
	}
	inline void Swap(Mat5Diag& mat2)
	{
		assert(_order==mat2._order);
		assert(_decomposed==mat2._decomposed);
		_data.swap(mat2._data);
	}
	inline int Order() const
	{
		return _order;
	}
	inline bool Decomposed() const
	{
		return _decomposed;
	}
	inline real& U2Diag(int i)
	{
		assert(i>=0);
		assert(i<_order-2);
		return _data[i*5+4];
	}
	inline real& UDiag(int i)
	{
		assert(i>=0);
		assert(i<_order-1);
		return _data[i*5+3];
	}
	inline real& Diag(int i)
	{
		assert(i>=0);
		assert(i<_order);
		return _data[i*5+2];
	}
	inline real& LDiag(int i)
	{
		assert(i>=0);
		assert(i<_order-1);
		return _data[i*5+1+5];
	}
	inline real& L2Diag(int i)
	{
		assert(i>=0);
		assert(i<_order-2);
		return _data[i*5+10];
	}

	inline const real& GetU2Diag(int i) const
	{
		assert(i>=0);
		assert(i<_order-2);
		return _data[i*5+4];
	}
	inline const real& GetUDiag(int i) const
	{
		assert(i>=0);
		assert(i<_order-1);
		return _data[i*5+3];
	}
	inline const real& GetDiag(int i) const
	{
		assert(i>=0);
		assert(i<_order);
		return _data[i*5+2];
	}
	inline const real& GetLDiag(int i) const
	{
		assert(i>=0);
		assert(i<_order-1);
		return _data[i*5+1+5];
	}
	inline const real& GetL2Diag(int i) const
	{
		assert(i>=0);
		assert(i<_order-2);
		return _data[i*5+10];
	}
};

#endif //_INCLUDED_L3_MAT5_H