matrix.h

矩阵运算的模板类

	friend Matrix<T> ArrayTranspose<T>(const Matrix<T>& A);
	friend Matrix<T> Conj<T>(const Matrix<T>& A);
	friend Matrix<T> Cov<T>(const Matrix<T>& A);
	friend Matrix<T> Orth<T>(const Matrix<T>& A);
	friend Matrix<T> Null<T>(const Matrix<T>& A);

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

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

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

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

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

#endif

	/////////////////////////////
	// member functions
	/////////////////////////////
	Matrix<T> t() const;
	Matrix<T> ct() const;
	Matrix<T> i() const;
	int p(string& name) const;
	int NormRows(int i=0);
	int NormCols(int j=0);
	Vector<T> Sum_(int idx=1) const;
	int Print(ostream& s, const char *name=0);
	int Print(ostream& s, string& name);
	int Read(istream& s, string& name);
	int Read(istream& s, const char *name=0);
	int Resize(int new_m, int new_n);
	Vector<T> OutputAsRow();
	Vector<T> OutputAsCol();
	Vector<T> Row(int m1) const;
	Matrix<T> Rows(int m1,int m2) const;
	Vector<T> Col(int n1) const;
	Matrix<T> Cols(int n1,int n2) const;
	Matrix<T> SubMatrix(int m1,int m2,int n1,int n2) const;
	T Scalar(const Matrix<T>& A) const;
	int SetRow(int rowX,const Vector<T>& A);
	int SetCol(int colX,const Vector<T>& A);
	int SetIdentity();
	int Save(const string& fname);
	int SaveMatlab(const string& fname);
	int Input(const Vector<T>& V);
	int Input(const Matrix<T>& A);

	bool IsSquare() const;
	bool IsSymmetric(double tolerance=1e-6) const;
	bool IsZero(double tolerance=1e-6) const;
	bool IsEye(double tolerance=1e-6) const;
	bool IsNormal(double tolerance=1e-6) const;
	bool IsOrthogonal(double tolerance=1e-6) const;
	bool IsDiagonal(double tolerance=1e-6) const;
	bool IsPD(double tolerance=1e-6) const;
	bool IsND(double tolerance=1e-6) const;
	bool IsLowerTri(double tolerance=1e-6) const;
	bool IsUpperTri(double tolerance=1e-6) const;
	bool IsTri(double tolerance=1e-6) const;
	bool IsUnitary(double tolerance=1e-6) const;
	bool IsHermitian(double tolerance=1e-6) const;
	bool IsReal() const;
	bool IsComplex() const;
	bool IsDouble() const;

	void FromInt(const Matrix<int>& A);
	void FromFloat(const Matrix<float>& A);
	void FromDouble(const Matrix<double>& A);
	Matrix<int> ToInt() const;
	Matrix<float> ToFloat() const;
	Matrix<double> ToDouble() const;
#ifndef DISABLE_COMPLEX
	void FromCDouble(const Matrix<complex<double> >& A);
	Matrix<complex<double> > ToCDouble() const;
#endif
	Matrix<double> Real() const;
	Matrix<double> Imag() const;

public:
	/////////////////////////////
	// exported static functions
	/////////////////////////////
	static Matrix<T> Randn(int m,int n);
	static Matrix<T> Eye(int m_,int n_);
	static Matrix<T> Zeros(int m_,int n_);
	static Matrix<T> CRandn(int m,int n);

	/////////////////////////////
	// local static functions
	/////////////////////////////
	static double normnxn(const Matrix<T>& A);
	static int UpperBackSubstitute(const Matrix<T>& A,Vector<T>& b);
	static int LowerBackSubstitute(const Matrix<T>& A,Vector<T>& b);
	static int cleanNearZero(Matrix<T>& A,double tolerance=1e-6);
	static int svdClean(Vector<T>& s, double tolerance=1e-6);
	static int pinv(Matrix<T>& U, Matrix<T>& S, Matrix<T>& V, Matrix<T>& pinvA);
	static int sortEigen(Matrix<T>& EVec,Matrix<T>& EVal);
	static int sortEigenCompact(Matrix<T>& EVec,Matrix<T>& EVal);
	static int Compact2Complex(Matrix<T>& EVec, Matrix<T>& EVal, Matrix<complex<T> >& CVec, Matrix<complex<T> >& CVal);
	static int Complex2Compact(Matrix<complex<T> >& CVec, Matrix<complex<T> >& CVal, Matrix<T>& EVec, Matrix<T>& EVal);
	static double pythagoras(const T a, const T b);
	static double pythagoras3(const T a, const T b, const T c);
	static int jacobiTransform(Matrix<T>& A, Vector<T>& d, Matrix<T>& v,int *rotN);
	static void jacobiRotate(Matrix<T>& a,int i,int j,int k,int l,double s,double tau);
	static Matrix<T> ExpmReal(const Matrix<T>& A);
	static Matrix<T> ExpmComplex(const Matrix<T>& A);
	static Matrix<T> LogmReal(const Matrix<T>& A);
	static Matrix<T> LogmComplex(const Matrix<T>& A);

	/////////////////////////////
	// From Meschach library
	/////////////////////////////
	// mes_eig.h: Sensitive to precision. Use double precision.
	static T mes_in_prod(const Vector<T>& a,const Vector<T>& b);
	static T mes_in_prod_offset(const Vector<T>& a,const Vector<T>& b,int i0);
	static double mes_norm_inf(const Vector<T>& x);
	static double mes_norm2(const Vector<T>& x);
	static double mes_norm2_offset(const Vector<T>& a,int i0);	
	static void	mes_get_row(const Matrix<T>& mat,int row,Vector<T>& vec);
	static void	mes_get_col(const Matrix<T>& mat,int col,Vector<T>& vec);
	static void	mes_set_row(Matrix<T>& mat,int row,Vector<T>& vec);
	static void	mes_set_col(Matrix<T>& mat,int col,Vector<T>& vec);	
	static void	mes_copy_offset(const Vector<T>& in,Vector<T>& out,int i0);
	static int mes_trieig(Vector<T>& a,Vector<T>& b,Matrix<T>& Q,int maxIter=100);
	static int mes_symmeig(const Matrix<T>& A,Matrix<T>& Q,Vector<T>& out);
	static int mes_hhldr3cols(Matrix<T>& A, int k, int j0, double beta, double nu1, double nu2, double nu3);
	static int mes_hhldr3rows(Matrix<T>& A, int k, int i0, double beta, double nu1, double nu2, double nu3);
	static int mes_schur(Matrix<T>& A,Matrix<T>& Q,int maxIter=100);
	static int mes_schur_evals(const Matrix<T>& t,Vector<T>& real_pt,Vector<T>& imag_pt);
	static int mes_schur_vecs(const Matrix<T>& t,const Matrix<T>& Q,Matrix<T>& X_re,Matrix<T>& X_im);
	static int mes_Hfactor(Matrix<T>& A, Vector<T>& diag, Vector<T>& beta);
	static int mes_makeHQ(const Matrix<T>& H, Vector<T>& diag, Vector<T>& beta, Matrix<T>& Qout);
	static int mes_makeH(const Matrix<T>& H,Matrix<T>& Hout);
	static int mes_hhvec(Vector<T>& vec,int i0,double* beta,Vector<T>& out,T* newval);
	static int mes_hhtrvec(Vector<T>& hh,double beta,int i0,Vector<T>& in,Vector<T>& out);
	static int mes_hhtrrows(Matrix<T>& M,int i0,int j0,Vector<T>& hh,double beta);
	static int mes_hhtrcols(Matrix<T>& M,int i0,int j0,Vector<T>& hh,double beta);
	static int mes_rot_vec(Vector<T>& x,int i,int k,double c,T s,Vector<T>& out);
	static int mes_rot_rows(Matrix<T>& mat,int i,int k,double c,T s,Matrix<T>& out);
	static int mes_rot_cols(Matrix<T>& mat,int i,int k,double c,T s,Matrix<T>& out);
	static void	mes_hhldr3(double x, double y, double z, double *nu1, double *beta, double *newval);
	static void	mes_givens(T x,T y,double* c,T* s);
	// for complex matrices
	static int mes_zHfactor(Matrix<T>& A, Vector<T>& diag);
	static int mes_zHQunpack(Matrix<T>& HQ,Vector<T>& diag,Matrix<T>& Q,Matrix<T>& H);
	static int mes_zschur(Matrix<T>& A,Matrix<T>& Q,int maxIter=100);
	static int mes_zschur_evals(const Matrix<T>& t,Vector<T>& EVal);
	static int mes_zschur_vecs(const Matrix<T>& t,const Matrix<T>& Q,Matrix<T>& EVec);

	// mes_rand.h
	static double mes_mrand();
	static int mes_mrandlist(T *a, int len);
	static void mes_smrand(int seed);
	static void	mes_v_rand(Vector<T>& x);
	// mes_svd.h: Sensitive to precision. Use double precision.
	static int mes_fixsvd(Vector<T>& dd, Matrix<T>& U, Matrix<T>& V);
	static int mes_bisvd(Vector<T>& dd, Vector<T>& ff, Matrix<T>& U, Matrix<T>& V,int maxIter=100);
	static int mes_svd(const Matrix<T>& A, Matrix<T>& U, Matrix<T>& V, Vector<T>& dd);
	static int mes_bifactor(Matrix<T>& A, Matrix<T>& U, Matrix<T>& V);
	static int mes_rotsvd(Vector<T>& dd, Matrix<T>& U, Matrix<T>& V);
	// mes_solve.h
	static int mes_CHfactor(Matrix<T>& A);
	static int mes_CHsolve(const Matrix<T>& A,const Vector<T>& b,Vector<T>& x);
	static int mes_LDLfactor(Matrix<T>& A);
	static int mes_LDLsolve(const Matrix<T>& LDL,const Vector<T>& b,Vector<T>& x);
	static int mes_Usolve(const Matrix<T>& matrix,const Vector<T>& b,Vector<T>& out,double diag);
	static int mes_Lsolve(const Matrix<T>& matrix,const Vector<T>& b,Vector<T>& out,double diag);
	static int mes_UTsolve(const Matrix<T>& U,const Vector<T>& b,Vector<T>& out,double diag);
	static int mes_Dsolve(const Matrix<T>& A,const Vector<T>& b,Vector<T>& x);
	static int mes_LTsolve(const Matrix<T>& L,const Vector<T>& b,Vector<T>& out,double diag);
	// mes_factor.h
	static int mes_LUfactor(Matrix<T>& A,Vector<int>& pivot);
	static int mes_LUsolve(Matrix<T>& A,Vector<int>& pivot,const Vector<T>& b,Vector<T>& x);
	static int mes_LUTsolve(Matrix<T>& LU,Vector<int>& pivot,const Vector<T>& b,Vector<T>& x);
	static int mes_inverse(const Matrix<T>& A,Matrix<T>& out);
	static double mes_LUcondest(Matrix<T>& LU,Vector<int>& pivot);
	static int mes_QRfactor(Matrix<T>& A,Vector<T>& diag);
	static int mes_QRCPfactor(Matrix<T>& A,Vector<T>& diag,Vector<int>& px);
	static int mes_Qsolve_(Matrix<T>& QR,Vector<T>& diag,const Vector<T>& b,Vector<T>& x,Vector<T>& tmp);
	static int mes_makeQ(Matrix<T>& QR,Vector<T>& diag,Matrix<T>& Qout);
	static int mes_makeR(Matrix<T>& QR,Matrix<T>& Rout);
	static int mes_QRsolve(Matrix<T>& QR,Vector<T>& diag,const Vector<T>& b,Vector<T>& x);
	static int mes_QRCPsolve(Matrix<T>& QR,Vector<T>& diag,Vector<int>& pivot,const Vector<T>& b,Vector<T>& x);
	static int mes_Umlt(const Matrix<T>& U,const Vector<T>& x,Vector<T>& out);
	static int mes_UTmlt(const Matrix<T>& U,const Vector<T>& x,Vector<T>& out);
	static int mes_QRTsolve(const Matrix<T>& A,Vector<T>& diag,Vector<T>& c,Vector<T>& sc);
	static double mes_QRcondest(const Matrix<T>& QR,int p=2);
	static void mes_interchange(Matrix<T>& A,int i,int j);
	static int mes_BKPfactor(Matrix<T>& A,Vector<int>& pivot,Vector<int>& blocks);
	static int mes_BKPsolve(const Matrix<T>& A,Vector<int>& pivot,Vector<int>& block,const Vector<T>& b,Vector<T>& x);
	// mes_perm.h
	static int mes_px_inv(Vector<int>& px,Vector<int>& out);
	static int mes_px_mlt(Vector<int>& px1,Vector<int>& px2,Vector<int>& out);
	static int mes_px_vec(Vector<int>& px,Vector<T>& vector,Vector<T>& out);
	static int mes_pxinv_vec(Vector<int>& px,Vector<T>& x,Vector<T>& out);
	static int mes_px_transp(Vector<int>& px,int i1,int i2);
	static int mes_myqsort(int *a,int num);
	static int mes_px_sign(Vector<int>& px);
	static int mes_px_cols(Vector<int>& px,Matrix<T>& A,Matrix<T>& out);
	static int mes_px_rows(Vector<int>& px,Matrix<T>& A,Matrix<T>& out);
	// mes_math.h
	static int mes_exp(Matrix<T>& A,double eps,Matrix<T>& out,int* q_out,int* j_out);

	/////////////////////////////
	// From Numerical Recipes in C
	/////////////////////////////
	static int nrc_choldcmp(Matrix<T>& A, Vector<T>& d);
	static int nrc_cholsolv(const Matrix<T>& L, const Vector<T>& d, const Vector<T>& b, Vector<T>& x);
	static int nrc_ludcmp(Matrix<T>& A,Vector<int>& index,double* d,Vector<int>& index2);
	static int nrc_lusolv(const Matrix<T>& LU,const Vector<int>& index,const Vector<T>& b,Vector<T>& x);
	static int nrc_qrdcmp(Matrix<T>& A,Vector<T>& c,Vector<T>& d);
	static int nrc_qrsolv(const Matrix<T>& QR,const Vector<T>& c,const Vector<T>& d,const Vector<T>& b, Vector<T>& x);
	static int nrc_svdcmp(Matrix<T>& A, Vector<T>& S, Matrix<T>& V,int maxIter=100);
	static int nrc_svsolv(Matrix<T>& U, Vector<T>& S, Matrix<T>& V,const Vector<T>& b,Vector<T>& x);
	static int nrc_tred2(Matrix<T>& A,Vector<T>& d,Vector<T>& e,int compEVec);
	static int nrc_tqli(Vector<T>& d,Vector<T>& e,Matrix<T>& Z,int compEVec);

};

typedef Matrix<float> FMatrix;
typedef Matrix<double> DMatrix;
typedef Matrix<FComplex> CFMatrix;
typedef Matrix<DComplex> CDMatrix;


// ------------------------------------------------------------------------
// Matrix implementation
// ------------------------------------------------------------------------
#ifndef local_max
#define local_max(a, b)  (((a) > (b)) ? (a) : (b))
#endif
#ifndef local_min
#define local_min(a, b)  (((a) < (b)) ? (a) : (b)) 
#endif


/////////////////////////////