#ifndef MATRIX_H

#define MATRIX_H

// 定義適合數(shù)學運算的實數(shù)matrix類，數(shù)據將存放在buffer類中

#include <iostream.h>
#include <iomanip.h>
#include <stdio.h>

// buffer.h包含實數(shù)緩存類buffer的定義
#include "buffer.h"

// 實數(shù)矩陣類，每個矩陣有rownum行，colnum列，總共rownum乘colnum個實數(shù)
// 下標一律從0開始起算，即存在第0行第0列
class matrix {
 public:
	buffer * buf;  // 緩存區(qū)指針，指向存放實數(shù)的地方，可以是任何媒體
		// buffer類是一抽象類，具體的子類可以將數(shù)據存放在內存或者
		// 臨時文件中，用戶還可以設計自己的緩存區(qū)子類
	size_t rownum, colnum;  // 矩陣的行數(shù)和列數(shù)
	unsigned	istrans:1;		// 轉置標志，0為缺省，1為轉置
	unsigned isneg:1;			// 為負標志，0為非負，1為取負
	unsigned issym:1;			// 對稱標志，0為非對稱，1為對稱，不具有約束力
	// 在下面的所有構造函數(shù)中，緩存區(qū)指針b通常不設，缺省為0，這樣程序將根據
	// 自動申請缺省的緩存區(qū)類型。如果用戶定義了特殊的緩顧慮區(qū)類型，可以
	// 先申請后在構造函數(shù)中指明
	matrix(buffer * b=0); // 缺省構造函數(shù)，產生0行0列空矩陣
	matrix(size_t n, buffer * b=0); // 產生nX1列向量
	matrix(size_t r, size_t c, buffer * b=0); // 構造函數(shù)，產生r行c列矩陣
	matrix(matrix& m); //拷貝構造函數(shù)，產生的矩陣內容同m相同，
		// 新產生的矩陣并不申請新的緩存，而是指向與m同樣的緩存，
		// 只是它們共同指向的緩存的引用數(shù)增1，這樣可以節(jié)省存儲資源
	matrix(const char * filename, buffer * b=0); /* 文件構造函數(shù)，
		使用數(shù)據文件對矩陣初使化，數(shù)據文件為文本文件，格式：
		以小寫的matrix關鍵詞打頭，然后是空格，再跟矩陣的行數(shù)，再跟
		空格，再跟列數(shù)，再跟空格，接著是內容，
		內容按行分，首先是行號，從1開始計數(shù)，行號后面緊跟著冒號：，
		然后是這行的數(shù)據，各個數(shù)據之間都用空格或者換行隔開 */
	matrix(void * data, size_t r, size_t c=1, buffer * b=0); /* 數(shù)據構造函數(shù)，
		行數(shù)為r,列數(shù)為c,data必須指向一存放實數(shù)的內存，按先列后行的
		順序進行存放。構造函數(shù)進行數(shù)據拷貝，因此在構造函數(shù)完成后,
		data指向的數(shù)據區(qū)可以釋放 */
	DOUBLE operator()(size_t r, size_t c); // 重載函數(shù)運算符()返回第r行第c列的實數(shù)，
	DOUBLE operator()(size_t r){return (*this)(r,0);}; // 取回第r行0列的實數(shù)
	virtual void set(size_t r, size_t c, DOUBLE v) {
		size_t l;
		v = isneg ? -v : v;
		l = istrans ? r+c*rownum : r*colnum+c;
		buf = buf->set(l, v); }; // 設置第r行第c列的值為v，
		// 在設置時如果緩存的引用數(shù)大于1，產生新的緩存
		// 因此每次設置返回的緩存區(qū)指針都要再賦給buf。
	void set(size_t r, DOUBLE v) {
		set(r,0,v);};	// 設置第r行0列的值為v,用作列向量
	virtual ~matrix(){ // 析構函數(shù)，如果緩存的引用數(shù)大于1，說明還有其它的
			// 矩陣在用此緩存，則不釋放緩存，只是將引用數(shù)減1
		buf->refnum--;
		if(!buf->refnum) delete buf;};
	matrix& operator=(matrix& m); // 重載賦值運算符，使矩陣的內容等于另一個矩陣
	matrix& operator=(DOUBLE a); // 通過賦值運算符將矩陣所有元素設為a
	matrix operator*(DOUBLE a);	// 重載乘法運算符，實現(xiàn)矩陣與數(shù)相乘，矩陣在前
				// 數(shù)在后
	matrix& operator*=(DOUBLE a);  // 重載*=運算符，功能同上，但*產生了一個新
			// 的矩陣，原矩陣不變，而本運算符修改了原矩陣
	friend matrix operator*(DOUBLE a, matrix& m);
			// 數(shù)乘矩陣，但數(shù)在前矩陣在后
	matrix operator+(matrix& m); // 矩陣相加，原二矩陣不變并產生一個新的和矩陣
	matrix& operator+=(matrix &m); // 矩陣相加，不產生新矩陣，有一矩陣的內容改變?yōu)楹?
	matrix operator+(DOUBLE a);	// 矩陣加常數(shù)，指每一元素加一固定的常數(shù)，產生
											// 新矩陣，原矩陣不變
	friend matrix operator+(DOUBLE a, matrix& m); // 常數(shù)加矩陣，產生新的和矩陣
	matrix& operator+=(DOUBLE a);	// 矩陣自身加常數(shù)，自身內容改變

	matrix operator*(matrix& m); // 矩陣相乘，產生出新的矩陣
	matrix& operator*=(matrix& m); // 矩陣相乘，結果修改原矩陣

	matrix operator-(); // 矩陣求負，產生新的負矩陣
	matrix & trans(){
		size_t r = rownum; rownum = colnum; colnum = r;
		istrans = !istrans; return (*this);}; // 矩陣自身轉置
	matrix t(); // 矩陣轉置，產生新的矩陣
	matrix & neg(){isneg = !isneg; return (*this);}; // 將自己求負

	matrix operator-(matrix& m); // 矩陣相減，產生新的矩陣
	matrix& operator-=(matrix& m); // 矩陣相減，結果修改原矩陣
	matrix operator-(DOUBLE a) 	// 矩陣減常數(shù)，指每一元素減一固定的常數(shù)，產生
											// 新矩陣，原矩陣不變
	{ return (*this)+(-a); };
	friend matrix operator-(DOUBLE a, matrix& m); // 常數(shù)減矩陣，產生新的矩陣
	matrix& operator-=(DOUBLE a)	// 矩陣自身減常數(shù)，自身內容改變
	{ return operator+=(-a);};

	int isnear(matrix& m, double e = defaulterr); // 檢查兩個矩陣是否近似相等,
								// 如近似相等則返回1，否則返回0，e為允許誤差
	int isnearunit(double e = defaulterr); // 檢查矩陣是否近似為單位矩陣
								// 如是則返回1，否則返回0，e為允許誤差

	matrix row(size_t r); // 提取第r行行向量
	matrix col(size_t c); // 提取第c列列向量

	void checksym();	// 檢查和嘗試調整矩陣到對稱矩陣

	void swapr(size_t r1, size_t r2, size_t k=0); // 交換兩行,只交換k列以上
	void swapc(size_t c1, size_t c2, size_t k=0); // 交換兩列,只交換k行以上
	void swap(size_t r1, size_t c1, size_t r2, size_t c2){ // 交換兩個元素
		DOUBLE a; a=(*this)(r1,c1);set(r1,c1,(*this)(r2,c2));
		set(r2,c2,a);};
	DOUBLE maxabs(size_t &r, size_t &c, size_t k=0);
		 // 求主元值和位置，即k行k列之后的最大元素，
		// 元素所在的行列將放在r,c中，返回此元素值
	size_t zgsxy(matrix & m, int fn=0);	// 主高斯消元運算
	matrix& operator/=(matrix m); // 用主高斯消元法求解線性方程的解
		// 矩陣本身為系數(shù)矩陣，m為常數(shù)向量，返回解向量
	matrix operator/(matrix m); // 左乘m的逆矩陣

	matrix& inv();		// 用全選主元高斯-約當法求逆矩陣
	matrix operator~();	// 求逆矩陣，但產生新矩陣
	friend matrix operator/(DOUBLE a, matrix& m); // 求逆矩陣再乘常數(shù)
	matrix& operator/=(DOUBLE a); // 所有元素乘a的倒數(shù)，自身改變
	matrix operator/(DOUBLE a); // 所有元素乘a的倒數(shù)，產生新的矩陣
	matrix cholesky(matrix& d);	// 用喬里斯基分解法求對稱正定陣的線性
		// 方程組ax=d返回方程組的解，本身為a，改變?yōu)榉纸怅噓,d不變

	DOUBLE det(DOUBLE err=defaulterr);		// 求行列式的值
	size_t rank(DOUBLE err=defaulterr);		// 求矩陣的秩

	void house(buffer & b, buffer & c); // 用豪斯荷爾德變換將對稱陣變?yōu)閷ΨQ三對
											// 角陣，b返回主對角線元素，c返回次對角線元素
	void trieigen(buffer& b, buffer& c, size_t l=600, DOUBLE eps=defaulterr);
						// 計算三對角陣的全部特征值與特征向量
	matrix eigen(matrix & eigenvalue, DOUBLE eps=defaulterr, size_t l=600);
		 // 計算對稱陣的全部特征向量和特征值
		// 返回按列排放的特征向量，而eigenvalue將返回一維矩陣，為所有特征值
		// 組成的列向量
	int ispositive();		// 判定矩陣是否對稱非負定陣，如是返回1，否則返回0

	void hessenberg();	// 將一般實矩陣約化為赫申伯格矩陣
	void qreigen(matrix & a, matrix & b, size_t l=600, DOUBLE eps=defaulterr);
		// 求一般實矩陣的所有特征根
		// a和b均返回rownum行一列的列向量矩陣，返回所有特征根的實部和虛部

	friend ostream& operator<<(ostream& o, matrix& m);
	friend istream& operator>>(istream& in, matrix& m);
	virtual DOUBLE value(size_t r, size_t c){ // 返回第r行c列的數(shù)值，不做行列檢查
		DOUBLE a;
		a = istrans ? (*buf)[r+c*rownum] : (*buf)[r*colnum + c];
		return isneg ? -a : a;
	};
 protected:
	DOUBLE detandrank(size_t & r, DOUBLE err);	// 求方陣的行列式和秩
};

inline matrix operator*(DOUBLE a, matrix& m) {
	return m*a;
};

inline matrix operator+(DOUBLE a, matrix& m) { // 常數(shù)加矩陣，產生新的和矩陣
	return m+a;
};

matrix operator-(DOUBLE a, matrix& m); // 常數(shù)減矩陣，產生新的矩陣
matrix operator/(DOUBLE a, matrix& m); // 求逆矩陣再乘常數(shù)


ostream& operator<<(ostream& o, matrix& m);
istream& operator>>(istream& in, matrix& m);

matrix unit(size_t n);	// 產生n階單位矩陣

DOUBLE gassrand(int rr=0);	// 返回一零均值單位方差的正態(tài)分布隨機數(shù)
			// rr為種子

class gassvector //返回零均值給定協(xié)方差陣的正態(tài)隨機向量的類
{
 public:
	matrix a;	// a是增益矩陣，由協(xié)方差矩陣算出
	gassvector(matrix & r);	//r必須是正定對稱陣，為正態(tài)隨機向量的協(xié)方差
	matrix operator()(matrix & r); // 返回給定協(xié)方差矩陣的正態(tài)隨機向量
	matrix operator()();	// 返回已設定的協(xié)方差矩陣的正態(tài)隨機向量
	void setdata(matrix & r); // 根據協(xié)方差矩陣設置增益矩陣
};

DOUBLE lineopt(matrix& a,matrix& b, matrix& c, matrix & x);
	// 線性規(guī)劃最優(yōu)點尋找程序，a為mXn不等式約束條件左端系數(shù)矩陣，b為不等式約束
	// 條件的右端項，為m維向量，c為目標函數(shù)系數(shù)，n維向量，x返回極小點，為n維向量

#endif // MATRIX_H