matrix.h

矩阵运算的模板类

template <typename T> Matrix<T> mtl_mpower(const Matrix<T>& A, double x){
	assert(A.m==A.n);
	if(x==0){
		return Matrix<T>::Eye(A.m,A.n);
	}
	else if(x>0){
		if(Abs(x-(int)x)<real(mtl_numeric_limits<T>::epsilon())){
			return mtl_mpower(A,(int)x);
		}
		Matrix<T> EVec, EVal;
		int isOK=Eig(A,EVec,EVal);
		Matrix<T> C(A.m,A.n); C=T();
		for(int i=1;i<=A.m;i++){
			if(A.IsReal() && real(EVal[i][i])<0){
				cerr << "ERROR: The A matrix has negative eigenvalues.\n";
				cerr << "       Use a complex matrix.\n";
				return Matrix<T>();
			}
			C[i][i]=pow(EVal[i][i],x);
		}
		C=EVec*C*Inv(EVec);
		return C;
	}
	else{
		return Powm(Inv(A),-x);
	}
}


template <typename T> Matrix<T> mtl_mpower(const Matrix<T>& A, int x){
	assert(A.m==A.n);
	if(x==0){
		return Matrix<T>::Eye(A.m,A.n);
	}
	else if(x>0){
		Matrix<T> C(A.m,A.n);
		C.SetIdentity();
		for(int i=1;i<=x;i++)
			C=C*A;
		return C;
	}
	else{
		int retCode;
		Matrix<T> C=Inv(A,&retCode);
		if(retCode==0) {
			cerr << "ERROR: Singular matrix.\n";
			assert(0);
			Matrix<T> Inf(A.m,A.n,mtl_numeric_limits<T>::max());
			return Inf;
		}
		return Powm(C,-x);
	}
}


template <typename T> ostream& mtl_ostream(ostream& s, const Matrix<T>& A) {
	for(int i=1; i<=A.m; i++) {
		s << "        ";
		for(int j=1; j<=A.n; j++) {
			s << Num2str(A[i][j]) << " ";
		}
		s << endl;
	}
	return s;
}


template <typename T> istream& mtl_istream(istream& s, Matrix<T>& A) {
	for (int i=1; i<=A.m; i++){
		for (int j=1; j<=A.n; j++)
			s >> A[i][j];
	}
	return s;
}



/////////////////////////////
// member functions
/////////////////////////////
template <typename T> int Matrix<T>::Resize(int new_m, int new_n) {
	if(m==new_m && n==new_n) return 1;
	// must preserve data
	int i;
	for(i=1;i<=local_min(m,new_m);i++) { mat[i].Resize(new_n); mat[i].SetType(ROW_VECTOR); }
	mat.Resize(new_m);
	for(i=m+1;i<=new_m;i++) { mat[i].Resize(new_n); mat[i].SetType(ROW_VECTOR); }
	m=new_m; n=new_n;
	return 1;
}


template <typename T> Vector<T> Matrix<T>::Sum_(int idx) const{
	assert(idx==1 || idx==2);
	Vector<T> C;
	int i,j;
	if(idx==1){
		C.SetType(ROW_VECTOR);
		C.Resize(n);
		for(j=1;j<=n;j++){
			T sum=0; for(i=1;i<=m;i++) sum+=mat[i][j]; C[j]=sum;
		}
	}
	else if(idx==2){
		C.SetType(COL_VECTOR);
		C.Resize(m);
		for(i=1;i<=m;i++){
			T sum=0; for(j=1;j<=n;j++) sum+=mat[i][j]; C[i]=sum;
		}
	}
	return C;
}


template <typename T> Vector<T> Matrix<T>::OutputAsRow(){
	Vector<T> C(m*n,ROW_VECTOR);
	int k=1;
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++)
		C[k++]=mat[i][j];
	return C;
}


template <typename T> Vector<T> Matrix<T>::OutputAsCol(){ 
	Vector<T> C(m*n,COL_VECTOR);
	int k=1;
	for(int j=1;j<=n;j++) for(int i=1;i<=m;i++) 
		C[k++]=mat[i][j];
	return C;
}


template <typename T> Vector<T> Matrix<T>::Row(int rowX) const {
	assert(rowX>=1 && rowX<=m);
	return mat[rowX];
}


template <typename T> Matrix<T> Matrix<T>::Rows(int m1,int m2) const {
	assert(m1>=1 && m1<=m && m2>=1 && m2<=m && m2>=m1);
	Matrix<T> C(m2-m1+1,n);
	for(int i=1;i<=m2-m1+1;i++) for(int j=1;j<=n;j++) 
		C[i][j]=mat[i+m1-1][j];
	return C;
}


template <typename T> Vector<T> Matrix<T>::Col(int colX) const {
	assert(colX>=1 && colX<=n);
	Vector<T> C(m,COL_VECTOR);
	for(int i=1;i<=m;i++) C[i]=mat[i][colX];
	return C;
}


template <typename T> Matrix<T> Matrix<T>::Cols(int n1,int n2) const {
	assert(n1>=1 && n1<=n && n2>=1 && n2<=n && n2>=n1);
	Matrix<T> C(m,n2-n1+1);
	for(int j=1;j<=n2-n1+1;j++) for(int i=1;i<=m;i++)
		C[i][j]=mat[i][j+n1-1];
	return C;
}


template <typename T> int Matrix<T>::SetRow(int rowX,const Vector<T>& A){
	assert(n==A.Size());
	mat(rowX)=A;
	return 1;
}


template <typename T> int Matrix<T>::SetCol(int colX,const Vector<T>& A){
	assert(m==A.Size());
	for(int i=1;i<=m;i++) mat[i][colX]=A[i];
	return 1;
}


template <typename T> 	int Matrix<T>::Print(ostream& s, const char *name) {
	if(name != 0 && name[0] != 0) s << name << " ";
	s << 2 << " " << m << " " << n << endl;
	s << (*this) << endl;
	return 1;
}


template <typename T> 	int Matrix<T>::Print(ostream& s, string& name) {
	return Print(s,name.c_str());
}


template <typename T> 	int Matrix<T>::Read(istream& s, const char *name) {
	string dummy;
	s >> dummy;
	if(name != 0 && name[0] != 0){
		assert(name==dummy);
	}
	int itmp, mm, nn;
	s >> itmp >> mm >> nn;
	assert(itmp==2);
	Resize(mm,nn);
	s >> (*this);
	return 1;
}


template <typename T> 	int Matrix<T>::Read(istream& s, string& name) {
	return Read(s,name.c_str());
}


template <typename T> int Matrix<T>::SetIdentity() {
	(*this)=T();
	for (int i=1; i<=local_min(m,n); i++) mat[i][i] = 1;
	return 1;
}


template <typename T> Matrix<T> Matrix<T>::SubMatrix(int m1,int m2, int n1, int n2) const{
	Matrix<T> C(m2-m1+1,n2-n1+1);
	for(int i=1;i<=m2-m1+1;i++){
		for(int j=1;j<=n2-n1+1;j++){
			C.mat[i][j]=mat[m1-1+i][n1-1+j];
		}
	}
	return C;
}


template <typename T> int Matrix<T>::NormCols(int j)
{
	if(j==0) for(int k=1;k<=n;k++) NormCols(k);	
	else{
		int i;
		double sqSum=0;
		for(i=1;i<=m;i++) sqSum+=real(mat[i][j]*conj(mat[i][j]));
		if(sqSum>0){
			double normWj=sqrt(sqSum);
			for(i=1;i<=m;i++) mat[i][j]/=normWj;
		}
	}
	return 1;
}


template <typename T> int Matrix<T>::NormRows(int i)
{
	if(i==0) for(int k=1;k<=n;k++) NormRows(k);	
	else{
		int j;
		double sqSum=0;
		for(j=1;j<=n;j++) sqSum+=real(mat[i][j]*conj(mat[i][j]));
		if(sqSum>0){
			double normWi=sqrt(sqSum);
			for(j=1;j<=n;j++) mat[i][j]/=normWi;
		}
	}
}


template <typename T> int Matrix<T>::Save(const string& fname){
	FILE* fp=fopen(fname.c_str(),"wb");
	if(!fp) return 0;
	fprintf(fp,"FMAT %d %d\n",m,n);
	for(int i=1;i<=m;i++){
		for(int j=1;j<=n;j++)
			fprintf(fp,"%f ",mat[i][j]);
		fprintf(fp,"\n");
	}
	fclose(fp);
	return m*n;
}


template <typename T> int Matrix<T>::SaveMatlab(const string& fname){
	FILE* fp=fopen(fname.c_str(),"wb");
	if(!fp) return 0;
	for(int i=1;i<=m;i++){
		for(int j=1;j<=n;j++)
			fprintf(fp,"%f ",mat[i][j]);
		fprintf(fp,"\n");
	}
	fclose(fp);
	return m*n;
}


template <typename T> int Matrix<T>::Input(const Vector<T>& V){
	int k=1;
	for(int i=1;i<=m;i++){
		for(int j=1;j<=n;j++)
			mat[i][j]=V[k++];
	}
	return m*n;
}


template <typename T> int Matrix<T>::Input(const Matrix<T>& A){
	int mm=local_min(m,A.m);
	int nn=local_min(n,A.n);
	for(int i=1;i<=mm;i++){
		for(int j=1;j<=nn;j++)
			mat[i][j]=A.mat[i][j];
	}
	return mm*nn;
}



/////////////////////////////
// basic friend functions
/////////////////////////////
template <typename T> inline Vector<int> Size(const Matrix<T>& A) {
	Vector<int> C(2,ROW_VECTOR);
	C(1)=A.m; C(2)=A.n;
	return C;
}


template <typename T> inline int Size(const Matrix<T>& A,int d) {
	assert(d==1 || d==2);
	return ((d==1) ? A.m : A.n);
}


template <typename T> inline int NumRows(const Matrix<T>& M) {
	return M.m;
}


template <typename T> inline int NumCols(const Matrix<T>& M) {
	return M.n;
}


template <typename T> inline double Scalar(const Matrix<T>& A) {
	assert(A.m==1 && A.n==1);
	return A[1][1];
}


// To remove memory allocation
template <typename T> T Multiply3(const Vector<T>& x, const Matrix<T>& B, const Vector<T>& y){
	assert(x.Type() == ROW_VECTOR);
	assert(y.Type() == COL_VECTOR);
	assert(x.Size()==B.m && B.n==y.Size());
	T sum=T();
	for(int j=1;j<=B.n;j++){
		T s=T();
		for(int i=1;i<=B.m;i++) s+=x[i]*B.mat[i][j];
		sum += s*y[j];
	}
	return sum;
}


// To remove memory allocation
template <typename T> T MultiplyXtAY(const Vector<T>& x, const Matrix<T>& B, const Vector<T>& y){
	assert(x.Type() == COL_VECTOR);
	assert(y.Type() == COL_VECTOR);
	assert(x.Size()==B.m && B.n==y.Size());
	T sum=T();
	for(int j=1;j<=B.n;j++){
		T s=T();
		for(int i=1;i<=B.m;i++) s+=x[i]*B.mat[i][j];
		sum += s*y[j];
	}
	return sum;
}


// Transpose: return transposed matrix of given matrix.
template <typename T> Matrix<T> Matrix<T>::t() const {
	int i,j;
	Matrix<T> B(n, m);
	for (i=1; i<=m; i++) for (j=1; j<=n; j++){
			B.mat[j][i]=(mat[i][j]);
	}
	return B;
}