#ifndef MATRIX_H

#define MATRIX_H

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

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

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

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

	matrix operator+(DOUBLE a);	// 矩陣加常數，指每一元素加一固定的常數，產生
											// 新矩陣，原矩陣不變
	friend matrix operator+(DOUBLE a, matrix& m); // 常數加矩陣，產生新的和矩陣
	matrix& operator+=(DOUBLE a);	// 矩陣自身加常數，自身內容改變

	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) 	// 矩陣減常數，指每一元素減一固定的常數，產生
											// 新矩陣，原矩陣不變
	{ return (*this)+(-a); };
	friend matrix operator-(DOUBLE a, matrix& m); // 常數減矩陣，產生新的矩陣
	matrix& operator-=(DOUBLE a)	// 矩陣自身減常數，自身內容改變
	{ 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); // 用主高斯消元法求解線性方程的解
		// 矩陣本身為系數矩陣，m為常數向量，返回解向量
	matrix operator/(matrix m); // 左乘m的逆矩陣

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

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

	void house(buffer & b, buffer & c); // 用豪斯荷爾德變換將對稱陣變為對稱三對
											// 角陣，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列的數值，不做行列檢查
		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) { // 常數加矩陣，產生新的和矩陣
	return m+a;
};

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


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

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

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

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

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

#endif // MATRIX_H