matrix.h

该文件包含《线性代数方程组的求解》、《矩阵运算》、《矩阵特征值与特征向量的计算》、《数值积分》等算法与C程序

/*********************************************************************
 Matrix.h: interface&implementation for the CMatrix class.
 Write by hujinshan@2003.6.26 Airforce Engineering University
*********************************************************************/

#if !defined(HJS_MATRIX_H)
#define HJS_MATRIX_H

#include <valarray>
using namespace std;

template<class T>
class Slice_iter{
	valarray<T>* v;
	slice s;
	size_t curr;//index of current element

	T& ref(size_t i)const{return (*v)[s.start()+i*s.stride()];}
public:
	Slice_iter(valarray<T>* vv,slice ss):v(vv),s(ss),curr(0){}
	Slice_iter end()
	{
		Slice_iter t=*this;
		t.curr=s.size();//index of last-plus-one element
		return t;
	}
	Slice_iter& operator++(){curr++;return *this;}
	Slice_iter operator++(int){Slice_iter t=*this;curr++;return t;}
	size_t size(){return v->size();}
	T& operator[](size_t i){return ref(curr=i);}//C-style subscript
	T& operator()(size_t i){return ref(curr=i);}//Fortran-style subscript
	T& operator*(){return ref(curr);}//current element
	//...
};

template<class T>
class Cslice_iter{
	valarray<T>* v;
	slice s;
	size_t curr;//index of current element

	T& ref(size_t i)const{return (*v)[s.start()+i*s.stride()];}
public:
	Cslice_iter(valarray<T>* vv,slice ss):v(vv),s(ss),curr(0){}
	Cslice_iter end()
	{
		Cslice_iter t=*this;
		t.curr=s.size();//index of last-plus-one element
		return t;
	}
	Cslice_iter& operator++(){curr++;return *this;}
	Cslice_iter operator++(int){Cslice_iter t=*this;curr++;return t;}
	size_t size(){return v->size();}
	T& operator[](size_t i){return ref(curr=i);}//C-style subscript
	T& operator()(size_t i){return ref(curr=i);}//Fortran-style subscript
	T& operator*(){return ref(curr);}//current element
	//...
};

class CMatrix  
{
	valarray<double>* v;
	size_t d1, d2;
public:
	CMatrix(size_t x,size_t y);
	CMatrix(const CMatrix&);
	CMatrix& operator=(const CMatrix&);

	size_t size()const{return d1*d2;}
	size_t dim1()const{return d1;}
	size_t dim2()const{return d2;}

	Slice_iter<double> row(size_t i);
	Cslice_iter<double> row(size_t i)const;
	Slice_iter<double> column(size_t i);	
	Cslice_iter<double> column(size_t i)const;

	double& operator()(size_t x,size_t y);//Fortran-style subscripts
	double operator()(size_t x,size_t y)const;

	Slice_iter<double> operator()(size_t i){return row(i);}//C-style subscripts
	Cslice_iter<double> operator()(size_t i)const{return row(i);}

	Slice_iter<double> operator[](size_t i){return row(i);}//C-style subscripts
	Cslice_iter<double> operator[](size_t i)const{return row(i);}

	CMatrix& operator*=(double);
	valarray<double>& array(){return *v;}
	virtual ~CMatrix();

	void ZeroMatrix();
};

CMatrix::CMatrix(size_t x,size_t y)
{
	d1=x;d2=y;
	v=new valarray<double>(d1*d2);
}
CMatrix::CMatrix(const CMatrix& mat)
{
	(*this)=mat;
}

CMatrix& CMatrix::operator=(const CMatrix& mat)
{
	d1=mat.dim1();
	d2=mat.dim2();
	v=new valarray<double>(d1*d2);
	for(int i=0;i<d1;i++)
		for(int j=0;j<d2;j++)
			(*this) [i][j]=mat[i][j];
	return *this;
}

CMatrix::~CMatrix()
{
	delete v;
}
double& CMatrix::operator() (size_t x,size_t y)
{
	return row(x)[y];
}
static double mul(Slice_iter<double>& v1,const valarray<double>& v2)
{
	double res=0;
	for(int i=0;i<v1.size();i++)
		res+=v1[i]*v2[i];
	return res;
}
static valarray<double> operator*(CMatrix& m,const valarray<double>& v)
{
	valarray<double> res(m.dim1());
	for(int i=0;i<m.dim1();i++)
		res[i]=mul(m.row(i),v);
	return res;
}
CMatrix& CMatrix::operator*=(double d)
{
	(*v) *=d;
	return *this;
}
inline Slice_iter<double> CMatrix::row(size_t i)
{
	return Slice_iter<double>(v,slice(i,d2,d1));
}
inline Cslice_iter<double> CMatrix::row(size_t i)const
{
	return Cslice_iter<double>(v,slice(i,d2,d1));
}
inline Slice_iter<double> CMatrix::column(size_t i)
{
	return Slice_iter<double>(v,slice(i*d2,d2,1));
}
inline Cslice_iter<double> CMatrix::column(size_t i)const
{
	return Cslice_iter<double>(v,slice(i*d2,d2,1));
}
void CMatrix::ZeroMatrix()
{
	for(int i=0;i<size();i++)
		(*v)[i]=0;
}
#endif //!defined(HJS_MATRIX_H)
//end of file