matrix.h

矩阵运算的模板类

// constructor and destructor
/////////////////////////////
template <typename T> Matrix<T>::Matrix(int m_, int n_, const T& v_) : m(m_), n(n_) {
	int i;
	mat.Resize(m);
	for(i=1;i<=m;i++) { mat[i].SetType(ROW_VECTOR); mat[i].Resize(n); }
	for(i=1;i<=m;i++) for(int j=1;j<=n;j++)  mat[i][j]=v_;
}


template <typename T> Matrix<T>::Matrix(int m_, int n_, const char* v_) : m(m_), n(n_) {
	int i,k;
	mat.Resize(m);
	for(i=1;i<=m;i++) { mat[i].SetType(ROW_VECTOR); mat[i].Resize(n); }
	Vector<string> num=mtl_split(v_);
	for(i=1,k=1;i<=m;i++) for(int j=1;j<=n;j++) {
		if(k<=num.Size()) { mat[i][j]=T(atof(num[k].c_str())); k++; }
		else mat[i][j]=T();
	}
}


template <typename T> Matrix<T>::Matrix(const Vector<T>& v_, int m_, int n_) : m(m_), n(n_) {
	int i,k;
	mat.Resize(m);
	for(i=1;i<=m;i++) { mat[i].SetType(ROW_VECTOR); mat[i].Resize(n); }
	for(i=1,k=1;i<=m;i++) for(int j=1;j<=n;j++) {
		if(k<=v_.Size()) { mat[i][j]=v_[k]; k++; }
		else mat[i][j]=T();
	}
}


//template <typename T> inline Matrix<T>::Matrix(const Matrix<T>& A) : mat(A.mat),m(A.m),n(A.n) {
//}


template <typename T> inline Matrix<T>::~Matrix() {
}


/////////////////////////////
// operator overloading
/////////////////////////////
//template <typename T> inline Matrix<T>& Matrix<T>::operator=(const Matrix<T>& A) {
//	if(this == &A) return *this;
//	if ( A.n != n || A.m != m) Resize(A.m, A.n);
//	mat = A.mat;
//	return *this;
//}


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


#ifdef ALLOW_ACCESS_BY_BRACKET
template <typename T> inline Vector<T>& Matrix<T>::operator[](int i) { 
	return mat[i]; 
}


template <typename T> inline const Vector<T>& Matrix<T>::operator[](int i) const { 
	return mat[i]; 
}
#endif


template <typename T> inline Vector<T>& Matrix<T>::operator()(int i) {
	return mat(i); 
}


template <typename T> inline const Vector<T>& Matrix<T>::operator()(int i) const {
	return mat(i); 
}


template <typename T> inline T& Matrix<T>::operator()(int i, int j) {
	return (mat(i))(j); 
}


template <typename T> inline const T& Matrix<T>::operator()(int i, int j) const {
	return (mat(i))(j); 
}


template <typename T> inline Matrix<T> Matrix<T>::operator+() {
	return (*this);
}


template <typename T> Matrix<T> Matrix<T>::operator-() {
	Matrix<T> C=(*this);
	C*=(-1);
	return C;
}


template <typename T> Matrix<T>& Matrix<T>::operator+=(const Matrix<T>& B) {
	assert(m == B.m && n == B.n);
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]+=B.mat[i][j];
	return (*this);
}


template <typename T> Matrix<T>& Matrix<T>::operator-=(const Matrix<T>& B) {
	assert(m == B.m && n == B.n);
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]-=B.mat[i][j];
	return (*this);
}


template <typename T> Matrix<T>& Matrix<T>::operator+=(const double b) {
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]+=b;
	return (*this);
}


template <typename T> Matrix<T>& Matrix<T>::operator-=(const double b) {
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]-=b;
	return (*this);
}


template <typename T> Matrix<T>& Matrix<T>::operator*=(const double b) {
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]*=b;
	return (*this);
}


template <typename T> Matrix<T>& Matrix<T>::operator/=(const double b) {
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]/=b;
	return (*this);
}


#ifndef DISABLE_COMPLEX
template <typename T> Matrix<T>& Matrix<T>::operator+=(const complex<T> b) {
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]+=b;
	return (*this);
}


template <typename T> Matrix<T>& Matrix<T>::operator-=(const complex<T> b) {
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]-=b;
	return (*this);
}


template <typename T> Matrix<T>& Matrix<T>::operator*=(const complex<T> b) {
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]*=b;
	return (*this);
}


template <typename T> Matrix<T>& Matrix<T>::operator/=(const complex<T> b) {
	for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) mat[i][j]/=b;
	return (*this);
}
#endif


/////////////////////////////
// friend operator overloading
/////////////////////////////
template <typename T> inline Matrix<T> mtl_plus(const Matrix<T>& A, const Matrix<T>& B) {
	Matrix<T> C=A; C+=B; return C;
}


template <typename T> inline Matrix<T> mtl_minus(const Matrix<T>& A, const Matrix<T>& B) {
	Matrix<T> C=A; C-=B; return C;
}


template <typename T> inline Matrix<T> mtl_plus(const Matrix<T>& A, const double b) {
	Matrix<T> C=A; C+=b; return C;
}


template <typename T> inline Matrix<T> mtl_minus(const Matrix<T>& A, const double b) {
	Matrix<T> C=A; C-=b; return C;
}


template <typename T> inline Matrix<T> mtl_times(const Matrix<T>& A, const double b) {
	Matrix<T> C=A; C*=b; return C;
}


template <typename T> inline Matrix<T> mtl_divide(const Matrix<T>& A, const double b) {
	Matrix<T> C=A; C/=b; return C;
}


template <typename T> inline Matrix<T> mtl_plus(const double b, const Matrix<T>& A) {
	Matrix<T> C=A; C+=b; return C;
}


template <typename T> inline Matrix<T> mtl_minus(const double b, const Matrix<T>& A) {
	Matrix<T> C=-A; C+=b; return C;
}


template <typename T> inline Matrix<T> mtl_times(const double b, const Matrix<T>& A) {
	Matrix<T> C=A; C*=b; return C;
}


template <typename T> Matrix<T> mtl_divide(const double b, const Matrix<T>& A) {
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=C.m;i++) for(int j=1;C.n;j++)
		C.mat[i][j]=b/A.mat[i][j];
	return C;
}


// Multiply Matrix * Matrix
template <typename T> Matrix<T> mtl_mtimes(const Matrix<T>& A, const Matrix<T>& B) {
	assert(A.n == B.m);
	Matrix<T> C(A.m, B.n);
	for(int i=1;i<=A.m;i++) {
		for(int j=1;j<=B.n;j++){
			T sum=T();
			for(int k=1;k<=A.n;k++)
				sum+=A.mat[i][k]*B.mat[k][j];
			C.mat[i][j]=sum;
		}
	}
	return C;
}


#ifndef DISABLE_COMPLEX
template <typename T> inline Matrix<T> mtl_plus(const Matrix<T>& A, const complex<T> b) {
	Matrix<T> C=A; C+=b; return C;
}


template <typename T> inline Matrix<T> mtl_minus(const Matrix<T>& A, const complex<T> b) {
	Matrix<T> C=A; C-=b; return C;
}


template <typename T> inline Matrix<T> mtl_times(const Matrix<T>& A, const complex<T> b) {
	Matrix<T> C=A; C*=b; return C;
}


template <typename T> inline Matrix<T> mtl_divide(const Matrix<T>& A, const complex<T> b) {
	Matrix<T> C=A; C/=b; return C;
}


template <typename T> inline Matrix<T> mtl_plus(const complex<T> b, const Matrix<T>& A) {
	Matrix<T> C=A; C+=b; return C;
}


template <typename T> inline Matrix<T> mtl_minus(const complex<T> b, const Matrix<T>& A) {
	Matrix<T> C=-A; C+=b; return C;
}


template <typename T> inline Matrix<T> mtl_times(const complex<T> b, const Matrix<T>& A) {
	Matrix<T> C=A; C*=b; return C;
}
#endif


// Multiply Matrix*Vector'
template <typename T> Vector<T> mtl_mvtimes(const Matrix<T>& M, const Vector<T>& V) {
	assert(V.Type() == COL_VECTOR);
	assert(V.Size() == M.n); 
	Vector<T> C(M.m,COL_VECTOR);
	for(int i=1; i<=M.m; ++i) {
		T sum = 0;
		for(int k=1; k<=M.n; ++k) sum += M.mat[i][k] * V[k];
		C[i] = sum;
	}
	return C;
}


// Multiply Vector*Matrix
template <typename T> Vector<T> mtl_vmtimes(const Vector<T>& V, const Matrix<T>& M) {
	assert(V.Type() == ROW_VECTOR);
	assert(V.Size() == M.m); 
	Vector<T> C(M.n,ROW_VECTOR);
	for(int j=1; j<=M.n; ++j) {
		T sum = 0;
		for(int k=1; k<=M.m; ++k) sum += V[k] * M.mat[k][j];
		C[j] = sum;
	}
	return C;
}


// Array Multiply
template <typename T> Matrix<T>  mtl_times(const Matrix<T>& A, const Matrix<T>& B){
	assert(A.m==B.m && A.n==B.n);
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C.mat[i][j]=A.mat[i][j]*B.mat[i][j];
	return C;
}


// Array right Division A(i,j)/B(i,j)
template <typename T> Matrix<T>  mtl_rdivide(const Matrix<T>& A, const Matrix<T>& B){
	assert(A.m==B.m && A.n==B.n);
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C.mat[i][j]=A.mat[i][j]/B.mat[i][j];
	return C;
}


// Array left Division B(i,j)/A(i,j)
template <typename T> Matrix<T>  mtl_ldivide(const Matrix<T>& A, const Matrix<T>& B){
	assert(A.m==B.m && A.n==B.n);
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C.mat[i][j]=B.mat[i][j]/A.mat[i][j];
	return C;
}


// Matrix right Division A\B = pinv(A)*b
template <typename T> Matrix<T>  mtl_mldivide(const Matrix<T>& A, const Matrix<T>& B){
	Matrix<T> C(A.n,B.n);
	Vector<T> b(B.m),x(A.n);
	for(int j=1;j<=B.n;j++){
		b=B.Col(j);
		int isOK=Solve(A,b,x);
		C.SetCol(j,x);
	}
	return C;
}


// Matrix left Division A/B = A*inv(B). More precisely, A/B = (B'\A')'
template <typename T> Matrix<T>  mtl_mrdivide(const Matrix<T>& A, const Matrix<T>& B){
	return Transpose(mtl_mldivide(Transpose(B),Transpose(A)));
}


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


template <typename T> Matrix<T> mtl_power(const Matrix<T>& A, int x){
	Matrix<T> C(A.m,A.n);
	for(int i=1;i<=A.m;i++)
		for(int j=1;j<=A.n;j++) C[i][j]=Ipow(A[i][j],x);
	return C;
}