写了一个模版矩阵类，学到了很多东西

    本来是写一个模版矩阵类，后来有添加了几个线性代数相关的函数，高斯消去法和LU分解；刚刚又加了一个QR分解函数，还想继续添加……

   

/************************************************************************/
/* 模版矩阵类，支持基本的数据类型，和stl的复数类，如果自定义的类型需	*/
/* 要支持运算符：+,-,*,/,+=,-=,*=,/=,=，常用的数学函数					*/
/************************************************************************/
/* By JJJ																*/
/* QQ：771588193														*/
/* 2013/11/23															*/
/************************************************************************/


#pragma once
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <assert.h>
using namespace std;


namespace myblas
{

/************************************************************************/
/*      简单的动态数组                                                  */
/************************************************************************/
template<typename Ty>
class basic_vector
{
public:
	basic_vector():v_data(NULL),v_len(0)	{}

	explicit
	basic_vector(int n):v_data(NULL),v_len(0)	{
		alloc(n);
	}

	basic_vector(const basic_vector &vec):v_data(NULL),v_len(0)	{
		assign(vec);
	}

	basic_vector(basic_vector &&_tmp):v_data(_tmp.v_data),v_len(_tmp.v_len)	{
		_tmp.v_data = NULL;		_tmp.v_len = 0;
	}

	~basic_vector()	{
		release();
	}

	//
	inline void alloc(int n)
	{
		v_data = new Ty[n]();
		if (NULL!=v_data)
			v_len = n;
	}

	inline void release()
	{
		if (NULL!=v_data)
		{
			delete[] v_data;
			v_data = NULL;
			v_len = 0;
		}
	}

	bool assign(const basic_vector &vec)
	{
		if (NULL!=vec.v_data)
		{
			alloc(vec.v_len);
			if (NULL!=v_data)
			{
				int i;
				for (i=0; i<v_len; i++)
					v_data[i] = vec.v_data[i];
				return true;
			}
		}
		return false;
	}

	inline int get_len()
	{
		return v_len;
	}

	//
	basic_vector& operator = (const basic_vector &vec)
	{
		assign(vec);
	}

	inline Ty &operator () (int i)
	{
#ifdef _CHECK_BOARDER_
		assert(0<=i);
		assert(i<v_len);
#endif
		return v_data[i];
	}

	//
protected:
	Ty		*v_data;
	int		v_len;
};

typedef basic_vector<int> array_int;
ostream& operator << (ostream &out,array_int &arr)
{
	int i,j;
	int n = arr.get_len();
	out<<endl;
	for (i=0; i<n; i++)
	{
		for (j=0; j<n; j++)
			if (arr(i)==j)
				out<<1<<'\t';
			else
				out<<0<<'\t';
		out<<endl;
	}
	return out;
}

//

/************************************************************************/
/*      矩阵模版类                                                      */
/************************************************************************/

typedef enum matrix_type	{
	NONE =0,	IDENTIFY,	ONEVAL,	UPPER_TRIANGULAR,	LOWER_TRIANGULAR
}Mat_Ty;

template<typename T>
class basic_matrix
{
public:
	//
	//construct and destruct
	basic_matrix():m_data(NULL),m_row(0),m_col(0)	{}

	explicit
	basic_matrix(int n,Mat_Ty mty =NONE,T val =T(1)):m_data(NULL),m_row(0),m_col(0)	{
		if (allocate(n,n))	set_type(mty,val);
	}

	explicit
	basic_matrix(int r,int c,Mat_Ty mty =NONE,T val =T(1)):m_data(NULL),m_row(0),m_col(0)	{
		if (allocate(r,c))	set_type(mty,val);
	}

	basic_matrix(const basic_matrix &src):m_data(NULL),m_row(0),m_col(0)	{
		assign(src);
	}

	explicit
	basic_matrix(const basic_matrix &src, int r1,int c1,int r2,int c2):m_data(NULL),m_row(0),m_col(0)	{
		assign(src,r1,c1,r2,c2);
	}

	basic_matrix( basic_matrix &&_right):m_data(_right.m_data),m_row(_right.m_row),m_col(_right.m_col)	{
		_right.m_data = NULL;	_right.m_row = 0;	_right.m_col = 0;
	}

	~basic_matrix()	{
		release();
	}


	//
	//manage memory and data

	inline bool allocate(int r, int c)
	{
		if (r>0 && c>0)
		{
			T *tmp = new T[r * c]();
			if (NULL!=tmp)
			{
				m_data = tmp;
				m_row = r;
				m_col = c;
				return true;
			}
		}
		return false;
	}

	inline void release()
	{
		if (NULL!=m_data)
		{
			delete[] m_data;
			m_data = NULL;
			m_row = m_col = 0;
		}
	}

	inline void resize(int r,int c,bool Keep_oldvalue =false)
	{
		if (r<=0 || c<=0)
			return;

		if (r*c>m_row*m_col)
		{
			T *tmp = new T[r*c]();
			if (NULL!=tmp)
			{
				if (Keep_oldvalue)
					for (int i=0; i<m_row*m_col; i++)
						tmp[i] = m_data[i];

				delete[] m_data;
				m_data = tmp;
				m_row = r;
				m_col = c;
			}
		}
		else
		{
			m_row = r;
			m_col = c;
		}
	}

	bool assign(const basic_matrix &other)
	{
		if (NULL!=other.m_data)
		{
			if (NULL==m_data && (m_row*m_col)<(other.m_row*other.m_col))
			{
				release();
				allocate(other.m_row,other.m_col);
			}
			else
			{
				m_row = other.m_row;
				m_col = other.m_col;
			}
			if (NULL!=m_data)
			{
				int i;
				for (i=0; i<m_row*m_col; ++i)
					m_data[i] = other.m_data[i];
				return true;
			}
		}
		return false;
	}

	bool assign(const basic_matrix &src,int r1,int c1,int r2,int c2)
	{
		bool test = ( NULL!=src.m_data && r1<r2 && c1<c2 &&
			r1>=0 && c1>=0 && r2<src.m_row && c2<src.m_col );
		if (test)
		{
			int nr = r2 - r1 +1;
			int nc = c2 - c1 +1;
			if (NULL==m_data && (m_row*m_col)<(nr*nc))
			{
				release();
				allocate(nr,nc);
			}
			else
			{
				m_row = nr;
				m_col = nc;
			}
			if (NULL!=m_data)
			{
				int s_col = src.m_col;
				const T *s_data = src.m_data;
				int i(0),j(0),k(0);
				for (i=r1-1; i<r2; ++i)
					for (j=c1-1; j<c2; ++j)
						m_data[k++] = s_data[i*s_col+j];
				return true;
			}
		}
		return false;
	}

	const T& max_abs()
	{
		int i,k;
		for (k=0,i=1; i<m_row*m_col; i++)
			if ( abs(m_data[k]) < abs(m_data[i]) )
				k = i;

		return m_data[k];
	}

	const T& min_abs()
	{
		int i,k;
		for (k=0,i=1; i<m_row*m_col; i++)
			if ( abs(m_data[k]) > abs(m_data[i]) )
				k = i;

		return m_data[k];
	}

	inline static void swap(basic_matrix &a, basic_matrix &b)
	{
		int t;
		t = a.m_row;	a.m_row = b.m_row; b.m_row = t;
		t = a.m_col;	a.m_col = b.m_col; b.m_col = t;
		T *td;
		td = a.m_data;	a.m_data = b.m_data;	b.m_data = td;
	}

	inline void swap_row(int r1,int r2)
	{
#ifdef _CHECK_BOARDER_
		assert(0<=r1 && r1<m_row);
		assert(0<=r2 && r2<m_row);
#endif
		if ( r1!=r2 )
		{
			int i;
			T tmp;
			for (i=0; i<m_col; i++)
			{
				tmp = m_data[r1*m_col+i];
				m_data[r1*m_col+i] = m_data[r2*m_col+i];
				m_data[r2*m_col+i] = tmp;
			}
		}
	}

	inline void swap_col(int c1,int c2)
	{
#ifdef _CHECK_BOARDER_
		assert(0<=c1 && c1<m_col);
		assert(0<=c2 && c2<m_col);
#endif
		if ( c1!=c2 )
		{
			int i;
			T tmp;
			for (i=0; i<m_row; i++)
			{
				tmp = m_data[i*m_col+c1];
				m_data[i*m_col+c1] = m_data[i*m_col+c2];
				m_data[i*m_col+c2] = tmp;
			}
		}
	}

	inline void set_one_val(T val)
	{
		int i;
		for (i=0; i<(m_row * m_col); ++i)
		{
			m_data[i] = val;
		}
	}

	inline void set_ident(T val)
	{
		int n = min(m_row,m_col);
		set_one_val(T(0));
		int i;
		for (i=0; i<n; i++)
		{
			m_data[i*m_col+i] = val;
		}
	}

	inline void set_type(Mat_Ty mty,T val =T(1))
	{
		switch (mty)
		{
		case NONE:	break;
		case ONEVAL:	set_one_val(val); break;
		case IDENTIFY:	set_ident(val); break;
		default:	break;
		}
	}
	inline int get_row()
	{
		return m_row;
	}

	inline int get_col()
	{
		return m_col;
	}

	inline int get_size()
	{
		return (m_row * m_col);
	}

	inline bool is_empty()
	{
		return (NULL==m_data);
	}

	//
	// operators
	inline T& operator () (int i, int j =0)
	{
#ifdef _CHECK_BOARDER_
		// just for debug
		assert(NULL!=m_data);
		assert(0<=i);
		assert(i<m_row);
		assert(0<=j);
		assert(j<m_col);
#endif // _CHECK_BORDER_

		return m_data[i*m_col+j];
	}

	basic_matrix &operator = (const basic_matrix &other)
	{
		assign(other);
		return *this;
	}

	basic_matrix &operator = (basic_matrix &&_right)
	{
		swap(*this,_right);
		return *this;
	}

	basic_matrix operator + (const basic_matrix &other)
	{
		return add_another(*this,other,true);
	}

	basic_matrix operator - (const basic_matrix &other)
	{
		return add_another(*this,other,false);
	}

	basic_matrix &operator += (const basic_matrix &other)
	{
		add_self(*this,other,true);
		return (*this);
	}

	basic_matrix &operator -= (const basic_matrix &other)
	{
		add_self(*this,other,false);
		return (*this);
	}

	basic_matrix operator * (const basic_matrix &other)
	{
		return mat_multi(*this, other);
	}

	basic_matrix operator * (const T &&alpha)
	{
		if (NULL!=m_data)
		{
			basic_matrix dest(m_row,m_col);
			if (NULL!=dest.m_data)
			{
				for (int i=0; i<m_row*m_col; ++i)
				{
					dest.m_data[i] = m_data[i] * alpha;
				}
			}
			return move(dest);
		}
		else
			return basic_matrix();
	}

	friend basic_matrix operator * (const T &&alpha, const basic_matrix &src)
	{
		if (NULL!=src.m_data)
		{
			int r = src.m_row;
			int c = src.m_col;
			basic_matrix dest(r,c);
			if (NULL!=dest.m_data)
			{
				for (int i=0; i<r*c; ++i)
				{
					dest.m_data[i] = alpha * src.m_data[i];
				}
			}
			return move(dest);
		}
		else
			return basic_matrix();
	}

	basic_matrix operator ^ (int n)
	{
		if (m_row != m_col)
		{
			cout<<"Inputs must be a scalar and a square matrix."<<endl;
			return basic_matrix();
		}
		if (0==n)
		{
			return basic_matrix(m_row,m_col,IDENTIFY);
		} 
		else if (1==n)
		{
			return (*this);
		}
		else if (-1==n)
		{
			return inverse(*this);
		}
		else if (1<n)
		{
			basic_matrix prod = mat_multi(*this,*this);
			for (n-=2; n; n--)
				prod = mat_multi(prod,*this);
			return move(prod);
		} 
		else // if (-1>n)
		{
			basic_matrix tmp = inverse(*this);
			basic_matrix prod = mat_multi(tmp,tmp);
			for (n+=2; n; n++)
				prod = mat_multi(prod,tmp);
			return move(prod);
		}
	}

	basic_matrix operator ~()
	{
		return transpose(*this);
	}

	//
	//math compute function

	static basic_matrix transpose(const basic_matrix &a)
	{
		if (NULL!=a.m_data)
		{
			if (1==a.m_row && 1==a.m_col)
			{
				return basic_matrix(a);
			}
			else
			{
				basic_matrix b(a.m_col,a.m_row);
				T *a_data = a.m_data;
				T *b_data = b.m_data;
				int I_loop = a.m_row;
				int J_loop = a.m_col;
				int i,j;
				for (i=0; i<I_loop; ++i)
					for (j=0; j<J_loop; ++j)
						b_data[j*I_loop+i] = a_data[i*J_loop+j];

				return move(b);
			}
		}
		else
			return basic_matrix();
	}

	// return a +/- b
	static basic_matrix add_another(const basic_matrix &a,const basic_matrix &b,bool OP_ADD)
	{
		bool test = (
			NULL!= a.m_data && 
			NULL!=b.m_data && 
			a.m_row==b.m_row && 
			a.m_col==b.m_col);
		if (test)
		{
			basic_matrix c(a.m_row,a.m_col);
			if (NULL!=c.m_data)
			{
				const T *a_data = a.m_data;
				const T *b_data = b.m_data;
				T *c_data = c.m_data;
				int i;
				if (OP_ADD)
				{
					for (i=0; i<(a.m_row * a.m_col); ++i)
					{
						c_data[i] = a_data[i] + b_data[i];
					}
				}
				else
				{
					for (i=0; i<(a.m_row * a.m_col); ++i)
					{
						c_data[i] = a_data[i] - b_data[i];
					}
				}
			}
			return move(c);
		}
		else
		{
			cout<<"The matrices dimensions must agree."<<endl;
			return basic_matrix();
		}
	}

	// a = a +/- b
	static bool add_self(basic_matrix &a,const basic_matrix &b,bool OP_ADD)
	{
		bool test = (
			NULL!= a.m_data && 
			NULL!=b.m_data && 
			a.m_row==b.m_row && 
			a.m_col==b.m_col);
		if (test)
		{
			T *a_data = a.m_data;
			const T *b_data = b.m_data;
			int i;
			if (OP_ADD)
			{
				for (i=0; i<(a.m_row * a.m_col); ++i)
					a_data[i] += b_data[i];
			}
			else
			{
				for (i=0; i<(a.m_row * a.m_col); ++i)
					a_data[i] -= b_data[i];
			}
		}
		else
			cout<<"The matrices dimensions must agree."<<endl;
		return test;
	}

	//	compute c = a * b
	static basic_matrix mat_multi(const basic_matrix &a,const basic_matrix &b)
	{
		bool test = (NULL!= a.m_data && NULL!=b.m_data && a.m_col==b.m_row );
		if (test)
		{
			basic_matrix c(a.m_row,b.m_col);
			c.set_one_val(T(0));
			if (NULL!=c.m_data)
			{
				const T *A = a.m_data;
				const T *B = b.m_data;
				T *C = c.m_data;
				int I_loop = a.m_row;
				int J_loop = b.m_col;
				int K_loop = a.m_col;
				int i,j,k;

				for (i=0; i<I_loop; ++i)
					for (k=0; k<K_loop; ++k)
						for (j=0; j<J_loop; ++j)
							C[i*J_loop+j] += A[i*K_loop+k] * B[k*J_loop+j];

			}
			return move(c);
		}
		else
		{
			cout<<"Inner matrix dimensions must agree."<<endl;
			return basic_matrix();
		}
	}


	//
protected:
	T		*m_data;		//data pointer
	int		m_row,m_col;	//matrix size
};




//

// Extraction operator for matrix elements.
template<typename T>
istream& operator >> 
	(istream& _is, basic_matrix<T>& _mat)
{
	int r = _mat.get_row();
	int c = _mat.get_col();
	int i,j;

	for (i=0; i<r; i++)
	{
		cout<<"input rows "<<i+1<<": "<<endl;
		for (j=0; j<r; j++)
		{
			_is >> _mat(i,j);
		}
	}
	return _is;
}

// Insertion operator for matrix elements.
template<typename T,class _Elem,class _Tr>
basic_ostream<_Elem, _Tr>& operator <<
	(basic_ostream<_Elem, _Tr>& _os, basic_matrix<T>& _mat)
{
	if (_mat.is_empty())
		return (_os << endl << "No data in this matrix" << endl );

	// display all elements of basic_matrix<_Ty>
	basic_ostringstream<_Elem, _Tr, allocator<_Elem> > _Sstr;
	_Sstr.flags(_os.flags());
	_Sstr.imbue(_os.getloc());
	_Sstr.precision(_os.precision());
	int r = _mat.get_row();
	int c = _mat.get_col();
	int i,j;
	_Sstr << '[' << endl;
	for (i=0; i<r; i++)
	{
		for (j=0; j<c; j++)
		{
			_Sstr << fixed << _mat(i,j) << '\t';
		}_Sstr << endl;
	}_Sstr << ']' <<endl;
	basic_string<_Elem, _Tr, allocator<_Elem> > _Str = _Sstr.str();
	return (_os << _Str.c_str());
}

template<typename T,class _Elem,class _Tr> inline
	basic_ostream<_Elem, _Tr>& operator <<
	(basic_ostream<_Elem, _Tr>& _os, basic_matrix<T>&& _mat)
{
	return (_os << forward< basic_matrix<T>& >(_mat) );
}

//
// out/in binary files
template<typename T>
ifstream& operator >> (ifstream &_ifs, basic_matrix<T> &_mat)
{
	// need the data type (T) with simple struct.
	int size = _mat.get_row() * _mat.get_col() * sizeof(T);
	_ifs.read((char *)&_mat(0),size);
	return _ifs;
}

template<typename T>
ofstream& operator << (ofstream &_ofs, basic_matrix<T> &_mat)
{
	// need the data type (T) with simple struct.
	int size = _mat.get_row() * _mat.get_col() * sizeof(T);
	_ofs.write((char *)&_mat(0),size);
	return _ofs;
}


/************************************************************************/
/* functions                                                            */
/************************************************************************/
// convert value type
template<typename T1,typename T2>
bool convert(basic_matrix<T2> &dst,basic_matrix<T1> &src)
{
	bool test = ( src.is_empty() || dst.get_row()==src.get_row() || dst.get_col()==src.get_col() );
	if (test)
	{
		int i,j;
		for (i=0; i<src.get_row(); i++)
			for (j=0; j<src.get_col(); j++)
				dst(i,j) = T2(src(i,j));cout<<"convertint..."<<endl;
	}
	return test;
}

//
template<typename T>
basic_matrix<T> inverse(basic_matrix<T> &a)
{
	if (!a.is_empty() && a.get_row()==a.get_col())
	{
		basic_matrix<T> E(a.get_row());
		E.set_ident(T(1));
		return gauss_jordan(a,E);
	}
	else
	{
		cout<<"matrix must be square"<<endl;
		return basic_matrix<T>();
	}
}
template<typename T>
basic_matrix<T> inverse(basic_matrix<T> &&a)
{
	return inverse( forward< basic_matrix<T>& > (a));
}


// Gauss eliminate
template<typename T>
basic_matrix<T> gauss_elimi(basic_matrix<T> A,basic_matrix<T> b)
{
	bool test = ( !A.is_empty() && A.get_row()==A.get_col() && A.get_row()==b.get_row() );
	if (test)
	{
		int n = A.get_row();
		int col_b = b.get_col();
		int i,j,k;
		T C,tmp;

		// P is permutation matrix
		array_int P(n);
		for (i=0; i<n; i++)
			P(i) = i;

		for (i=0; i<n; ++i)
		{
			// choose the pivot in columns
			k = i;
			for (j=i+1; j<n; j++)
				if ( abs(A(P(j),i)) > abs(A(P(k),i)) )
					k = j;
			if ( abs(A(P(k),i))==abs(T(0)) )
				return basic_matrix<T>();
			if (k!=i)
				swap(P(i),P(k));

			for (j=i+1; j<n; j++)
			{
				C = A(P(j),i) / A(P(i),i);

				for (k=i; k<n; k++)
					A(P(j),k) -= C * A(P(i),k);

				for (k=0; k<col_b; k++)
					b(P(j),k) -= C * b(P(i),k);
			}
		}

		basic_matrix<T> X(b.get_row(),b.get_col(),ONEVAL,T(0));
		for (i=n-1; i>=0; --i)
		{
			C = T(1) / A(P(i),i);
			for (j=0; j<col_b; j++)
			{
				tmp = T(0);
				for (k=i+1; k<n; k++)
					tmp += A(P(i),k) * X(k,j);
				X(i,j) = C * ( b(P(i),j) - tmp );
			}
		}
		return move(X);
	}
	else
	{
		cout<<"matrix dimensions must agree."<<endl;
		return basic_matrix<T>();
	}
}

// Gauss-Jordan elimination
template<typename T>
basic_matrix<T> gauss_jordan(basic_matrix<T> A,basic_matrix<T> X)
{
	bool test = ( !A.is_empty() && A.get_row()==A.get_col() && A.get_row()==X.get_row() );
	if (test)
	{
		int r = A.get_row();
		int col_a = A.get_col();
		int col_x = X.get_col();
		array_int P(r);
		int i,j,k,index;
		T C;

		// 下三角变换
		for (i=0; i<r-1; ++i)
		{
			// choose the pivot in rows
			index = 0;
			for (k=1; k<col_a; k++)
			{
				if ( abs(A(i,k)) > abs(A(i,index)) )
					index = k;
			}
			P(i) = index;

			// then elimination
			for (j=i+1; j<r; ++j)
			{
				C = A(j,index) / A(i,index);

				for (k=0; k<col_a; ++k)
					A(j,k) -= C * A(i,k);

				for (k=0; k<col_x; ++k)
					X(j,k) -= C * X(i,k);
			}
		}
		index = 0;
		for (k=1; k<col_a; k++)
		{
			if ( abs(A(i,k)) > abs(A(i,index)) )
				index = k;
		}
		P(i) = index;
		
		// up elimination
		for (i=r-1; i>0; --i)
		{
			index = P(i);
			for (j=i-1; j>=0; --j)
			{
				C = A(j,index) / A(i,index);
				for (k=0; k<col_a; ++k)
					A(j,k) -= C * A(i,k);

				for (k=0; k<col_x; ++k)
					X(j,k) -= C * X(i,k);
			}
		}
		
		for (i=0; i<r; ++i)
		{
			C = T(1) / A(i,P(i));
			for (j=0; j<col_x; ++j)
				X(i,j) *= C;
		}
		//  reorder
		basic_matrix<T> Y(r,col_x);
		for (i=0; i<r; i++)
		{
			for (j=0; j<col_x; j++)
			{
				Y(P(i),j) = X(i,j);
			}
		}
		return move(Y);
	}
	else
	{
		cout<<"In function gauss_jordan, matrix dimensions must agree."<<endl;
		return basic_matrix<T>();
	}
}

// LU: Doolittle's method
template<typename T>
array_int doolittle(basic_matrix<T> &a, basic_matrix<T> &L, basic_matrix<T> &U)
{
	assert(!a.is_empty() && a.get_row()==a.get_col());

	int n = a.get_row();
	if (U.get_row()!=n || U.get_col()!=n)
		U.resize(n,n);
	if (L.get_row()!=n || L.get_col()!=n)
		L.resize(n,n);
	U.set_one_val(T(0));
	L.set_one_val(T(0));

	int r,i,j,k;
	T tmp;
	
	array_int P(n);	// P(n): permutation matrix
	for (i=0; i<n; i++)
		P(i) = i;
	basic_vector<T> S(n);	// S(n): maximum element in each row
	for (i=0; i<n; i++)
	{
		tmp = a(i,0);
		for (j=1; j<n; j++)
			if ( abs(tmp) < abs(a(i,j)) )
				tmp = a(i,j);
		S(i) = tmp;
	}

	for (r=0; r<n; r++)
	{
		// choose pivot
		k = r;
		for (j=r+1; j<n; j++)
			if ( (abs(a(P(k),r))/abs(S(k))) < (abs(a(P(j),r))/abs(S(j))) )
				k = j;
		if (k!=r)
			swap(P(k),P(r));

		// U is an upper triangular matrix
		for (j=r; j<n; j++)
		{
			tmp = T(0);
			for (k=0; k<r; k++)
				tmp += L(P(r),k) * U(k,j);

			U(r,j) = a(P(r),j) - tmp;
		}

		// L is a lower triangular matrix
		L(P(r),r) = T(1);
		for (i=r+1; i<n; i++)
		{
			tmp = T(0);
			for (k=0; k<r; k++)
				tmp += L(P(i),k) * U(k,r);

			L(P(i),r) = ( a(P(i),r) - tmp ) / U(r,r);
		}
	}

	/*
	int r,k,i,j;
	T s;
	for (r=0; r<n; r++)
	{
		//
		for (j=r; j<n; j++)
		{
			s = T(0);
			for (k=0; k<r; k++)
				s += L(r,k) * U(k,j);

			U(r,j) = a(r,j) - s;
		}
		//
		for (i=r+1; i<n; i++)
		{
			s = T(0);
			for (k=0; k<r; k++)
				s += L(i,k) * U(k,r);

			L(i,r) = (a(i,r) - s) / U(r,r);
		}
	}
	*/

	return move(P);
}
//

// LU: method from Doolittle, 
template<typename T>
array_int doolittle(basic_matrix<T> &a)
{
	assert(!a.is_empty());
	assert(a.get_row()==a.get_col());

	int n = a.get_row();
	array_int P(n);
	int r,k,i,j;
	T tmp;

	for (i=0; i<n; i++)
		P(i) = i;

	for (r=0; r<n; r++)
	{
		k = r;
		for (j=r+1; j<n; j++)
		{
			if ( abs(a(j,r)) < abs(a(k,r)))
				k = j;
		}
		if (k!=r)
			swap(P(k),P(r));
		a.swap_row(k,r);
		//
		for (j=r; j<n; j++)
		{
			tmp = T(0);
			for (k=0; k<r; k++)
				tmp += a(r,k) * a(k,j);

			a(r,j) = a(r,j) - tmp;
		}
		//
		for (i=r+1; i<n; i++)
		{
			tmp = T(0);
			for (k=0; k<r; k++)
				tmp += a(i,k) * a(k,r);

			a(i,r) = (a(i,r) - tmp) / a(r,r);
		}
	}

	return move(P);
}


// deter of matrix
template<typename T>
T deter(basic_matrix<T> a)
{
	if (!a.is_empty() && a.get_row() == a.get_col())
	{
		int n = a.get_row();
		array_int P(n);
		bool sig = false;
		int i,j,k,index;
		T C;
		for (i=0; i<n; i++)
			P(i) = i;
		//
		for (i=0; i<n-1; i++)
		{
			index = i;
			for (k=i+1; k<n; k++)
				if ( abs(a(i,index)) < abs(a(i,k)) )
					index = k;

			if (index!=i)
			{
				swap(P(i),P(index));
				sig = !sig;
			}

			for (j=i+1; j<n; j++)
			{
				C = a(P(j),i) / a(P(i),i);

				for (k=0; k<n; k++)
					a(P(j),k) -= C * a(P(i),k);
			}
		}

		C = sig?T(-1):T(1);
		for (int i=0; i<n; ++i)
		{
			C *= a(P(i),i);
		}

		return C;
	}
	else
	{
		cout<<"matrix must be square!"<<endl;
		return T(0);
	}
}

// solve equation
template<typename T>
basic_matrix<T> solve(basic_matrix<T> &A,basic_matrix<T> &b)
{
	assert(!A.is_empty());
	assert(A.get_row()==A.get_col());
	assert(A.get_row()==b.get_row());

	int n = A.get_row();
	int col_b = b.get_col();
	int i,j,k;
	T C,tmp;
	basic_matrix<T> L(n),U(n);
	basic_matrix<T> Y(n,b.get_col());
	basic_matrix<T> X(n,b.get_col());
	array_int P = doolittle(A,L,U);

	for (i=0; i<n; i++)
	{
		for (j=0; j<col_b; j++)
		{
			tmp = T(0);
			for (k=0; k<i; k++)
				tmp += L(P(i),k) * Y(k,j);

			Y(i,j) = b(P(i),j) - tmp ;
		}
	}
	
	for (i=n-1; i>=0; --i)
	{
		C = U(i,i);
		for (j=0; j<col_b; j++)
		{
			tmp = T(0);
			for (k=i+1; k<n; k++)
				tmp += U(i,k) * X(k,j);

			X(i,j) = (Y(i,j) - tmp) / C;
		}
	}

	return move(X);
}

// solve equation
template<typename T>
basic_matrix<T> solve2(basic_matrix<T> a,basic_matrix<T> b)
{
	int n = a.get_row();
	int col_b = b.get_col();
	int i,j,k;
	T C,tmp;
	basic_matrix<T> Y(n,col_b);

	array_int P = doolittle(a);

	for (i=0; i<n; i++)
		for (j=0; j<col_b; j++)
			Y(i,j) = b(P(i),j);

	for (i=0; i<n; i++)
	{
		for (j=0; j<col_b; j++)
		{
			tmp = T(0);
			for (k=0; k<i; k++)
				tmp += a(i,k) * b(k,j);

			b(i,j) = Y(i,j) - tmp ;
		}
	}

	for (i=n-1; i>=0; --i)
	{
		C = a(i,i);
		for (j=0; j<col_b; j++)
		{
			tmp = T(0);
			for (k=i+1; k<n; k++)
				tmp += a(i,k) * Y(k,j);

			Y(i,j) = (b(i,j) - tmp) / C;
		}
	}

	return move(Y);
}


//
// QR: the data type used operator <,>,<=,>= and function sqrt(..)
template<typename T>
bool QR(basic_matrix<T> A,basic_matrix<T> &Q,basic_matrix<T> &R)
{
	if (A.is_empty())
		return false;
	if (Q.is_empty())
		Q.resize(A.get_row(),A.get_row());
	Q.set_ident(T(1));

	int a_row = A.get_row();
	int a_col = A.get_col();
	int N = min(a_row, a_col);
	int i,j,k,t;
	T sum,dr,cr,hr;
	basic_vector<T> ur(a_row), pr(a_col), wr(a_row);

	for (t=0; t<N; t++)
	{
		sum	= T(0);
		for (k=t; k<a_row; k++)
			sum += A(k,t) * A(k,t);
		dr = sqrt(sum);
		cr = A(t,t)>T(0)?(-dr):dr;
		hr = cr *( cr - A(t,t) );

		// ur
		for (i=0; i<t; i++)
			ur(i) = T(0);
		ur(t) = A(t,t) - cr;
		for (i=t+1; i<a_row; i++)
			ur(i) = A(i,t);

		// wr
		for (i=0; i<a_row; i++)
		{
			sum = T(0);
			for (j=0; j<a_row; j++)
				sum += Q(i,j) * ur(j);
			wr(i) = sum;
		}

		// Q
		for(i=0; i<a_row; i++)
			for(j=0; j<a_row; j++)
				Q(i,j) -= wr(i) * ur(j) / hr;

		// pr
		for(i=0; i<a_col; i++)
		{
			sum=0;
			for(j=0; j<a_row; j++)
				sum += A(j,i) * ur(j);
			pr(i) = sum / hr;
		}

		// R
		for(i=0;i<a_row;i++)
			for(j=0;j<a_col;j++)
				A(i,j) -= ur(i) * pr(j);

	}
	basic_matrix<T>::swap(A,R);
	return true;
}

}
//

 

    我测试用的cpp文件，比较乱……

#include <stdlib.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include <complex>
#include <cmath>
using namespace std;
#ifndef _CHECK_BOARDER_
#define _CHECK_BOARDER_
#endif // !_CHECK_BORDER_
#include "basic_matrix.h"
using namespace myblas;
typedef double realnum;
typedef realnum data_t;
typedef basic_matrix<data_t> matrix;

//
void rand(matrix &a)
{
	srand(clock());
	int i,j;
	data_t C = data_t(1<<9);
	for (i=0; i<a.get_row(); i++)
		for (j=0; j<a.get_col(); j++)
			a(i,j) = data_t(rand()&((1<<10)-1))/C;
}


realnum test(matrix &a, matrix &b)
{
	if (a.get_row()!=b.get_row() && a.get_col()!=b.get_col())
		return -1;
	realnum esp = 0.0;
	for (int i = 0; i < a.get_row(); i++)
	{
		for (int j = 0; j < a.get_col(); j++)
		{
			esp = max(esp,abs(a(i,j)-b(i,j)));
		}
	}
	return esp;
}

realnum test(matrix &&a, matrix &&b){   return test(forward<matrix&>(a),forward<matrix&>(b));    }
realnum test(matrix &a, matrix &&b){   return test(forward<matrix&>(a),forward<matrix&>(b));    }
realnum test(matrix &&a, matrix &b){   return test(forward<matrix&>(a),forward<matrix&>(b));    }

matrix rand(int n,int m)
{
	matrix a(n,m);
	rand(a);
	return move(a);
}

matrix rand(int n)
{
	matrix a(n);
	rand(a);
	return move(a);
}

int main(int argc,char *argv[])
{
	int m,n;
	if (2==argc)
	{
		argv++;
		m=atoi(argv[1]);
		n = m;
	}
	else if (2<argc)
	{
		m=atoi(argv[1]);
		n=atoi(argv[2]);
	}
	else
		cin>>m>>n;
	{
		matrix a=rand(m,n),q,r;

/*
		ifstream fin;
		fin.open("sample",ios::binary | ios::in);
		fin>>a;
		fin>>b;
		fin>>x0;
		fin.close();
*/

/*
		x = solve(a,b);
		cout<<"solve 1: "<<abs((a*x-b).max_abs())<<endl;
		x = solve2(a,b);
		cout<<"solve 2: "<<abs((a*x-b).max_abs())<<endl;
		x = gauss_elimi(a,b);
		cout<<"solve 3: "<<abs((a*x-b).max_abs())<<endl;
		x = gauss_jordan(a,b);
		cout<<"solve 4: "<<abs((a*x-b).max_abs())<<endl;
*/

		cout.precision(10);
/*		
		cout<<"a:"<<a<<endl;
		cout<<inverse(a*a*a*a)<<endl;
		cout<<(a^-4)<<endl;
*/
		if (QR(a,q,r))
		{
			int t = 16;
			
			cout<<"q:"<<q<<endl;
			cout<<"r:"<<r<<endl;
			cout<<"q*r:"<<q*r<<endl;
			cout<<q*r<<endl;
			for (int i=0; i<t; i++)
				QR(r*q,q,r);
			cout<<"r*q: "<<r*q<<endl;
		}


/*
		ofstream fout;
		fout.open("data",ios::binary | ios::out);
		fout<<a;
		fout<<b;
		fout<<x;
		fout.close();
*/

	}
	
//	system("pause");
	return 0;
}