#ifndef _INCLUDED_L3_MAT_H
#define _INCLUDED_L3_MAT_H

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

// *LU_Decomposition - multiple solution of system with Decompose, Solve, SolvePack*

#include <string.h>
#include <math.h>
#include <stdlib.h>
#include "l3_array.h"
#include "l3_str.h"
#include "l3_tool.h"

/////////////////////////////////////////////////////////////////////////////

bool SolveMat2x2(double coeff[], double * x, double * y) // coeff[6]
{
	// Solves 2x2 linear system using single (row) pivoting Gauss LU decomposition:
	// Ax+By=E
	// Cx+Dy=F
	assert(coeff!=NULL);
	assert(x!=NULL);
	assert(y!=NULL);
	double & A=coeff[0];
	double & B=coeff[1];
	double & C=coeff[2];
	double & D=coeff[3];
	if(A*D==B*C)
	{
		return false;
	}
	double& E=coeff[4];
	double& F=coeff[5];
	if(A==0)
	{
		const double temp1=A;
		const double temp2=B;
		const double temp3=E;
		A=C;
		B=D;
		E=F;
		C=temp1;
		D=temp2;
		F=temp3;
	}
	if(C!=0)
	{
		C/=A;
		D-=C*B;
		// LU decomposed.
		F-=C*E;
	}
	*y=F/D;
	*x=(E-B*(*y))/A;

	return true;
}

/////////////////////////////////////////////////////////////////////////////

/*
// Disabled because of bug in dmc Digital Mars
bool SolveMat3x3(double coeff[], double * x, double * y, double * z)// coeff[12]
{
	// Solves 3x3 linear system using single (row) pivoting Gauss LU decomposition:
	// Ax+By+Cz=J
	// Dx+Ey+Fz=K
	// Gx+Hy+Iz=L
	assert(coeff!=NULL);
	assert(x!=NULL);
	assert(y!=NULL);
	assert(z!=NULL);
	double & A=coeff[0];
	double & B=coeff[1];
	double & C=coeff[2];
	double & D=coeff[3];
	double & E=coeff[4];
	double & F=coeff[5];
	double & G=coeff[6];
	double & H=coeff[7];
	double & I=coeff[8];

	if(A*E*I+B*F*G+C*D*H==A*F*H+C*E*G+B*D*I)
	{
		return false;
	}
	double& J=coeff[9];
	double& K=coeff[10];
	double& L=coeff[11];

	if(A==0)
	{
		const double temp1=A;
		const double temp2=B;
		const double temp3=C;
		const double temp4=J;
		if(D!=0)
		{
			A=D;
			B=E;
			C=F;
			J=K;
			D=temp1;
			E=temp2;
			F=temp3;
			K=temp4;
		}
		else // G!=0
		{
			assert(G!=0);
			A=G;
			B=H;
			C=I;
			J=L;
			G=temp1;
			H=temp2;
			I=temp3;
			L=temp4;
		}
	}
	if(D!=0)
	{
		D/=A;
		E-=D*B;
		F-=D*C;
	}
	if(G!=0)
	{
		G/=A;
		H-=G*B;
		I-=G*C;
	}
	if(E==0)
	{
		const double temp1=D;
		const double temp2=E;
		const double temp3=F;
		const double temp4=K;
		D=G;
		E=H;
		F=I;
		K=L;
		G=temp1;
		H=temp2;
		I=temp3;
		L=temp4;
	}
	if(H!=0)
	{
		H/=E;
		I-=H*F;
	}

	// LU decomposed.
	K-=D*J;
	L-=(G*J+H*K);
	assert(I!=0);
	*z=L/I;
	assert(E!=0);
	*y=(K-F*(*z))/E;
	assert(A!=0);
	*x=(J-B*(*y)-C*(*z))/A;
	return true;
}
*/

/////////////////////////////////////////////////////////////////////////////

template<class real>
class Single_LU_Decomposition
{	// After DecomposeAndSolve() matrix content is changed,
	// and contains only U part and diagonal.
	// DecomposeAndSolve() is available only once.
	int _order;
	basearray<real> _data;
	basearray<int> _rowindex;
	bool _decomposed;
	basearray<int> _resultorder;
	basearray<real> _resultsorted;

public:
	explicit Single_LU_Decomposition(int order);
	inline int Order() const
	{
		return _order;
	}
	inline void Erase()
	{	// Optional in use, sets all matrix elements to 0.
		// If not initialized, content of matrix is garbage!
		_data.erase();
	}
	inline void Reset()
	{	// You should fill whole matrix and rhs vector with
		// new values before next DecomposeAndSolve() call.
		// If not initialized, content of matrix is garbage!
		const int order_p1=_order+1;
		for(int i=0, k=0; i<_order; i++, k+=order_p1)
		{
			_rowindex[i]=k;
		}
		_decomposed=false;
	}
	inline real& operator()(int x, int y)
	{
		assert(x>=0);
		assert(x<_order);
		assert(y>=0);
		assert(y<_order);
		return _data[x+_rowindex[y]];
	}
	inline real operator()(int x, int y) const
	{
		assert(x>=0);
		assert(x<_order);
		assert(y>=0);
		assert(y<_order);
		return _data[x+_rowindex[y]];
	}
	inline real& RHS(int y)
	{
		assert(y>=0);
		assert(y<_order);
		return _data[_order+_rowindex[y]];
	}
	inline real RHS(int y) const
	{
		assert(y>=0);
		assert(y<_order);
		return _data[_order+_rowindex[y]];
	}
	inline void SetRHS(const basearray<real>& rhs)
	{
		assert(rhs.size()==_order);
		for(int y=0; y<_order; y++)
		{
			_data[_order+_rowindex[y]]=rhs[y];
		}
	}
	void ApplyDiagonalPreconditioner();
	int DecomposeAndSolve(basearray<real>& result);
	void Show() const;
};

template<class real>
Single_LU_Decomposition<real>::Single_LU_Decomposition(int order)
	: _order(order),
	_data(order*order+order),
	_rowindex(order),
	_decomposed(false),
	_resultorder(order),
	_resultsorted(order)
{
	assert(order>=1);
	Reset();
}
template<class real>
void Single_LU_Decomposition<real>::ApplyDiagonalPreconditioner()
{	// Forget it with full pivoting! This is better to not rescale at all,
	// because full pivoting will always choose the best pivot.
	// Applying diagonal Jacobi preconditioner.
	// Use after assembling the matrix, before DecomposeAndSolve().
	// Also known as explicit scaling.
	for(int y=0; y<_order; y++)
	{
		real sumabs=0;
		for(int i=0; i<_order; i++)
		{
			sumabs+=fabs(operator()(i, y));
		}
		//const real preconditioner_diagonal_value=1./exp(- floor(1+(log(sumabs)/log(2))) *log(2));
			// can get out of range for exp or log
		const real preconditioner_diagonal_value=operator()(y, y);
		//const real preconditioner_diagonal_value=1;
		if(preconditioner_diagonal_value!=0)
	 	{
			for(int x=0; x<_order; x++)
			{
				operator()(x, y)/=preconditioner_diagonal_value;
			}
			RHS(y)/=preconditioner_diagonal_value;
		}
	}
}
template<class real>
int Single_LU_Decomposition<real>::DecomposeAndSolve(basearray<real>& result)
{
	// Some tests for 100% dense matrices:
	// CeleronA416Mhz, 83Mhz bus: (LU+Solve)
	//Microsoft VC++5: 0.01s for 100, 2.2s for 400
	//GNU gcc: 1.3s for 400
	//KAI KCC: 1.3s for 400
	assert(_decomposed==false);
	assert(result.size()==_order);
	_decomposed=true;
	int i;
	for(i=0; i<_resultorder.size(); i++)
	{
		_resultorder[i]=i;
	}
	for(i=0; i<_order; i++)
	{
		// Full pivoting:
		int bestpx=i;
		int bestpy=i;

		real bestpivot=0;
		for(int py=i; py<_order; py++)
		{
			for(int px=i; px<_order; px++)
			{
				const real testval=fabs(_data[_rowindex[py]+px]);
				if(testval>bestpivot)
				{
					bestpx=px;
					bestpy=py;
					bestpivot=testval;
				}
			}
		}
		// Fast column pivoting:
		if(bestpy!=i)
		{
			const int temp=_rowindex[i];
			_rowindex[i]=_rowindex[bestpy];
			_rowindex[bestpy]=temp;
		}
		// Slow row pivoting (but Please don't erase this part,
		// because 'fast' pivoting implementation depends on this):
		if(bestpx!=i)
		{
			const int temp=_resultorder[i];
			_resultorder[i]=_resultorder[bestpx];
			_resultorder[bestpx]=temp;
			for(int t=0; t<_order; t++)
			{
				const real tempval=_data[_rowindex[t]+i];
				_data[_rowindex[t]+i]=_data[_rowindex[t]+bestpx];
				_data[_rowindex[t]+bestpx]=tempval;
			}
		}
		//

		const int ii=_rowindex[i];
		const real pivot=_data[i+ii];
		//assert(pivot!=0);
		if(pivot==0)// Defensive.
		{
			result.erase();
			return 1;
		}
		for(int j=i+1; j<_order; j++)
		{
			const int jj=_rowindex[j];
			real elimin=_data[i+jj];

			const real temp=elimin/pivot;
			real *dat=&_data[i+1];
			for(int k=i+1; k<=_order; k++, dat++)
				dat[jj]-=dat[ii]*temp;
		}
	}
	for(i=_order-1; i>=0; i--)
	{
		real tmp=0;
		real * const dat=&_data[_rowindex[i]];
		for(int k=i+1; k<_order; k++)
		{
			tmp+=dat[k]*result[k];
		}
		result[i]=(dat[_order]-tmp)/dat[i];
	}
	for(i=0; i<_order; i++)
	{
		_resultsorted[_resultorder[i]]=result[i];
	}
	_resultsorted.swap(result);
	return 0;
}
template<class real>
void Single_LU_Decomposition<real>::Show() const
{
	for(int y=0; y<_order; y++)
	{
		cout.precision(4);
		for(int x=0; x<_order+1; x++)
		{
			cout<<_data[x+_rowindex[y]]<<'\t';
		}
		//cout<<"\tInd="<<_rowindex[y];
		cout<<endl;
	}
	for(int i=0; i<_order; i++)
	{
		cout<<_rowindex[i]<<',';
	}
	cout<<endl<<endl;
}

/////////////////////////////////////////////////////////////////////////////

template<class real>
class LU_Decomposition
{	// This LU is optimized for narrow-band, mostly sparse matrix with a lot of null elements
	// but does fine with ANY reversible matrix, being not worse than 20% than classic
	// LU for completely filled matrix. Memory is allocated for the whole matrix,
	// there is no cure and often no problem with this. Outperforms practically all well-known
	// sparse matrices (based on inefficiently managed lists) thanks to inherent simplicity
	// of code. Sparse version (for really large matrices, larger than RAM, but slower for
	// the smaller ones) based on array-based list (but who uses such a fast list?) coming soon...
	// f.ex. for 1.62%-filled 402x402 elements narrow-band sparse matrix this LU is 36 times faster
	// than classic LU, now takes 0.036s for LU (CeleronA 416MHz, 83bus).
	// Solve and SolvePack are optimized for very sparse vectors, but who knows the timing
	// for almost completely filled vectors. Theoretically I don't see there any penalty.
	// For finite element applications one could try to write GaussLU using the fact that
	// positions of nonnull elements are symmetric in such matrices (but thos matrices are not
	// symmetric in mathematical sense, transpose(A)!=A). This will complicate a little some of
	// the loops, rendering whole decomposition slower! You bet.
	int _order;
	int _orderpow2;
	basearray<real> _data;// 2D array of matrix elements, accessed linearly
	basearray<int> _rowindex;// linear indexes of following rows, can change during LU
	basearray<int> _rhs_rowindex;// corresponds to _rowindex but contains vector composants, not rows of matrix
	basearray<int> _firstnz;// first and nonnull element for every column
	basearray<int> _lastnz;// first and nonnull element for every column
	bool _decomposed;// true: _data contains decomposed L and U matrices,
					//false: contains original matrix; changed by LU, ReverseLU
	int _swap_rows;// numbers of row swaps, updated during LU
	mutable real _determinant;// updated by Determinant, after LU; otherwise is null
public:
	LU_Decomposition()
	{
	}
	explicit LU_Decomposition(int order);
	#ifdef _INCLUDED_L3_MAT3_H
		LU_Decomposition(const Mat3Diag<real>& m3);
		LU_Decomposition(const Mat3DiagSym<real>& m3);
		LU_Decomposition(const Mat3DiagCycl<real>& m3);
	#endif // _INCLUDED_L3_MAT3_H
	#ifdef _INCLUDED_L3_MAT5_H
		LU_Decomposition(const Mat5Diag<real>& m5);
	#endif // _INCLUDED_L3_MAT5_H
	inline void Erase()
	{	// Optional in use, sets all matrix elements to 0.
		// If not initialized, content of matrix is garbage!
		for(int i=0; i<_data.size(); i++)
			_data[i]=0;
	}
	inline bool Decomposed() const
	{
		return _decomposed;
	}
	inline real& operator()(int x, int y)
	{
		// Direct, fast access to matrix element, makes good abstract sense before LU
		// After LU should use L() or U() or even operator()(), but carefully.
		// x==row, y==column, using C++ array indexing
		// Changes of element ARE ALLOWED using this function definition
		assert(x>=0);
		assert(x<_order);
		assert(y>=0);
		assert(y<_order);
		return _data[x+_rowindex[y]];
	}
	inline real operator()(int x, int y) const
	{	// Direct, fast access to matrix element, makes good abstract sense before LU
		// After LU one should use L() or U() or even faster: use operator(), but carefully:
		// L has 1 on diagonal, so diagonal elements are those of U matrix
		// x==row, y==column, using C++ array indexing
		// Changes of element ARE NOT POSSIBLE using this function definition
		assert(x>=0);
		assert(x<_order);
		assert(y>=0);
		assert(y<_order);
		return _data[x+_rowindex[y]];
	}
	inline real L(int x, int y) const
	{	// It makes sense to call it only after LU to get access to elements
		// of lower-diagonal decomposed matrix
		assert(_decomposed==true);
		assert(x>=0);
		assert(x<_order);
		assert(y>=0);
		assert(y<_order);
		if(x<y)
			return _data[x+_rowindex[y]];
		if(x>y)
			return 0;
		return 1;
	}
	inline real U(int x, int y) const
	{	// It makes sense to call it only after LU to get access to elements
		// of upper-diagonal decomposed matrix
		assert(_decomposed==true);
		assert(x>=0);
		assert(x<_order);
		assert(y>=0);
		assert(y<_order);
		if(x>=y)
			return _data[x+_rowindex[y]];
		return 0;
	}
	inline basearray<real> GetMatrix() const
	{	// Raw, fast access to matrix elements.
		assert(_decomposed==false);
		return _data;
	}
	inline void Show() const
	{	// Prints matrix to the screen, and rowindex too.
		for(int y=0; y<_order; y++)
		{
			cout.precision(4);
			for(int x=0; x<_order; x++)
				cout<<_data[x+_rowindex[y]]<<'\t';
			cout<<endl;
		}
		for(int i=0; i<_order; i++)
			cout<<_rowindex[i]<<',';
		cout<<endl<<endl;
	}
	inline int Order() const
	{	// Returns order (and width and height) of matrix.
		return _order;
	}
	inline int MemoryUsage() const
	{	// Returns reliable, a few bytes underestimate of overall memory usage of this object, in bytes
		return _orderpow2*sizeof(real)+sizeof(int)*4*_order+sizeof(*this);
	}

	void ApplyAndGetDiagonalPreconditioner(basearray<real>& preconditioner);
#ifdef PLATFORM_HAS_FILESYSTEM
	void LoadTXT(const TextString& matrix_filename);
	void SaveTXT(const TextString& matrix_filename) const;
#endif // PLATFORM_HAS_FILESYSTEM
	int Decompose();
	real Determinant() const;
	real Condition() const;
	void Solve(basearray<real>& result, const basearray<real>& rhs) const;
	void Solve_IterativeImprovement(
		const LU_Decomposition<real>& original,
		const basearray<real>& rhs,
		basearray<real>& solution
		) const;
	void SolvePack(
		autoarray<basearray<real> >& resultpack,
		const autoarray<basearray<real> >& rhspack
		) const;
	bool IsSymmetric() const;
	void ReverseDecompose(LU_Decomposition& target) const;
	void Compare(const LU_Decomposition& error_reference) const;
	real TestAnySolution(const basearray<real>& result,
		const basearray<real>& rhs,
		real eps=0
		) const;
	void TestSolve(const basearray<real>& rhs, real eps=0) const;
	void TestSolvePack(
		const autoarray<basearray<real> >& rhspack,
		real eps=0
		) const;
	void TestDecompose() const;
	void TestMatrix() const;
	void SaveGIF(const TextString& gif_filename) const;
};

template<class real>
LU_Decomposition<real>::LU_Decomposition(int order)
	:_order(order),
	_orderpow2(order*order),
	_data(order*order),
	_rowindex(order),
	_rhs_rowindex(order),
	_firstnz(order),
	_lastnz(order),
	_decomposed(false),
	_swap_rows(0),
	_determinant(0)
{	// Matrix of linear system, once constructed, will not be resized.
	// only initial order is required, then manual initiation with operator(), then Decompose(),
	// then Solve() or SolvePack()
	assert(order>=2);
	int i;
	int k;
	for(i=0, k=0; i<_order; i++, k+=_order)
	{
		_rowindex[i]=k;
		_rhs_rowindex[i]=i;
		_firstnz[i]=0;
		_lastnz[i]=_order-1;
	}
}
#ifdef _INCLUDED_L3_MAT3_H
template<class real>
LU_Decomposition<real>::LU_Decomposition(const Mat3Diag<real>& m3)
	:_order(m3.Order()), _orderpow2(_order*_order), _data(_order*_order),
	_rowindex(_order), _rhs_rowindex(_order), _firstnz(_order), _lastnz(_order),
	_decomposed(m3.Decomposed()), _swap_rows(0), _determinant(0)
{	// You can copy values from Mat3Diag class and use test functions on it.
	assert(_order>=21);
	Erase();
	const int om1=_order-1;
	int i;
	int k;
	for(i=0, k=0; i<_order; i++, k+=_order)
	{
		_rowindex[i]=k;
		_rhs_rowindex[i]=i;
		operator()(i, i)=m3.GetDiagonal()[i];
		_firstnz[i]=0;
		_lastnz[i]=om1;
	}
	for(i=0; i<om1; i++)
	{
		operator()(i+1, i)=m3.GetUpperDiagonal()[i];
		operator()(i, i+1)=m3.GetLowerDiagonal()[i];
	}
}
template<class real>
LU_Decomposition<real>::LU_Decomposition(const Mat3DiagSym<real>& m3)
	:_order(m3.Order()), _orderpow2(_order*_order), _data(_order*_order),
	_rowindex(_order), _rhs_rowindex(_order), _firstnz(_order), _lastnz(_order),
	_decomposed(m3.Decomposed()), _swap_rows(0), _determinant(0)
{	// You can copy values from Mat3DiagSym class and use test functions on it.
	assert(_order>=2);
	Erase();
	const int om1=_order-1;
	int i;
	int k;
	for(i=0, k=0; i<_order; i++, k+=_order)
	{
		_rowindex[i]=k;
		_rhs_rowindex[i]=i;
		operator()(i, i)=m3.GetDiagonal()[i];
		_firstnz[i]=0;
		_lastnz[i]=om1;
	}
	for(i=0; i<om1; i++)
	{
		const real temp=GetUpperDiagonal()[i];
		operator()(i+1, i)=temp;
		operator()(i, i+1)=temp;
	}
}
template<class real>
LU_Decomposition<real>::LU_Decomposition(const Mat3DiagCycl<real>& m3)
	:_order(m3.Order()), _orderpow2(_order*_order), _data(_order*_order),
	_rowindex(_order), _rhs_rowindex(_order), _firstnz(_order), _lastnz(_order),
	_decomposed(m3.Decomposed()), _swap_rows(0), _determinant(0)
{	// You can copy values from Mat3DiagCycl class and use test functions on it.
	assert(_order>=2);
	Erase();
	const int om1=_order-1;
	int i;
	int k;
	for(i=0, k=0; i<_order; i++, k+=_order)
	{
		_rowindex[i]=k;
		_rhs_rowindex[i]=i;
		operator()(i, i)=m3.GetDiagonal()[i];
		_firstnz[i]=0;
		_lastnz[i]=om1;
	}
	for(i=0; i<om1; i++)
	{
		operator()(i+1, i)=m3.GetUpperDiagonal()[i];
		operator()(i, i+1)=m3.GetLowerDiagonal()[i];
	}
	const int om2=_order-2;
	for(i=0; i<om2; i++)
	{
		operator()(_order-1, i)=m3.GetRightColumn()[i];
		operator()(i, _order-1)=m3.GetLowerRow()[i];
	}
}
#endif // _INCLUDED_L3_MAT3_H
#ifdef _INCLUDED_L3_MAT5_H
template<class real>
LU_Decomposition<real>::LU_Decomposition(const Mat5Diag<real>& m5)
	:_order(m5.Order()), _orderpow2(_order*_order), _data(_order*_order),
	_rowindex(_order), _rhs_rowindex(_order), _firstnz(_order), _lastnz(_order),
	_decomposed(m5.Decomposed()), _swap_rows(0), _determinant(0)
{	// You can copy values from Mat5Diag class and use test functions on it
	assert(_order>=2);
	Erase();
	const int om1=_order-1;
	int i;
	int k;
	for(i=0, k=0; i<_order; i++, k+=_order)
	{
		_rowindex[i]=k;
		_rhs_rowindex[i]=i;
		operator()(i, i)=m5.GetDiag(i);
		_firstnz[i]=0;
		_lastnz[i]=om1;
	}
	for(i=0; i<om1; i++)
	{
		operator()(i+1, i)=m5.GetUDiag(i);
		operator()(i, i+1)=m5.GetLDiag(i);
	}
	const int om2=_order-2;
	for(i=0; i<om2; i++)
	{
		operator()(i+2, i)=m5.GetU2Diag(i);
		operator()(i, i+2)=m5.GetL2Diag(i);
	}
}
#endif // _INCLUDED_L3_MAT5_H

template<class real>
void LU_Decomposition<real>::ApplyAndGetDiagonalPreconditioner(basearray<real>& preconditioner)
{	// Applying diagonal Jacobi preconditioner.
	// Use after assembling the matrix, before DecomposeAndSolve().
	// Also known as explicit scaling.
	// Preconditioner is stored in preconditioner array.
	// Don't forget to divise rhs by preconditioner before
	// each Solve() call!
	assert(preconditioner.size()==_order);
	for(int y=0; y<_order; y++)
	{
		const real preconditioner_diagonal_value=operator()(y, y);
		if(preconditioner_diagonal_value!=0)
		{
			preconditioner[y]=preconditioner_diagonal_value;
			for(int x=0; x<_order; x++)
			{
				if(x==y)
				{
					operator()(x, y)=1;
					continue;
				}
				operator()(x, y)/=preconditioner_diagonal_value;
			}
			// RHS(y)/=preconditioner_diagonal_value;
			// Don't forget to divise rhs by preconditioner before each Solve() call!
		}
		else
		{
			preconditioner[y]=1;
		}
	}
}

#ifdef PLATFORM_HAS_FILESYSTEM

template<class real>
void LU_Decomposition<real>::LoadTXT(const TextString& matrix_filename)
{	// Loads whole matrix from a file (not its LU decomposed form!).
	// Width and Height numbers which are preceeding the data are tested for conformance.
	ifstream datfile(matrix_filename.c_str(), ios::in|ios::nocreate);
	PFileGood(datfile, matrix_filename.c_str());
	int temp;
	datfile>>temp;
	assert(!datfile.fail());
	*this=LU_Decomposition(temp);
	for(int y=0; y<_order; y++)
	{
		for(int x=0; x<_order; x++)
		{
			PFileGood(datfile, matrix_filename.c_str());
			datfile>>_data[x+_rowindex[y]];
		}
	}

	/*
	if(!datfile.fail())
	{// Still some unreliable, false warnings (remaining whitespaces etc.).
		cerr<<"WARNING in LU_Decomposition::LoadTXT()\n";
		cerr<<"There are still more data after matrix elements in file "<<matrix_filename<<endl;
	}*/
	datfile.close();
}

template<class real>
void LU_Decomposition<real>::SaveTXT(const TextString& matrix_filename) const
{	// Saves whole matrix (or its LU decomposed form)
	// to a file as plain text of tabularized numbers, row-major order
	// Width and Height numbers are preceeding the data, each row ends by end-of-line.
	ofstream datfile(matrix_filename.c_str());
	PFileGood(datfile, matrix_filename.c_str());
	datfile<<_order<<endl;
	datfile.precision(15);
	for(int y=0; y<_order; y++)
	{
		for(int x=0; x<_order; x++)
		{
			PFileGood(datfile, matrix_filename.c_str());
			datfile<<_data[x+_rowindex[y]]<<'\t';
		}
		PFileGood(datfile, matrix_filename.c_str());
		datfile<<endl;
	}
	PFileGood(datfile, matrix_filename.c_str());
	datfile.close();
}

#endif // PLATFORM_HAS_FILESYSTEM

template<class real>
int LU_Decomposition<real>::Decompose()
{	// LU decomposition of matrix of linear system using single pivoting,
	// doesn't tests for reversibility of matrix (assumes to be reversible).
	// To be called only once after having filled desired elements of the matrix,
	// please CLEAR all null elements manually before LU! see operator()().
	// Use calls of Solve, SolvePack to obtain resolutions of linear system.
	// vs. double: long double +140% slower, with float +20% faster
	assert(_decomposed==false);
	assert(_swap_rows==0);
	_decomposed=true;
	_determinant=0;
	int i;
	for(i=0; i<_order; i++)
	{
		int j;

		int imax=i;

		//7% slower for sparse (402x402, 1.62% fill) matrix due to single pivoting...
		// Single pivoting:
		real maxval=_data[_rowindex[imax]+i];
		if(maxval<0)
			maxval=-maxval;
		for(j=i+1; j<_order; j++)
		{
			real testval=_data[_rowindex[j]+i];
			if(testval!=0)
			{
				if(testval<0)
					testval=-testval;
				if(testval>maxval)
				{
					imax=j;
					maxval=testval;
				}
			}
		}

		/*
		#ifndef NDEBUG
		real totalmax=0;
		for(int jj=i; jj<_order; jj++)
		{
			for(int ii=i; ii<_order; ii++)
			{
				const real testmax=fabs(_data[_rowindex[jj]+ii]);
				if(testmax>totalmax)
				{
					totalmax=testmax;
				}
			}
		}
		if(totalmax>maxval)
		{
			cout<<'*'<<totalmax/maxval<<endl;
		}
		#endif
		*/

		if(imax!=i)
		{
			int temp=_rowindex[i];
			_rowindex[i]=_rowindex[imax];
			_rowindex[imax]=temp;
			temp=_rhs_rowindex[i];
			_rhs_rowindex[i]=_rhs_rowindex[imax];
			_rhs_rowindex[imax]=temp;
			_swap_rows++;
		}

		real pivot=_data[_rowindex[i]+i];
		assert(pivot!=0);// This condition fails always for singular matrices
		if(pivot==0)// Defensive.
		{
			return 1;
		}

		const int ii=_rowindex[i];
		for(j=i+1; j<_order; j++)
		{
			const int jj=_rowindex[j];
			real& pivot=_data[i+jj];
			if(pivot!=0)
			{
				pivot/=_data[i+ii];
				const real temp=pivot;
				if(temp==1)//1-3% slower due to 2 additional conditions: 1 and -1
				{
					real * dat=&_data[i+1];
					for(int k=i+1; k<_order; k++, dat++)
						dat[jj]-=dat[ii];
				}
				else if(temp==-1)
				{
					real *dat=&_data[i+1];
					for(int k=i+1; k<_order; k++, dat++)
						dat[jj]+=dat[ii];
				}
				else
				{
					real * dat=&_data[i+1];
					for(int k=i+1; k<_order; k++, dat++)
							dat[jj]-=dat[ii]*temp;
				}
			}
		}
	}
	const real * dat=&_data[0];
	for(int x=0; x<_order; x++, dat++)
	{
		int y;

		for(y=0; y<_order; y++)
		{
			if(dat[_rowindex[y]]!=0)
			{
				_firstnz[x]=y;
				break;
			}
		}
		for(y=_order-1; y>=0; y--)
		{
			if(dat[_rowindex[y]]!=0)
			{
				_lastnz[x]=y;
				break;
			}
		}
	}//17% slower due to non-zero element determination
	return 0;
}

template<class real>
real LU_Decomposition<real>:: Determinant() const
{	// Determinant of matrix can be calculated after LU decomposition.
	// Reversible matrix (only which can pass LU) has always det!=0.
	assert(_decomposed==true);
	if(_determinant!=0)
	{
		return _determinant;
	}
	_determinant=operator()(0, 0);
	for(int i=1; i<_order; i++)
	{
		_determinant*=operator()(i, i);
	}
	assert(_determinant!=0);
	if( (_swap_rows&1) !=0 )// even number of row swaps
	{
		_determinant=-_determinant;
	}
	return _determinant;
}

template<class real>
real LU_Decomposition<real>:: Condition() const
{	// Call it before Decomposition(), since needs original matrix!
	// Condition is an estimate of the condition of matrix A.
	// For the linear system A*x=b, changes in A and b
	// may cause changes cond times as large in x.
	// If Condition+1=Condition, A is singular to working precision.
	// Condition=(1-norm of A)*(an estimate of 1-norm of A-inverse)
	// estimate obtained by one step of inverse iteration for the
	// small singular vector. This involves solving two systems
	// of equations, (A_transpose)*y=e and A*z=y where e
	// is a vector of +1 chosen to cause growth in y.
	// Condition_estimate=(1-norm of z)/(1-norm of y)
	assert(_decomposed==false);
	LU_Decomposition copy(*this);
	real norm_a=0;
	for(int j=0; j<_order; j++)
	{
		real t=0;
		for(int i=0; i<_order; i++)
		{
			real val=operator()(i, j);
			if(val<0)
				val=-val;
			t+=val;
		}
		if(t>norm_a)
			norm_a=t;
	}
	copy.Decompose();
	basearray<real> e(_order);
	e.fill(1);
	basearray<real> y(_order);
	copy.Solve(y, e);
	basearray<real> z(_order);
	copy.Solve(z, y);
	real norm_y=0;
	real norm_z=0;
	for(int i=0; i<_order; i++)
	{
		real valy=y[i];
		if(valy<0)
			valy=-valy;
		norm_y+=valy;
		real valz=z[i];
		if(valz<0)
			valz=-valz;
		norm_z+=valz;
	}
	assert(norm_y!=0);
	real condition=norm_a*norm_z/norm_y;
	if(condition<1)
		condition=1;
	return condition;
}

template<class real>
void LU_Decomposition<real>:: Solve(basearray<real>& result, const basearray<real>& rhs) const
{	// Solves Ax=b, rhs is right-hand side (b) on input, result contains x on output.
	// A, it is matrix this linear system, must be already decomposed (LU).
	// Use Solvepack to solve multiple vectors at once, much faster.
	// Speed vs double: +150% slower using long double, +50% faster using float.
	// result and rhs can be the same vector.
	assert(_decomposed==true);
	assert(rhs.size()==_order);
	assert(result.size()==_order);
	//assert(&result != &rhs);
	/*
	int i;
	for(i=0; i<_order; i++)
		result[i]=rhs[_rhs_rowindex[i]];
	for(i=0; i<_order; i++)
	{
		const real temp=result[i];
		if(temp!=0)
		{
			const real * const dat=&_data[i];
			const int jend=_lastnz[i];
			for(int j=i+1; j<=jend; j++)
				result[j]-=dat[_rowindex[j]]*temp;
		}
	}
	for(i=_order-1; i>=0; i--)
	{
		const real * const dat=&_data[i];
		if(result[i]!=0)
		{
			result[i]/=dat[_rowindex[i]];
			const real temp=result[i];
			const int jend=_firstnz[i];
			for(int j=i-1; j>=jend; j--)
				result[j]-=dat[_rowindex[j]]*temp;
		}
	}
	*/
	int i;

	const int orderm1=_order-1;
	i=0;
	while(i<orderm1)
	{
		result[i]=rhs[_rhs_rowindex[i]];
		i++;
		result[i]=rhs[_rhs_rowindex[i]];
		i++;
	}
	if(i<_order)
	{
		result[i]=rhs[_rhs_rowindex[i]];
	}
	for(i=0; i<_order; i++)
	{
		const real temp=result[i];
		//if(temp!=0)
		{
			const real * const dat=&_data[i];
			const int jend=_lastnz[i];
			/*
			for(int j=i+1; j<=jend; j++)
				result[j]-=dat[_rowindex[j]]*temp;
			*/
			int j=i+1;
			const int jendm1=jend-1;
			while(j<=jendm1)
			{
				result[j]-=dat[_rowindex[j]]*temp;
				j++;
				result[j]-=dat[_rowindex[j]]*temp;
				j++;
			}
			if(j<=jend)
			{
				result[j]-=dat[_rowindex[j]]*temp;
			}
		}
	}
	for(i=_order-1; i>=0; i--)
	{
		const real * const dat=&_data[i];
		//if(result[i]!=0)
		{
			result[i]/=dat[_rowindex[i]];
			const real temp=result[i];
			const int jend=_firstnz[i];
			/*
			for(int j=i-1; j>=jend; j--)
				result[j]-=dat[_rowindex[j]]*temp;
			*/
			int j=i-1;
			const int jendp1=jend+1;
			while(j>=jendp1)
			{
				result[j]-=dat[_rowindex[j]]*temp;
				j--;
				result[j]-=dat[_rowindex[j]]*temp;
				j--;
			}
			if(j>=jend)
			{
				result[j]-=dat[_rowindex[j]]*temp;
			}
		}
	}
}

template<class real>
void LU_Decomposition<real>:: Solve_IterativeImprovement(
	const LU_Decomposition<real>& original,
	const basearray<real>& rhs,
	basearray<real>& solution
	) const
{	// It is recommended to call it at least once.
	// To do it you need to make a copy of matrix before decomposition
	// and to use it together with right-hand side
	// and almost exact solution provided by 'Solve' function member.
	// The goal is at least one order better accuracy.
	assert(_decomposed==true);
	assert(original._decomposed==false);
	assert(rhs.size()==_order);
	assert(original.Order()==_order);
	assert(solution.size()==_order);
	basearray<real> local_rhs(_order);
	for(int y=0; y<_order; y++)
	{
		local_rhs[y]=0;
		for(int x=0; x<_order; x++)
		{
			local_rhs[y]+=original(x, y)*solution[x];
		}
		local_rhs[y]-=rhs[y];
	}
	basearray<real> local_dx(_order);
	Solve(local_dx, local_rhs);
	for(int i=0; i<_order; i++)
	{
		solution[i]-=local_dx[i];
	}
}

template<class real>
void LU_Decomposition<real>:: SolvePack(
	autoarray<basearray<real> >& resultpack,
	const autoarray<basearray<real> >& rhspack
	) const
{	// Similiar to solve, but solves a pack of rhs vectors at once, doing its job faster.
	// Vectors are stored in columns of rhspack.
	assert(_decomposed==true);
	const int packsize=rhspack.size();
	int p;
	int i;
	for(i=0; i<_order; i++)
	{
		const int ri=_rhs_rowindex[i];
		for(p=0; p<packsize; p++)
		{
			assert(rhspack[p].size()==_order);
			resultpack[p][i]=rhspack[p][ri];
		}
	}
	for(i=0; i<_order; i++)
	{
		for(p=0; p<packsize; p++)
		{
			basearray<real>& resultp=resultpack[p];
			const real temp=resultp[i];
			if(temp!=0)
			{
				const real * const dat=&_data[i];
				const int jend=_lastnz[i];
				for(int j=i+1; j<=jend; j++)
					resultp[j]-=dat[_rowindex[j]]*temp;
			}
		}
	}
	for(i=_order-1; i>=0; i--)
	{
		for(p=0; p<packsize; p++)
		{
			basearray<real>& resultp=resultpack[p];
			if(resultp[i]!=0)
			{
				const real * const dat=&_data[i];
				resultp[i]/=dat[_rowindex[i]];
				const real temp=resultp[i];
				const int jend=_firstnz[i];
				for(int j=i-1; j>=jend; j--)
					resultp[j]-=dat[_rowindex[j]]*temp;
			}
		}
	}
}

// Test functions, make sure you got desired accuracy: //////////////////

template<class real>
bool LU_Decomposition<real>:: IsSymmetric() const
{
	for(int y=0; y<_order; y++)
	{
		for(int x=y+1; x<_order; x++)
		{
			if(operator()(x, y)!=operator()(y, x))
				return false;
		}
	}
	return true;
}
template<class real>
void LU_Decomposition<real>:: ReverseDecompose(LU_Decomposition& target) const
{	// Could be much faster without L(x,y) and U(x,y), but is only for testing purposes
	assert(_decomposed==true);
	assert(_order==target._order);
	for(int y=0; y<_order; y++)
	{
		for(int x=0; x<_order; x++)
		{
			real sum=0;
			for(int i=0; i<_order; i++)
				sum+=L(i, y)*U(x, i);
			target(x, _rhs_rowindex[y])=(real)sum;
		}
	}
}
template<class real>
void LU_Decomposition<real>:: Compare(const LU_Decomposition& error_reference) const
{	// Compares all matrix elements with error_reference matrix
	// f.ex to be called after LU...ReverseLU with original matrix
	// as parameter to get LU error estimation
	assert(_order==error_reference._order);
	real max_local_error=0;
	real sumpow2=0;
	for(int i=0; i<_orderpow2; i++)
	{
		real local_error=_data[i]-error_reference._data[i];
		if(local_error<0)
			local_error=-local_error;
		if(local_error>max_local_error)
			max_local_error=local_error;
		sumpow2+=local_error*local_error;
	}
	const double mean_square_error=
		sqrt( static_cast<double>(sumpow2/_orderpow2) );
	cout<<"Mean square error="<<mean_square_error<<", max. local error="<<max_local_error<<endl;
}
template<class real>
real LU_Decomposition<real>:: TestAnySolution(const basearray<real>& result,
	const basearray<real>& rhs,
	real eps // =0
	) const
{	// Shows accuracy of any given solution.
	// Should be called with non-LUdecomposed matrix.
	// If max. local error is smaller than eps, there is no screen output.
	// Use eps<0 to disable screen output.
	assert(_decomposed==false);
	assert(result.size()==_order);
	assert(rhs.size()==_order);
	assert(eps>=0);
	cout<<"Testing given result...\n";
	cout.flush();

	basearray<real> local_rhs(_order);
	int y;
	for(y=0; y<_order; y++)
	{
		real mulsum=0;
		for(int x=0; x<_order; x++)
			mulsum+=operator()(x, y)*result[x];
		local_rhs[y]=mulsum;
	}
	real max_local_error=0;
	real sumpow2=0;
	for(y=0; y<_order; y++)
	{
		real local_error=rhs[y]-local_rhs[y];
		if(local_error<0)
			local_error=-local_error;
		if(local_error>max_local_error)
			max_local_error=local_error;
		sumpow2+=local_error*local_error;
	}
	if(eps>=0 && max_local_error>=eps)
	{
		const double mean_square_error=sqrt( static_cast<double>(sumpow2/_order) );
		cout<<" TestSolve: Mean square error="<<mean_square_error<<", max. local error="<<max_local_error<<endl;
	}
	return sqrt( static_cast<double>(sumpow2/_order) );
}
template<class real>
void LU_Decomposition<real>:: TestSolve(const basearray<real>& rhs, real eps) const
{	// Shows accuracy of Solve member function only for given right-hand side vector.
	// Should be called with non-LUdecomposed matrix.
	// If max. local error is smaller than eps, there is no screen output.
	// Use eps<0 to disable screen output.
	assert(_decomposed==false);
	assert(rhs.size()==_order);
	assert(eps>=0);
	cout<<"Testing 'Solve'...\n";
	cout.flush();
	LU_Decomposition<real> copy(*this);
	copy.Decompose();
	basearray<real> test_result(rhs.size());
	copy.Solve(test_result, rhs);
	basearray<real> test_rhs(_order);
	int y;
	for(y=0; y<_order; y++)
	{
		real mulsum=0;
		for(int x=0; x<_order; x++)
			mulsum+=operator()(x, y)*test_result[x];
		test_rhs[y]=mulsum;
	}
	real max_local_error=0;
	real sumpow2=0;
	for(y=0; y<_order; y++)
	{
		real local_error=rhs[y]-test_rhs[y];
		if(local_error<0)
			local_error=-local_error;
		if(local_error>max_local_error)
			max_local_error=local_error;
		sumpow2+=local_error*local_error;
	}
	if(eps>=0 && max_local_error>=eps)
	{
		const double mean_square_error=
			sqrt( static_cast<double>(sumpow2/_order) );
		cout<<" TestSolve: Mean square error="<<mean_square_error<<", max. local error="<<max_local_error<<endl;
	}
}
template<class real>
void LU_Decomposition<real>:: TestSolvePack(
	const autoarray<basearray<real> >& rhspack,
	real eps
	) const
{	// Shows accuracy of SolvePack member function only for given right-hand side pack
	// of vectors.
	// Should be called with non-LUdecomposed matrix.
	// If max. local error is smaller than eps, there is no screen output.
	// Use eps<0 to disable screen output.
	assert(_decomposed==false);
	assert(rhspack[0].size()==_order);
	assert(eps>=0);
	cout<<"Testing 'SolvePack'...\n";
	cout.flush();
	LU_Decomposition<real> copy(*this);
	copy.Decompose();
	autoarray<basearray<real> > test_result(rhspack.size());
	int y;
	for(y=0; y<rhspack.size(); y++)
		test_result[y].setsize(rhspack[0].size());
	copy.SolvePack(test_result, rhspack);

	autoarray<basearray<real> > test_rhs(rhspack.size());
	for(y=0; y<rhspack.size(); y++)
		test_rhs[y].setsize(rhspack[0].size());

	int p;
	for(p=0; p<rhspack.size(); p++)
	{
		for(y=0; y<_order; y++)
		{
			real mulsum=0;
			for(int x=0; x<_order; x++)
				mulsum+=operator()(x, y)*test_result[p][x];
			test_rhs[p][y]=mulsum;
		}
	}

	real max_local_error=0;
	real sumpow2=0;
	for(p=0; p<rhspack.size(); p++)
	{
		for(y=0; y<_order; y++)
		{
			real local_error=rhspack[p][y]-test_rhs[p][y];
			if(local_error<0)
				local_error=-local_error;
			if(local_error>max_local_error)
				max_local_error=local_error;
			sumpow2+=local_error*local_error;
		}
	}
	if(eps>=0 && max_local_error>=eps)
	{
		const double mean_square_error=
			sqrt( static_cast<double>( sumpow2/(_order*rhspack.size()) ) );
		cout<<" TestSolvePack: Mean square error="<<mean_square_error<<", max. local error="<<max_local_error<<endl;
	}
}
template<class real>
void LU_Decomposition<real>:: TestDecompose() const
{	// Shows accuracy of LU decomposition and Determinant, memorized null elements.
	// Should be called with non-LUdecomposed matrix.
	assert(_decomposed==false);
	cout<<"Testing 'Decompose'...\n";
	cout.flush();
	LU_Decomposition<real> copy(*this);
	copy.Decompose();
	int unwanted_null_elements=0;
	for(int y=0; y<_order; y++)
		for(int x=0; x<_order; x++)
			if(_data[x+_rowindex[y]]==0)
				if(y>=copy._firstnz[x] && y<=copy._lastnz[x])
					unwanted_null_elements++;
	cout<<" After LU there is "<<unwanted_null_elements<<" unwanted null elements (";
	cout<<100.*unwanted_null_elements/(_order*_order)<<"%)"<<" of total "<<_order*_order<<".\n";

	const real det=copy.Determinant();
	cout<<" Determinant: "<<det<<endl;
	LU_Decomposition<real> lu_reversed(_order);
	copy.ReverseDecompose(lu_reversed);
	cout<<" LU test: ";
	Compare(lu_reversed);
	cout<<" LU swapped rows: "<<_swap_rows<<endl;
}
template<class real>
void LU_Decomposition<real>:: TestMatrix() const
{	// Tries to classify the matrix by its real, diplays matrix statistics
	assert(_decomposed==false);
	cout<<"Testing matrix...\n";
	cout.flush();
	// Order, memory usage
	// Has nonnull diagonal?
	// Number/percent of nonnull elements (is sparse?)
	// Matrix elements statistics: min, max, minabs, maxabs, median, variance
	// Symmetric placement
	// Symmetric values
	// Condition
	/*
	//TODO: Positive definite?
	//TODO: Diagonally dominant by row?
	//TODO: Diagonally dominant by column?
	//TODO: Diagonally dominant (by row and column)?
	*/
	cout<<" Order: "<<_order<<", n. of elements: "<<_order*_order<<", memory usage: "<<MemoryUsage()<<"B or ";
	cout<<MemoryUsage()/1024.0<<"KB\n";

	int diagonal_null_elements=0;
	int i;
	for(i=0; i<_order; i++)
		if(operator()(i, i)==0)
			diagonal_null_elements++;
	cout<<" Diagonal null elements: "<<diagonal_null_elements<<endl;
	int nonnull=0;
	real d0abs=(_data[0]>=0 ? _data[0] : -_data[0]);
	real min=_data[0];
	real max=_data[0];
	real minabs=d0abs;
	real maxabs=d0abs;
	real sum=_data[0];
	real sumabs=d0abs;
	real sumpow2=_data[0]*_data[0];
	for(i=0; i<_orderpow2; i++)
	{
		const real dat=_data[i];
		if(dat!=0)
			nonnull++;
		sum+=dat;
		sumpow2+=dat*dat;
		if(dat>max)
			max=dat;
		if(dat<min)
			min=dat;
		const real databs=(dat>=0 ? dat : -dat);
		sumabs+=databs;
		if(databs>maxabs)
			maxabs=databs;
		if(databs<minabs)
			minabs=databs;
	}
	const real average=sum/_orderpow2;
	const real averageabs=sumabs/_orderpow2;
	const real variance=sumpow2/_orderpow2 - average*average;
	assert(variance>=0);
	cout<<" Nonzero elements: "<<nonnull<<" of "<<_orderpow2<<" ("<<100.0*nonnull/_orderpow2<<"%)\n";
	cout<<" Matrix elements statistics: min="<<min<<" max="<<max<<endl;
	cout<<"						  minabs="<<minabs<<" maxabs="<<maxabs<<endl;
	cout<<" average="<<average<<" averageabs="<<averageabs<<" variance="<<variance<<endl;
	bool symmetric_nonzero_elements_matrix=true;
	bool symmetric_matrix=true;
	for(int y=0; y<_order; y++)
	{
		for(int x=y+1; x<_order; x++)
		{
			if(operator()(x, y)==0 && operator()(y, x)!=0)
			{
				symmetric_nonzero_elements_matrix=false;
				symmetric_matrix=false;
				break;
			}
			if(operator()(x, y)!=0 && operator()(y, x)==0)
			{
				symmetric_nonzero_elements_matrix=false;
				symmetric_matrix=false;
				break;
			}
			if(operator()(x, y)!=operator()(y, x))
				symmetric_matrix=false;
		}
	}
	cout<<" Matrix symmetry: ";
	if(symmetric_nonzero_elements_matrix==true)
		cout<<"has sym. nonnull elements, ";
	else
		cout<<"hasn't sym. nonnull elements, ";
	if(symmetric_matrix==true)
		cout<<"is sym. (A_transp==A)\n";
	else
		cout<<"is not sym. (A_transp!=A)\n";

	cout<<" Matrix condition: ";
	cout.flush();
	cout<<Condition()<<endl;
}
#ifdef _INCLUDED_L3_IMAGE_H
template<class real>
void LU_Decomposition<real>:: SaveGIF(const TextString& gif_filename) const
{	// Saves image of nonnull elements of matrix to GIF file.
	// One must include l3_image.h manually before using this function.
	// l3_image.h may not be freely available now.
	RGBImage img(_order+2, _order+2);
	img.Clear();
	img.SetColor(.5, .8, .8);
	img.HLineRel(0, 0, _order+2);
	img.VLineRel(0, 0, _order+2);
	img.VLineRel(_order+1, 0, _order+2);
	img.HLineRel(0, _order+1, _order+2);
	real minvalue=(real)fabs(_data[0]);
	real maxvalue=minvalue;
	for(int i=1; i<_orderpow2; i++)
	{
		const real testvalue=(real)fabs(_data[i]);
		if(testvalue<minvalue)
			minvalue=testvalue;
		else if(testvalue>maxvalue)
			maxvalue=testvalue;
	}
	const real rescale=1/(maxvalue-minvalue);

	for(int y=0; y<_order; y++)
	{
		for(int x=0; x<_order; x++)
		{
			const real value=_data[x+_rowindex[y]];
			if(value!=0)
			{
				const real col=static_cast<real>(.8-fabs(value)*rescale*0.8);
				assert(col>=0);
				assert(col<=1);
				img.SetColor(1, col, col);
			}
			else if(_decomposed==true && y>=_firstnz[x] && y<=_lastnz[x])
				img.SetColor(.8, .9, .9);
			else
				img.SetColor(1, 1, 1);
			img.RawPoint(x+1, y+1);
		}
	}
	img.SaveGIF(gif_filename, 0);
}
#endif

/////////////////////////////////////////////////////////////////////////////
/*
template<class real>
void ExplicitScaling(
			basearray<real>& solution,
			LU_Decomposition<real> mat,
			basearray<real> rhs
			)
{	// Scales a system (before LU decomposition) so that row norms would be small,
	// without (!) introducing rounding errors. Could improve matrix condition considerably.
	for(int y=0; y<mat.Order(); y++)
	{
		real maxrow=mat(0, y);
		if(maxrow<0)
			maxrow=-maxrow;
		int x;
		for(x=1; x<mat.Order(); x++)
		{
			real temp=mat(x, y);
			if(temp<0)
				temp=-temp;
			if(temp>maxrow)
				maxrow=temp;
		}
		//const real scaling_factor=1./maxrow;// intuitive
		const real scaling_factor=// no rounding errors, number close to intuitive approach
			exp(-static_cast<int>(1+log(maxrow)/log(2)) * log(2));
		assert(scaling_factor>0);
		//cout<<scaling_factor<<"|";
		for(x=0; x<mat.Order(); x++)
			mat(x, y)*=scaling_factor;
		rhs[y]*=scaling_factor;
	}
	mat.Solve(solution, rhs);
}
*/

/////////////////////////////////////////////////////////////////////////////

#endif // _INCLUDED_L3_MAT_H