// Matrix Template Class  ver.0.45
// Copyright 1998-2008, Foota Software (Noriyuki Futatsugi)

#ifndef TMATRIX_HEADER
#define TMATRIX_HEADER

#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>
#include <vector>
#include <cmath>

namespace tmatrix {
const double EPS      = 1e-12;  // maximum permissible error of non-diagonal component
const int    MAX_ITER = 100;    // maximum iterations
}

template<class T> class TMatrix;

/////////////////////////////////////////////////////////////////////
// vector template class
template<class T>
class TVector {
private:
	T *m_data;
	int m_num;

public:
	TVector(int num = 1) : m_data(new T[num]), m_num(num) { }
	TVector(const std::vector<T>&);
//	TVector(const TVector&);
	TVector(const TVector<T>&);
	~TVector() { delete [] m_data; }

	T& operator[](int pos) { if (pos < 0) pos = 0;  return m_data[pos]; }
	const T& operator[](int pos) const { return m_data[pos]; }
	T& operator()(int pos) { if (pos < 0) pos = 0;  return m_data[pos]; }
	const T& operator()(int pos) const { return m_data[pos]; }
	void Set(int pos, T val = 0) { m_data[pos] = val; }
	T Get(const int& pos) const { return m_data[pos]; }
	void SetArray(std::vector<T>& vec) { *this = vec; }
	std::vector<T> GetArray();
	void Print(int width = 3, int pre = 0, int col = 0, int base = 1);
	void Resize(int);
	int GetNum() const { return m_num; }

	TVector& operator=(const TVector<T>&);
	TVector& operator=(const std::vector<T>&);
	TVector& operator=(const T);
	TVector& operator=(const int);
//	TVector& operator=(const char*);
	TVector& operator=(const std::string&);
	TVector& operator=(const TMatrix<T>&);

	T Householder();
};

// constructor (std::std::vector<T>)
template<class T>
TVector<T>::TVector(const std::vector<T>& vec) : m_num(vec.size())
{
	m_data = new T[m_num];
	for (int i = 0; i < m_num; i++) m_data[i] = vec[i];
}

// copy constructor
template<class T>
TVector<T>::TVector(const TVector<T>& vec) : m_num(vec.GetNum())
{
	m_data = new T[m_num];
	for (int i = 0; i < m_num; i++) m_data[i] = vec[i];
}

// Get std::vector<T> format
template<class T>
std::vector<T> TVector<T>::GetArray()
{
	std::vector<T> vec;

	for (int i = 0; i < m_num; i++) vec.push_back(m_data[i]);

	return vec;
}

// output TVector
template<class T>
void TVector<T>::Print(int width, int pre, int col, int base)
{
	std::ios_base::fmtflags ios_org;
	int pre_org;
	ios_org = std::cout.setf(std::ios_base::fixed | std::ios_base::right);
	pre_org = std::cout.precision(pre);

	if (col > 0) {
		for (int j = 0; j < m_num; j += col) {
			for (int i = j; i < j + col && i < m_num; i++)
				std::cout << std::setw(width) << i + base;
			std::cout << std::endl;
			for (int i = j; i < j + col && i < m_num; i++) {
				std::cout << std::setw(width) << m_data[i];
				if ((i + 1) % col == 0 && (i + 1) != m_num) std::cout << std::endl;
			}
			std::cout << std::endl;
		}
	} else if (col < 0) {
		col = -col;
		for (int i = 0; i < m_num; i++) {
			std::cout << std::setw(width) << m_data[i];
			if ((i + 1) % col == 0 && (i + 1) != m_num) std::cout << std::endl;
		}
		std::cout << std::endl;
	} else {
		for (int i = 0; i < m_num; i++) std::cout << std::setw(width) << m_data[i];
		std::cout << std::endl;
	}

	std::cout.flags(ios_org);
	std::cout.precision(pre_org);
}

// resize memories
template<class T>
void TVector<T>::Resize(int num)
{
	delete [] m_data;

	m_data = new T[num];
	m_num = num;
}

// substitute TVector
template<class T>
TVector<T>& TVector<T>::operator=(const TVector<T>& vec)
{
	if (&vec == this) return *this;

	if (m_num != vec.GetNum()) {
		delete [] m_data;
		m_num = vec.GetNum();
		m_data = new T[m_num];
	}

	for (int i = 0; i < m_num; i++) m_data[i] = vec[i];

	return *this;
}

// substitute std::vector
template<class T>
TVector<T>& TVector<T>::operator=(const std::vector<T>& vec)
{
	if (m_num != static_cast<int>(vec.size())) {
		delete [] m_data;
		m_num = static_cast<int>(vec.size());
		m_data = new T[m_num];
	}

	for (int i = 0; i < m_num; i++) m_data[i] = vec[i];

	return *this;
}

// substitute T
template<class T>
TVector<T>& TVector<T>::operator=(const T val)
{
	for (int i = 0; i < m_num; i++) m_data[i] = val;

	return *this;
}

// substitute int
template<class T>
TVector<T>& TVector<T>::operator=(const int val)
{
	return operator=(static_cast<T>(val));
}

// substitute std::string
template<class T>
TVector<T>& TVector<T>::operator=(const std::string& data)
{
	std::stringstream ss;
	double value;
	std::vector<T> vec;

	std::string::size_type index_first = 0;
	std::string::size_type index_last = 0;
	for (int i = 0; i < m_num; i++) {
		index_first = data.find_first_of("-+.01234567890", index_last);
		index_last = data.find_first_not_of("-+.01234567890", index_first);
		ss.str(data.substr(index_first, index_last - index_first));
		ss >> value;
		ss.clear();
		vec.push_back(value);
	}
	if (m_num != vec.size()) {
		std::cerr << "Error: illegal format in TMatrix." << std::endl;
		exit(1);
	}
	for (std::string::size_type i = 0; i < vec.size(); i++) m_data[i] = vec[i];

	return *this;
}

/*
template<class T>
TVector<T>& TVector<T>::operator=(const char *data)
{
//	char buf[m_num][50];
	char **buf;
	char *ptr;
	int i;

	buf = new char*[m_num];
	for (i = 0; i < m_num; i++) buf[i] = new char[50];
	if ((ptr = strtok(const_cast<char*>(data), "\n\r\t, ")) != NULL)
		strcpy(buf[0], ptr);
	for (i = 0; i < m_num; i++)
		if ((ptr = strtok(NULL, "\n\r\t, ")) != NULL) strcpy(buf[i+1], ptr);

	for (i = 0; i < m_num; i++) m_data[i] = atof(buf[i]);

	for (i = 0; i < m_num; i++) delete [] buf[i];
	delete [] buf;

	return *this;
}
*/

// substitute diagonal component of TMatrix
template<class T>
TVector<T>& TVector<T>::operator=(const TMatrix<T>& mat)
{
	if (mat.GetRow() != mat.GetCol()) return *this;

	if (m_num != mat.GetRow()) {
		delete [] m_data;
		m_num = mat.GetRow();
		m_data = new T[m_num];
	}

	for (int i = 0; i < m_num; i++) m_data[i] = mat[i][i];

	return *this;
}

// summation of TVector
template<class T>
TVector<T> operator+(const TVector<T>& x, const TVector<T>& y)
{
	TVector<T> z(x.GetNum());

	for (int i = 0; i < z.GetNum(); i++) z[i] = x[i] + y[i];

	return z;
}

// subtraction of TVector
template<class T>
TVector<T> operator-(const TVector<T>& x, const TVector<T>& y)
{
	TVector<T> z(x.GetNum());

	for (int i = 0; i < z.GetNum(); i++) z[i] = x[i] - y[i];

	return z;
}

// product of TVector
template<class T>
T operator*(const TVector<T>& x, const TVector<T>& y)
{
	T s;
	int n5;

	if (x.GetNum() != y.GetNum()) return 0;

	s = 0;  n5 = x.GetNum() % 5;
	for (int i = 0; i < n5; i++) s += x[i] * y[i];
	for (int i = n5; i < x.GetNum(); i += 5)
		s += x[i] * y[i] + x[i+1] * y[i+1] + x[i+2] * y[i+2]
			+ x[i+3] * y[i+3] + x[i+4] * y[i+4];

	return s;
}

/////////////////////////////////////////////////////////////////////
// matrix template class
template<class T>
class TMatrix {
private:
	T *m_data;
	int m_nrow, m_ncol;

public:
	template<class U>
	class CProxyMatrix {
	private:
		U *m_px_data;
		int m_px_col;

	public:
		CProxyMatrix(U *p, int col) : m_px_data(p), m_px_col(col) { }

		U& operator[](int index) { return m_px_data[index]; }
		const U& operator[](int index) const
			{ return m_px_data[index]; }
	};

	TMatrix(int num = 1)
		: m_data(new T[num*num]), m_nrow(num), m_ncol(num) { }
	TMatrix(int row, int col)
		: m_data(new T[row*col]), m_nrow(row), m_ncol(col) { }
	TMatrix(const std::vector<std::vector<T> >&);
//	TMatrix(const TMatrix&);
	TMatrix(const TMatrix<T>&);
	~TMatrix() { delete [] m_data; }

	CProxyMatrix<T> operator[](int row)
		{ return CProxyMatrix<T>(m_data + row * m_ncol, m_ncol); }
	const CProxyMatrix<T> operator[](int row) const
		{ return CProxyMatrix<T>(m_data + row * m_ncol, m_ncol); }

	T& operator()(int row, int col) { return m_data[row*m_ncol+col]; }
	const T& operator()(int row, int col) const
		{ return m_data[row*m_ncol+col]; }
	void Set(int row, int col, T num = 0)
		{ m_data[row*m_ncol+col] = num; }
	T Get(const int& row, const int& col) const
		{ return m_data[row*m_ncol+col]; }
	void SetArray(std::vector<std::vector<T> >& vec) { *this = vec; }
	std::vector<std::vector<T> > GetArray();
	void Print(int width = 3, int pre = 0, int col = 0, int base = 1);
	void Resize(int num) { Resize(num, num); }
	void Resize(int, int);
	int GetRow() const { return m_nrow; }
	int GetCol() const { return m_ncol; }
//	TVector<T> SliceRow(int);
//	TVector<T> SliceRow(int) const;
//	TVector<T> SliceCol(int);
//	TVector<T> SliceCol(int) const;
	void InsertRow(int, const TVector<T>&);
	void InsertCol(int, const TVector<T>&);
	std::vector<T> SliceRow(int);
	std::vector<T> SliceRow(int) const;
	std::vector<T> SliceCol(int);
	std::vector<T> SliceCol(int) const;
	void InsertRow(int, const std::vector<T>&);
	void InsertCol(int, const std::vector<T>&);
	TMatrix& LU(T&, int*);
	TMatrix& LUInverse(TMatrix&, T&);
	TMatrix& SetInverse(T&);
	TMatrix Inverse(T&);
	TMatrix Inverse();
	T Determinant();
	TMatrix Trans();
	T Trace();

	TMatrix& operator=(const TMatrix<T>&);
	TMatrix& operator=(const T);
	TMatrix& operator=(const int);
//	TMatrix& operator=(const char*);
	TMatrix& operator=(const std::string&);
	TMatrix& operator=(const TVector<T>&);

	TMatrix Diag(TVector<T>&);
	TMatrix& Tridiagonalize(TVector<T>&, TVector<T>&);
	TMatrix& Eigen(TVector<T>&);
};

// constructor (std::std::vector<std::std::vector<T> >)
template<class T>
TMatrix<T>::TMatrix(const std::vector<std::vector<T> >& mat)
	: m_nrow(mat.size()), m_ncol(mat[0].size())
{
	m_data = new T[m_nrow*m_ncol];

	for (int i = 0; i < m_nrow; i++)
		for (int j = 0; j < m_ncol; j++)
			m_data[i*m_ncol+j] = mat[i][j];
}

// copy constructor
template<class T>
TMatrix<T>::TMatrix(const TMatrix<T>& mat)
	: m_nrow(mat.GetRow()), m_ncol(mat.GetCol())
{
	m_data = new T[m_nrow*m_ncol];

	for (int i = 0; i < m_nrow; i++)
		for (int j = 0; j < m_ncol; j++)
			m_data[i*m_ncol+j] = mat[i][j];
}

// get std::vector<std::vector<T> > format
template<class T>
std::vector<std::vector<T> > TMatrix<T>::GetArray()
{
	std::vector<std::vector<T> > vec(m_nrow);
	for (int i = 0; i < m_nrow; i++) {
		vec[i].resize(m_ncol);
		for (int j = 0; j < m_ncol; j++) vec[i].push_back(m_data[i*m_ncol+j]);
	}

	return vec;
}

// output TMatrix
// width: width of value
// pre:   precision of value
// col:   column of value
// base:  basement of adscript of value
template<class T>
void TMatrix<T>::Print(int width, int pre, int col, int base)
{
	std::ios_base::fmtflags ios_org;
	int pre_org;
	ios_org = std::cout.setf(std::ios_base::fixed | std::ios_base::right);
	pre_org = std::cout.precision(pre);
	const int width_col = 5;

	if (col > 0) {
		for (int k = 0; k < m_ncol; k += col) {
			std::cout << std::setw(width_col) << "";
			for (int i = k; i < k + col && i < m_ncol; i++)
				std::cout << std::setw(width) << i + base;
			std::cout << std::endl;
			for (int i = 0; i < m_nrow; i++) {
				std::cout << std::setw(width_col) << i + base;
				for (int j = k; j < m_ncol; j++) {
					std::cout << std::setw(width) << m_data[i*m_ncol+j];
					if ((j + 1) % col == 0) break;
				}
				std::cout << std::endl;
			}
			if (k + col < m_ncol) std::cout << std::endl;
		}
	} else if (col < 0) {
		col = -col;
		for (int i = 0; i < m_nrow; i++) {
			for (int j = 0; j < m_ncol; j++) {
				std::cout << std::setw(width) << m_data[i*m_ncol+j];
				if ((j + 1) % col == 0 || (j + 1) == m_ncol) std::cout << std::endl;
			}
			if (i + 1 < m_nrow) std::cout << std::endl;
		}
	} else {
		for (int i = 0; i < m_nrow; i++) {
			for (int j = 0; j < m_ncol; j++)
				std::cout << std::setw(width) << m_data[i*m_ncol+j];
			std::cout << std::endl;
		}
	}

	std::cout.flags(ios_org);
	std::cout.precision(pre_org);
}

// resize memories
template<class T>
void TMatrix<T>::Resize(int row, int col)
{
	delete [] m_data;

	m_data = new T[row*col];
	m_nrow = row;
	m_ncol = col;
}

// extract any rows
template<class T>
std::vector<T> TMatrix<T>::SliceRow(int row)
{
	std::vector<T> vec(m_ncol);

	for (int i = 0; i < m_ncol; i++) vec[i] = m_data[row*m_ncol+i];

	return vec;
}

// extract any columns
template<class T>
std::vector<T> TMatrix<T>::SliceCol(int col)
{
	std::vector<T> vec(m_nrow);

	for (int i = 0; i < m_nrow; i++) vec[i] = m_data[i*m_ncol+col];

	return vec;
}

// insert any rows
template<class T>
void TMatrix<T>::InsertRow(int row, const TVector<T>& vec)
{
	for (int i = 0; i < m_ncol; i++) m_data[row*m_ncol+i] = vec[i];
}

template<class T>
void TMatrix<T>::InsertRow(int row, const std::vector<T>& vec)
{
	for (int i = 0; i < m_ncol; i++) m_data[row*m_ncol+i] = vec[i];
}

// insert any columns
template<class T>
void TMatrix<T>::InsertCol(int col, const TVector<T>& vec)
{
	for (int i = 0; i < m_nrow; i++) m_data[i*m_ncol+col] = vec[i];
}

template<class T>
void TMatrix<T>::InsertCol(int col, const std::vector<T>& vec)
{
	for (int i = 0; i < m_nrow; i++) m_data[i*m_ncol+col] = vec[i];
}

// substitute TMatrix
template<class T>
TMatrix<T>& TMatrix<T>::operator=(const TMatrix<T>& mat)
{
	if (&mat == this) return *this;

	if (m_nrow != mat.GetRow() || m_ncol != mat.GetCol()) {
		delete [] m_data;
		m_nrow = mat.GetRow();
		m_ncol = mat.GetCol();
		m_data = new T[m_nrow*m_ncol];
	}

	for (int i = 0; i < m_nrow; i++)
		for (int j = 0; j < m_ncol; j++)
			m_data[i*m_ncol+j] = mat[i][j];

	return *this;
}

// substitute T
template<class T>
TMatrix<T>& TMatrix<T>::operator=(const T num)
{
	for (int i = 0; i < m_nrow; i++)
		for (int j = 0; j < m_ncol; j++)
			m_data[i*m_ncol+j] = num;

	return *this;
}

// substitute int
template<class T>
TMatrix<T>& TMatrix<T>::operator=(const int num)
{
	return operator=(T(num));
}

// substitute std::string
template<class T>
TMatrix<T>& TMatrix<T>::operator=(const std::string& data)
{
	std::stringstream ss;
	double value;
	std::vector<T> vec;

	std::string::size_type index_first = 0;
	std::string::size_type index_last = 0;
	for (int i = 0; i < m_nrow * m_ncol; i++) {
		index_first = data.find_first_of("-+.01234567890", index_last);
		index_last = data.find_first_not_of("-+.01234567890", index_first);
		ss.str(data.substr(index_first, index_last - index_first));
		ss >> value;
		ss.clear();
		vec.push_back(value);
	}
	if (m_nrow * m_ncol != vec.size()) {
		std::cerr << "Error: illegal format in TMatrix." << std::endl;
		exit(1);
	}
	for (std::string::size_type i = 0; i < vec.size(); i++) m_data[i] = vec[i];

	return *this;
}

/*
template<class T>
TMatrix<T>& TMatrix<T>::operator=(const char *data)
{
//	char buf[m_nrow*m_ncol][50];
	char **buf;
	char *ptr;
	int i, j;

	buf = new char*[m_nrow*m_ncol];
	for (i = 0; i < m_nrow * m_ncol; i++) buf[i] = new char[50];
	if ((ptr = strtok(const_cast<char*>(data), "\n\r\t, ")) != NULL)
		strcpy(buf[0], ptr);

	for (i = 0; i < m_nrow; i++)
		for (j = 0; j < m_ncol; j++)
			if ((ptr = strtok(NULL, "\n\r\t, ")) != NULL)
				strcpy(buf[i*m_ncol+j+1], ptr);
	for (i = 0; i < m_nrow; i++)
		for (j = 0; j < m_ncol; j++)
			m_data[i*m_ncol+j] = atof(buf[i*m_ncol+j]);

	for (i = 0; i < m_nrow * m_ncol; i++) delete [] buf[i];
	delete [] buf;

	return *this;
}
*/

// substitute TVector (diagonal component of square matrix)
template<class T>
TMatrix<T>& TMatrix<T>::operator=(const TVector<T>& vec)
{
	if (m_nrow != m_ncol || m_nrow != vec.GetNum()) {
		delete [] m_data;
		m_nrow = m_ncol = vec.GetNum();
		m_data = new T[m_nrow*m_ncol];
	}

	for (int i = 0; i < m_nrow; i++) {
		for (int j = 0; j < m_ncol; j++) m_data[i*m_ncol+j] = 0;
		m_data[i*m_ncol+i] = vec[i];
	}

	return *this;
}

// summation of TMatrix
template<class T>
TMatrix<T> operator+(const TMatrix<T>& x, const TMatrix<T>& y)
{
	TMatrix<T> z(x.GetRow(), x.GetCol());

	for (int i = 0; i < z.GetRow(); i++)
		for (int j = 0; j < z.GetCol(); j++)
			z[i][j] = x[i][j] + y[i][j];

	return z;
}

// subtraction of TMatrix
template<class T>
TMatrix<T> operator-(const TMatrix<T>& x, const TMatrix<T>& y)
{
	TMatrix<T> z(x.GetRow(), x.GetCol());

	for (int i = 0; i < z.GetRow(); i++)
		for (int j = 0; j < z.GetCol(); j++)
			z[i][j] = x[i][j] - y[i][j];

	return z;
}

// product of TMatrix
template<class T>
TMatrix<T> operator*(const TMatrix<T>& x, const TMatrix<T>& y)
{
	TMatrix<T> z(x.GetRow(), y.GetCol());
	T s;

	for (int i = 0; i < z.GetRow(); i++) {
		for (int j = 0; j < z.GetCol(); j++) {
			s = 0;
			for (int k = 0; k < x.GetCol(); k++) s += x[i][k] * y[k][j];
			z[i][j] = s;
		}
	}

	return z;
}

// multiply TMatrix by TVector
template<class T>
TVector<T> operator*(const TMatrix<T>& x, const TVector<T>& y)
{
	TVector<T> z(x.GetCol());
	T s;

	for (int i = 0; i < x.GetRow(); i++) {
		s = 0;
		for (int j = 0; j < x.GetCol(); j++) s += x[i][j] * y[j];
		z[i] = s;
	}

	return z;
}

// ------------------------------------------------------------
// algorithms for matrix
// ------------------------------------------------------------

// inverse matrix and determinant (simple version)
template<class T>
TMatrix<T>& TMatrix<T>::SetInverse(T& det)
{
	T t, u;

	if (m_nrow != m_ncol) return *this;

	det = 1;
	for (int k = 0; k < m_nrow; k++) {
		t = m_data[k*m_ncol+k];  det *= t;
		if (t == 0) return *this = 0;  // determinant is 0.
		for (int i = 0; i < m_nrow; i++) m_data[k*m_ncol+i] /= t;
		m_data[k*m_ncol+k] = 1 / t;
		for (int j = 0; j < m_nrow; j++) {
			if (j != k) {
				u = m_data[j*m_ncol+k];
				for (int i = 0; i < m_nrow; i++) {
					if (i != k) m_data[j*m_ncol+i] -= m_data[k*m_ncol+i] * u;
					else	   m_data[j*m_ncol+i] = -u / t;
				}
			}
		}
	}

	return *this;
}

// LU decomposition
template<class T>
TMatrix<T>& TMatrix<T>::LU(T& det, int* ip)
{
	T t, u;
	int j, ii, ik;
//	int i, j, k, ii, ik;

	if (m_nrow != m_ncol) return *this;

	j = 0;
	TVector<T> weight(m_nrow);
	det = 0;  // determinant
	for (int k = 0; k < m_nrow; k++) {  // each rows
		ip[k] = k;  // defaults for information of exchange rows
		u = 0;  // element of max absolute value in the row
		for (int j = 0; j < m_nrow; j++) {
			t = fabs(m_data[k*m_ncol+j]);  if (t > u) u = t;
		}
		if (u == 0) goto EXIT;  // if u == 0, disable LU decomposition
		weight[k] = 1 / u;  // reciprocal number of max absolute value
	}
	det = 1;  // default of determinant
	for (int k = 0; k < m_nrow; k++) {  // each rows
		u = -1;
		for (int i = k; i < m_nrow; i++) {  // above of the row
			ii = ip[i];  // (weight x absolute) to find max row
			t = fabs(m_data[ii*m_ncol+k]) * weight[ii];
			if (t > u) {  u = t;  j = i;  }
		}
		ik = ip[j];
		if (j != k) {
			ip[j] = ip[k];  ip[k] = ik;  // exchange row-numbers
			det = -det;  // change sign to exchange rows
		}
		u = m_data[ik*m_ncol+k];  det *= u;  // diagonal component
		if (u == 0) goto EXIT;  // if u == 0, disable LU decomposition
		for (int i = k + 1; i < m_nrow; i++) {  // Gaussian elimination
			ii = ip[i];
			t = (m_data[ii*m_ncol+k] /= u);
			for (int j = k + 1; j < m_nrow; j++)
				m_data[ii*m_ncol+j] -= t * m_data[ik*m_ncol+j];
		}
	}

  EXIT:

	return *this;
}

// inverse matrix and determinant (LU decomposition)
template<class T>
TMatrix<T>& TMatrix<T>::LUInverse(TMatrix<T>& mat_inv, T& det)
{
	T t;
	int *ip;  // information to exchange rows
	int ii;
//	int i, j, k, ii;

	if (m_nrow != m_ncol) return *this;

	if ((ip = new int[m_nrow]) == NULL) {
		std::cerr << "Error: memory fault (LUInverse)"<< std::endl;
		exit(1);
	}

	LU(det, ip);
	if (det != 0) {
		for (int k = 0; k < m_nrow; k++) {
			for (int i = 0; i < m_nrow; i++) {
				ii = ip[i];  t = (ii == k);
				for (int j = 0; j < i; j++)
					t -= m_data[ii*m_ncol+j] * mat_inv[j][k];
				mat_inv[i][k] = t;
			}
			for (int i = m_nrow - 1; i >= 0; i--) {
				t = mat_inv[i][k];  ii = ip[i];
				for (int j = i + 1; j < m_nrow; j++)
					t -= m_data[ii*m_ncol+j] * mat_inv[j][k];
				mat_inv[i][k] = t / m_data[ii*m_ncol+i];
			}
		}
	}

	delete [] ip;

	return *this;
}

// inverse matrix and determinant
template<class T>
TMatrix<T> TMatrix<T>::Inverse(T& det)
{
	TMatrix<T> mat(*this), mat_inv(m_nrow, m_ncol);

	mat.LUInverse(mat_inv, det);  // using LU decomposition

//	return mat.SetInverse(det);  // (simple version)
	return mat_inv;
}

// inverse matrix
template<class T>
TMatrix<T> TMatrix<T>::Inverse()
{
	T det;
	TMatrix<T> mat(*this), mat_inv(m_nrow, m_ncol);

	mat.LUInverse(mat_inv, det);  // using LU decomposition

	return mat_inv;
}

// determinant
template<class T>
T TMatrix<T>::Determinant()
{
	T det;
	TMatrix<T> mat(*this), mat_inv(m_nrow, m_ncol);

	mat.LUInverse(mat_inv, det);  // using LU decomposition

	return det;
}

// transposed matrix
template<class T>
TMatrix<T> TMatrix<T>::Trans()
{
	TMatrix<T> mat(m_ncol, m_nrow);

	for (int i = 0; i < m_nrow; i++)
		for (int j = 0; j < m_ncol; j++)
			mat[j][i] = m_data[i*m_ncol+j];

	return mat;
}

// diagonal component of square matrix (trace)
template<class T>
T TMatrix<T>::Trace()
{
	T tr = 0;

	if (m_nrow != m_ncol) return tr;

	for (int i = 0; i < m_nrow; i++)
		tr += m_data[i*m_ncol+i];

	return tr;
}

// diagonalization
template<class T>
TMatrix<T> TMatrix<T>::Diag(TVector<T>& d)
{
	TMatrix<T> mat(*this);

	return mat.Eigen(d);  // QR algorithm
}

// Householder transformation
template<class T>
T TVector<T>::Householder()
{
	T s, t;

	s = sqrt(*this * *this);  // squared innner product

	if (s != 0) {
		if (m_data[0] < 0) s = -s;
		m_data[0] += s;  t = 1 / sqrt(m_data[0] * s);
		for (int i = 0; i < m_num; i++) m_data[i] *= t;
	}

	return -s;
}

// tridiagonalization
template<class T>
TMatrix<T>& TMatrix<T>::Tridiagonalize(TVector<T>& d, TVector<T>& e)
{
	T s, t, p, q;
	TVector<T> v, w;
	TVector<T> v_temp, d_temp, w_temp;
//	int i, j, k, l;

	if (m_nrow != m_ncol) return *this;

	for (int k = 0; k < m_nrow - 2; k++) {
		v = SliceRow(k);  d[k] = v[k];
		v_temp.Resize(m_nrow - k - 1);
		for (int i = 0; i < m_nrow - k - 1; i++) v_temp[i] = v[i+k+1];
		e[k] = v_temp.Householder();
		for (int i = 0; i < m_nrow - k - 1; i++) v[i+k+1] = v_temp[i];
		if (e[k] == 0) continue;
		for (int i = k + 1; i < m_nrow; i++) {
			s = 0;
			for (int j = k + 1; j < i; j++) s += m_data[j*m_ncol+i] * v[j];
			for (int j = i; j < m_nrow; j++) s += m_data[i*m_ncol+j] * v[j];
			d[i] = s;
		}
		v_temp.Resize(m_nrow - k - 1);
		d_temp.Resize(m_nrow - k - 1);
		for (int i = 0; i < m_nrow - k - 1; i++) 
			{ v_temp[i] = v[i+k+1];  d_temp[i] = d[i+k+1]; }
		t = v_temp * d_temp / 2;
		for (int i = m_nrow - 1; i > k; i--) {
			p = v[i];  q = d[i] - t * p;  d[i] = q;
			for (int j = i; j < m_nrow; j++)
				m_data[i*m_ncol+j] -= p * d[j] + q * v[j];
		}
		InsertRow(k, v);
	}
	if (m_nrow >= 2) {
		d[m_nrow-2] = m_data[(m_nrow-2)*m_ncol+(m_nrow-2)];
		e[m_nrow-2] = m_data[(m_nrow-2)*m_ncol+(m_nrow-1)];
	}
	if (m_nrow >= 1)   d[m_nrow-1] = m_data[(m_nrow-1)*m_ncol+(m_nrow-1)];
	for (int k = m_nrow - 1; k >= 0; k--) {
		v = SliceRow(k);
		if (k < m_nrow - 2) {
			for (int i = k + 1; i < m_nrow; i++) {
				w = SliceRow(i);
				v_temp.Resize(m_nrow - k - 1);
				w_temp.Resize(m_nrow - k - 1);
				for (int l = 0; l < m_nrow - k - 1; l++) 
				{ v_temp[l] = v[l+k+1];  w_temp[l] = w[l+k+1]; }
				t = v_temp * w_temp;
				for (int j = k + 1; j < m_nrow; j++) w[j] -= t * v[j];
				InsertRow(i, w);
			}
		}
		for (int i = 0; i < m_nrow; i++) v[i] = 0;
		v[k] = 1;
		InsertRow(k, v);
	}

	return *this;
}

// diagonalization (QR algorithm)
template<class T>
TMatrix<T>& TMatrix<T>::Eigen(TVector<T>& d)
{
	T c, s, t, w, x, y;
	TVector<T> v;
	TVector<T> e;
	TVector<T> e_temp;
	int j, iter;
//	int i, j, k, h, iter;

	if (m_nrow != m_ncol) return *this;

	v.Resize(m_nrow);
	e.Resize(m_nrow);
	e_temp.Resize(m_nrow);
	Tridiagonalize(d, e_temp);
	e[0] = 0;
	for (int i = 0; i < m_nrow - 1; i++) e[i+1] = e_temp[i];
	for (int h = m_nrow - 1; h > 0; h--) {
		j = h;
		while (fabs(e[j]) > tmatrix::EPS * (fabs(d[j-1]) + fabs(d[j]))) j--;
		if (j == h) continue;
		iter = 0;
		do {
			if (++iter > tmatrix::MAX_ITER) return *this;
			w = (d[h-1] - d[h]) / 2;
			t = e[h] * e[h];
			s = sqrt(w * w + t);  if (w < 0) s = -s;
			x = d[j] - d[h] + t / (w + s);  y = e[j+1];
			for (int k = j; k < h; k++) {
				if (fabs(x) >= fabs(y)) {
					t = -y / x;  c = 1 / sqrt(t * t + 1);
					s = t * c;
				} else {
					t = -x / y;  s = 1 / sqrt(t * t + 1);
					c = t * s;
				}
				w = d[k] - d[k+1];
				t = (w * s + 2 * c * e[k+1]) * s;
				d[k] -= t;  d[k+1] += t;
				if (k > j) e[k] = c * e[k] - s * y;
				e[k+1] += s * (c * w - 2 * s * e[k+1]);
				for (int i = 0; i < m_nrow; i++) {
					x = m_data[k*m_ncol+i];  y = m_data[(k+1)*m_ncol+i];
					m_data[k*m_ncol+i] = c * x - s * y;
					m_data[(k+1)*m_ncol+i] = s * x + c * y;
				}
				if (k < h - 1) {
					x = e[k+1];  y = -s * e[k+2];
					e[k+2] *= c;
				}
			}
		} while (fabs(e[h]) > tmatrix::EPS * (fabs(d[h-1]) + fabs(d[h])));
	}
	// eigen values and eigen vectors
	// sort by descending order
/*
	int h;
	for (int k = 0; k < m_nrow - 1; k++) {
		h = k;  t = d[h];
		for (int i = k + 1; i < m_nrow; i++)
			if (d[i] < t) { h = i;  t = d[h]; }
		d[h] = d[k];  d[k] = t;
		v = SliceRow(h);
		InsertRow(h, SliceRow(k));
		InsertRow(k, v);
	}
*/
	for (int k = 0; k < m_nrow - 1; k++) {
		j = k;  t = d[j];
		for (int i = k + 1; i < m_nrow; i++)
			if (d[i] > t) { j = i;  t = d[j]; }
		d[j] = d[k];  d[k] = t;
		v = SliceRow(j);
		InsertRow(j, SliceRow(k));
		InsertRow(k, v);
	}

	return *this;
}

#endif // TMATRIX_HEADER