matrix.h

矩阵运算的模板类

// -*- c++ -*-
///////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2001 Oh-Wook Kwon, all rights reserved. [e-mail]
//
//                          Easy Matrix Template Library
// 
// This Easy Matrix Template Library is provided "as is" without any express 
// or implied warranty of any kind with respect to this software. 
// In particular the authors shall not be liable for any direct, 
// indirect, special, incidental or consequential damages arising 
// in any way from use of the software.
// 
// Everyone is granted permission to copy, modify and redistribute this
// Easy Matrix Template Library, provided:
//  1.  All copies contain this copyright notice.
//  2.  All modified copies shall carry a notice stating who
//      made the last modification and the date of such modification.
//  3.  No charge is made for this software or works derived from it.  
//      This clause shall not be construed as constraining other software
//      distributed on the same medium as this software, nor is a
//      distribution fee considered a charge.
//
///////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
// Filename: Matrix.h
//
// Version: 1.0 (November 5, 2001)
//
// Programmer: Oh-Wook Kwon ([e-mail])
//
// Installed Operating Systems/Compilers:
//     Windows 32/Visual C++ 6.0
//     FreeBSD/GNU C++ 2.95.2
//
// Other Information:
//     Refer to README.txt in this directory.
//
///////////////////////////////////////////////////////////////////////////////

#ifndef _MATRIX_H_
#define _MATRIX_H_

// Indicate that ezMTL is used.
#ifndef __EZMTL__
#define __EZMTL__
#endif

#include "./Vector.h"

template <typename T> class Vector;
template <typename T> class Matrix;

// Template functions with default arguments
#ifdef __GNUC__
template <typename T> Vector<T> Sum(const Matrix<T>& A,int idx=1);
template <typename T> Vector<T> Mean(const Matrix<T>& A,int idx=1);
template <typename T> double Norm(const Matrix<T>& A, int p=2);
template <typename T> double NormF(const Matrix<T>& A);
template <typename T> Vector<T> Var(const Matrix<T>& A,int k=0);
template <typename T> double Cond(const Matrix<T>& A,int p=2);
template <typename T> Matrix<T> Inv(const Matrix<T>& A,int* retCode=0);
template <typename T> Matrix<T> PinvRecipe(const Matrix<T>& A,int* retCode=0);
template <typename T> Matrix<T> PinvMeschach(const Matrix<T>& A, int* retCode=0);
template <typename T> Matrix<T> Pinv(const Matrix<T>& A,int* retCode=0);
#else
template <typename T> Vector<T> Sum(const Matrix<T>& A,int idx);
template <typename T> Vector<T> Mean(const Matrix<T>& A,int idx);
template <typename T> double Norm(const Matrix<T>& A, int p);
template <typename T> double NormF(const Matrix<T>& A);
template <typename T> Vector<T> Var(const Matrix<T>& A,int k);
template <typename T> double Cond(const Matrix<T>& A,int p);
template <typename T> Matrix<T> Inv(const Matrix<T>& A,int* retCode);
template <typename T> Matrix<T> PinvRecipe(const Matrix<T>& A,int* retCode);
template <typename T> Matrix<T> PinvMeschach(const Matrix<T>& A, int* retCode);
template <typename T> Matrix<T> Pinv(const Matrix<T>& A,int* retCode);
#endif

// friend template functions
template <typename T> Matrix<T> mtl_plus(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> mtl_minus(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> mtl_plus(const Matrix<T>& A, const double b);
template <typename T> Matrix<T> mtl_minus(const Matrix<T>& A, const double b);
template <typename T> Matrix<T> mtl_times(const Matrix<T>& A, const double b);
template <typename T> Matrix<T> mtl_divide(const Matrix<T>& A, const double b);
template <typename T> Matrix<T> mtl_plus(const double b, const Matrix<T>& A);
template <typename T> Matrix<T> mtl_minus(const double b, const Matrix<T>& A);
template <typename T> Matrix<T> mtl_times(const double b, const Matrix<T>& A);
template <typename T> Matrix<T> mtl_divide(const double b, const Matrix<T>& A);
template <typename T> Matrix<T> mtl_times(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> mtl_mtimes(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Vector<T> mtl_mvtimes(const Matrix<T>& M, const Vector<T>& V);
template <typename T> Vector<T> mtl_vmtimes(const Vector<T>& V, const Matrix<T>& M);
template <typename T> Matrix<T> mtl_mrdivide(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> mtl_rdivide(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> mtl_mldivide(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> mtl_ldivide(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> mtl_power(const Matrix<T>& A, double x);
template <typename T> Matrix<T> mtl_power(const Matrix<T>& A, int x);
template <typename T> Matrix<T> mtl_mpower(const Matrix<T>& A, double x);
template <typename T> Matrix<T> mtl_mpower(const Matrix<T>& A, int x);
template <typename T> ostream& mtl_ostream(ostream& s, const Matrix<T>& A);
template <typename T> istream& mtl_istream(istream& s, Matrix<T>& A);

#ifndef DISABLE_COMPLEX
template <typename T> Matrix<T> mtl_plus(const Matrix<T>& A, const complex<T> b);
template <typename T> Matrix<T> mtl_minus(const Matrix<T>& A, const complex<T> b);
template <typename T> Matrix<T> mtl_times(const Matrix<T>& A, const complex<T> b);
template <typename T> Matrix<T> mtl_divide(const Matrix<T>& A, const complex<T> b);
template <typename T> Matrix<T> mtl_plus(const complex<T> b, const Matrix<T>& A);
template <typename T> Matrix<T> mtl_minus(const complex<T> b, const Matrix<T>& A);
template <typename T> Matrix<T> mtl_times(const complex<T> b, const Matrix<T>& A);
template <typename T> Matrix<T> mtl_divide(const complex<T> b, const Matrix<T>& A);
#endif

template <typename T> Vector<int> Size(const Matrix<T>& A);
template <typename T> int Size(const Matrix<T>& A,int d);
template <typename T> int NumRows(const Matrix<T>& M);
template <typename T> int NumCols(const Matrix<T>& M);
template <typename T> T Multiply3(const Vector<T>& x, const Matrix<T>& B, const Vector<T>& y);
template <typename T> T MultiplyXtAY(const Vector<T>& x, const Matrix<T>& B, const Vector<T>& y);
template <typename T> Matrix<T> ArrayMultiply(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> ArrayDivide(const Matrix<T>& A, const Matrix<T>& B);
template <typename T> Matrix<T> Abs(const Matrix<T>& A);
template <typename T> Matrix<T> Sign(const Matrix<T>& A);
template <typename T> Matrix<T> Pow(const Matrix<T>& A, double x);
template <typename T> Matrix<T> Pow(const Matrix<T>& A, int x);
template <typename T> Matrix<T> Powm(const Matrix<T>& A, double x);
template <typename T> Matrix<T> Powm(const Matrix<T>& A, int x);
template <typename T> Matrix<T> Exp(const Matrix<T>& A);
template <typename T> Matrix<T> Expm(const Matrix<T>& A);
template <typename T> Matrix<T> Log(const Matrix<T>& A);
template <typename T> Matrix<T> Logm(const Matrix<T>& A);
template <typename T> Matrix<T> Sqrt(const Matrix<T>& A);
template <typename T> Matrix<T> Sqrtm(const Matrix<T>& A);
template <typename T> Matrix<int> Int(const Matrix<T>& A);
template <typename T> Matrix<float> Float(const Matrix<T>& A);
template <typename T> Matrix<double> Double(const Matrix<T>& A);

template <typename T> int Eig(const Matrix<T>& A, Matrix<T>& EVec, Matrix<T>& EVal);
template <typename T> int EigS(const Matrix<T>& A, Matrix<T>& EVec, Matrix<T>& EVal);
template <typename T> int EigJacobi(const Matrix<T>& A, Matrix<T>& EVec, Matrix<T>& EVal);
template <typename T> int EigHHolder(const Matrix<T>& A, Matrix<T>& EVec, Matrix<T>& EVal);
template <typename T> int EigReal(const Matrix<T>& A, Matrix<T>& EVec, Matrix<T>& EVal);
template <typename T> int EigHermitian(const Matrix<T>& A, Matrix<T>& EVec, Matrix<T>& EVal);
template <typename T> int EigComplex(const Matrix<T>& A, Matrix<T>& EVec, Matrix<T>& EVal);
template <typename T> int SimDiag(const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& E, Matrix<T>& D);
template <typename T> Matrix<T> Transpose(const Matrix<T>& A);
template <typename T> Matrix<T> ArrayTranspose(const Matrix<T>& A);
template <typename T> Matrix<T> Conj(const Matrix<T>& A);
template <typename T> Matrix<T> Cov(const Matrix<T>& A);
template <typename T> int Svd(const Matrix<T> A, Matrix<T>& U, Matrix<T>& S, Matrix<T>& V);
template <typename T> int SvdRecipe(const Matrix<T> A, Matrix<T>& U, Matrix<T>& S, Matrix<T>& V);
template <typename T> int SvdMeschach(const Matrix<T>& A, Matrix<T>& U, Matrix<T>& S, Matrix<T>& V);
template <typename T> int SvdSolve(const Matrix<T>& A, const Vector<T>& b, Vector<T>& x);
template <typename T> int Schur(const Matrix<T>& A, Matrix<T>& t, Matrix<T>& Q);
template <typename T> int SchurReal(const Matrix<T>& A, Matrix<T>& t, Matrix<T>& Q);
template <typename T> int SchurComplex(const Matrix<T>& A, Matrix<T>& t, Matrix<T>& Q);
template <typename T> int Hess(const Matrix<T>& A, Matrix<T>& H, Matrix<T>& P);
template <typename T> Matrix<T> Orth(const Matrix<T>& A);
template <typename T> Matrix<T> Null(const Matrix<T>& A);
template <typename T> int Lu(const Matrix<T>& A, Matrix<T>& L, Matrix<T>& U, Matrix<T>& P);
template <typename T> int Lu(const Matrix<T>& A, Matrix<T>& L, Matrix<T>& U);
template <typename T> int LuSolve(const Matrix<T>& A, const Vector<T>& b, Vector<T>& x);
template <typename T> int Chol(const Matrix<T>& A, Matrix<T>& L);
template <typename T> int CholSolve(const Matrix<T>& A, const Vector<T>& b, Vector<T>& x);
template <typename T> int Ldl(const Matrix<T>& A, Matrix<T>& L, Matrix<T>& D);
template <typename T> int LdlSolve(const Matrix<T>& A, const Vector<T>& b, Vector<T>& x);
template <typename T> T Det(const Matrix<T>& A);
template <typename T> double LogAbsDet(const Matrix<T>& A);
template <typename T> int Qr(const Matrix<T>& A, Matrix<T>& Q, Matrix<T>& R);
template <typename T> int QrSolve(const Matrix<T>& A, const Vector<T>& b, Vector<T>& x);
template <typename T> int Qrcp(const Matrix<T>& A, Matrix<T>& QRCP, Vector<T>& diag, Vector<int>& pivot);
template <typename T> int QrcpSolve(const Matrix<T>& A, const Vector<T>& b, Vector<T>& x);
template <typename T> int Bkp(const Matrix<T>& A, Matrix<T>& BKP, Vector<int>& pivot, Vector<int>& blocks);
template <typename T> int BkpSolve(const Matrix<T>& A, const Vector<T>& b, Vector<T>& x);
template <typename T> double CondLU(const Matrix<T>& A);
template <typename T> double CondQR(const Matrix<T>& A,int p);
template <typename T> Vector<T> Diag(const Matrix<T>& A);
template <typename T> T Trace(const Matrix<T>& A);
template <typename T> Vector<T> Moment(const Matrix<T>& A, int order);
template <typename T> Vector<T> Std(const Matrix<T>& A);
template <typename T> Vector<T> Skewness(const Matrix<T>& A);
template <typename T> Vector<T> Kurtosis(const Matrix<T>& A);
template <typename T> Matrix<T> Reshape(const Matrix<T>& A);
template <typename T> int Solve(const Matrix<T>& A, const Vector<T>& b, Vector<T>& x);
template <typename T> int Rank(const Matrix<T>& A);


// ------------------------------------------------------------------------
// Matrix declaration
// ------------------------------------------------------------------------
template <typename T>
class Matrix {
private:
	// I chose to implement a matrix by a Vector of Vector because it is
	// much simpler than other methods. For sparse matrices, the method
	// wastes memory.
	int m;
	int n;
	Vector<Vector<T> > mat;
public:
	Matrix(int m_=0,int n_=0,const T& v=T());
	Matrix(int m_, int n_, const char* v_);
	Matrix(const Vector<T>& v, int m_,int n_);
	//Matrix(const Matrix<T>& A);
	virtual ~Matrix();
	//Matrix<T>&  operator=(const Matrix<T>& A);
	Matrix<T>&  operator=(const T& a);
#ifdef ALLOW_ACCESS_BY_BRACKET
	// The [] operator does not check if the arguments are in the valid range.
	// I recommend to use it only in implementation of the Matrix class.
	Vector<T>&			operator[](int i);
	const Vector<T>&	operator[](int i) const;
#endif
	Vector<T>&			operator()(int i);
	const Vector<T>&	operator()(int i) const;
	T&			operator()(int i,int j);
	const T&	operator()(int i,int j) const;
	Matrix<T>   operator+();
	Matrix<T>   operator-();
	Matrix<T>&  operator+=(const Matrix<T>& B);
	Matrix<T>&  operator-=(const Matrix<T>& B);
	Matrix<T>&  operator+=(const double b);
	Matrix<T>&  operator-=(const double b);
	Matrix<T>&  operator*=(const double b);
	Matrix<T>&  operator/=(const double b);
#ifndef DISABLE_COMPLEX
	Matrix<T>&  operator+=(const complex<T> b);
	Matrix<T>&  operator-=(const complex<T> b);
	Matrix<T>&  operator*=(const complex<T> b);
	Matrix<T>&  operator/=(const complex<T> b);
#endif
	/////////////////////////////
	// friend operator overloading
	/////////////////////////////
	// GCC-2.95-2 did not allow friend operator overloading definition outside class declaration.
	friend ostream& operator<<(ostream& s, const Matrix<T>& A) {
		return mtl_ostream(s,A);
	}
	friend ostream& operator<<(ostream& s, Matrix<T>& A) {
		return mtl_ostream(s,A);
	}
	friend istream& operator>>(istream& s, Matrix<T>& A) {
		return mtl_istream(s,A);
	}
	friend Matrix<T> operator+(const Matrix<T>& A, const Matrix<T>& B) {
		return mtl_plus(A,B);
	}
	friend Matrix<T> operator-(const Matrix<T>& A, const Matrix<T>& B) {
		return mtl_minus(A,B);
	}
	friend Matrix<T> operator+(const Matrix<T>& A, const double b) {
		return mtl_plus(A,b);
	}
	friend Matrix<T> operator-(const Matrix<T>& A, const double b) {
		return mtl_minus(A,b);
	}
	friend Matrix<T> operator*(const Matrix<T>& A, const double b) {
		return mtl_times(A,b);
	}
	friend Matrix<T> operator/(const Matrix<T>& A, const double b) {
		return mtl_divide(A,b);
	}
	friend Matrix<T> operator+(const double b, const Matrix<T>& A) {
		return mtl_plus(b,A);
	}
	friend Matrix<T> operator-(const double b, const Matrix<T>& A) {
		return mtl_minus(b,A);
	}
	friend Matrix<T> operator*(const double b, const Matrix<T>& A) {
		return mtl_times(b,A);
	}
	friend Matrix<T> operator/(const double b, const Matrix<T>& A) {
		return mtl_divide(b,A);
	}
	// Multiply Matrix * Matrix
	friend Matrix<T> operator*(const Matrix<T>& A, const Matrix<T>& B) {
		return mtl_mtimes(A,B);
	}
	// Multiply Matrix*Vector'
	friend Vector<T> operator*(const Matrix<T>& M, const Vector<T>& V) {
		return mtl_mvtimes(M,V);
	}
	// Multiply Vector*Matrix
	friend Vector<T> operator*(const Vector<T>& V, const Matrix<T>& M) {
		return mtl_vmtimes(V,M);
	}

#ifndef DISABLE_COMPLEX
	friend Matrix<T> operator+(const Matrix<T>& A, const complex<T> b) {
		return mtl_plus(A,b);
	}
	friend Matrix<T> operator-(const Matrix<T>& A, const complex<T> b) {
		return mtl_minus(A,b);
	}
	friend Matrix<T> operator*(const Matrix<T>& A, const complex<T> b) {
		return mtl_times(A,b);
	}
	friend Matrix<T> operator/(const Matrix<T>& A, const complex<T> b) {
		return mtl_divide(A,b);
	}
	friend Matrix<T> operator+(const complex<T> b, const Matrix<T>& A) {
		return mtl_plus(b,A);
	}
	friend Matrix<T> operator-(const complex<T> b, const Matrix<T>& A) {
		return mtl_minus(b,A);
	}
	friend Matrix<T> operator*(const complex<T> b, const Matrix<T>& A) {
		return mtl_times(b,A);
	}
	friend Matrix<T> operator/(const complex<T> b, const Matrix<T>& A) {
		return mtl_divide(b,A);
	}
#endif

	/////////////////////////////
	// friend functions
	/////////////////////////////
#ifdef __GNUC__
	friend Matrix<T> mtl_plus<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> mtl_minus<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> mtl_plus<T>(const Matrix<T>& A, const double b);
	friend Matrix<T> mtl_minus<T>(const Matrix<T>& A, const double b);
	friend Matrix<T> mtl_times<T>(const Matrix<T>& A, const double b);
	friend Matrix<T> mtl_divide<T>(const Matrix<T>& A, const double b);
	friend Matrix<T> mtl_plus<T>(const double b, const Matrix<T>& A);
	friend Matrix<T> mtl_minus<T>(const double b, const Matrix<T>& A);
	friend Matrix<T> mtl_times<T>(const double b, const Matrix<T>& A);
	friend Matrix<T> mtl_divide<T>(const double b, const Matrix<T>& A);
	friend Matrix<T> mtl_times<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> mtl_mtimes<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Vector<T> mtl_mvtimes<T>(const Matrix<T>& M, const Vector<T>& V);
	friend Vector<T> mtl_vmtimes<T>(const Vector<T>& V, const Matrix<T>& M);
	friend Matrix<T> mtl_mrdivide<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> mtl_rdivide<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> mtl_mldivide<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> mtl_ldivide<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> mtl_power<T>(const Matrix<T>& A, double x);
	friend Matrix<T> mtl_power<T>(const Matrix<T>& A, int x);
	friend Matrix<T> mtl_mpower<T>(const Matrix<T>& A, double x);
	friend Matrix<T> mtl_mpower<T>(const Matrix<T>& A, int x);
	friend ostream& mtl_ostream<T>(ostream& s, const Matrix<T>& A);
	friend istream& mtl_istream<T>(istream& s, Matrix<T>& A);

#ifndef DISABLE_COMPLEX
	friend Matrix<T> mtl_plus<T>(const Matrix<T>& A, const complex<T> b);
	friend Matrix<T> mtl_minus<T>(const Matrix<T>& A, const complex<T> b);
	friend Matrix<T> mtl_times<T>(const Matrix<T>& A, const complex<T> b);
	friend Matrix<T> mtl_divide<T>(const Matrix<T>& A, const complex<T> b);
	friend Matrix<T> mtl_plus<T>(const complex<T> b, const Matrix<T>& A);
	friend Matrix<T> mtl_minus<T>(const complex<T> b, const Matrix<T>& A);
	friend Matrix<T> mtl_times<T>(const complex<T> b, const Matrix<T>& A);
	friend Matrix<T> mtl_divide<T>(const complex<T> b, const Matrix<T>& A);
#endif

	friend Vector<int> Size<T>(const Matrix<T>& A);
	friend int Size<T>(const Matrix<T>& A,int d);
	friend int NumRows<T>(const Matrix<T>& M);
	friend int NumCols<T>(const Matrix<T>& M);
	friend T Multiply3<T>(const Vector<T>& x, const Matrix<T>& B, const Vector<T>& y);
	friend T MultiplyXtAY<T>(const Vector<T>& x, const Matrix<T>& B, const Vector<T>& y);
 	friend Matrix<T> ArrayMultiply<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> ArrayDivide<T>(const Matrix<T>& A, const Matrix<T>& B);
	friend Matrix<T> Abs<T>(const Matrix<T>& A);
	friend Matrix<T> Sign<T>(const Matrix<T>& A);
	friend Matrix<T> Pow<T>(const Matrix<T>& A, double x);
	friend Matrix<T> Pow<T>(const Matrix<T>& A, int x);
	friend Matrix<T> Powm<T>(const Matrix<T>& A, double x);
	friend Matrix<T> Powm<T>(const Matrix<T>& A, int x);
	friend Matrix<T> Exp<T>(const Matrix<T>& A);
	friend Matrix<T> Expm<T>(const Matrix<T>& A);
	friend Matrix<T> Log<T>(const Matrix<T>& A);
	friend Matrix<T> Logm<T>(const Matrix<T>& A);
	friend Matrix<T> Sqrt<T>(const Matrix<T>& A);
	friend Matrix<T> Sqrtm<T>(const Matrix<T>& A);
	friend Matrix<int> Int<T>(const Matrix<T>& A);
	friend Matrix<float> Float<T>(const Matrix<T>& A);
	friend Matrix<double> Double<T>(const Matrix<T>& A);
	friend Matrix<T> Transpose<T>(const Matrix<T>& A);