#include <stdio.h>
#include <math.h>
#include <string.h>
#include <fstream.h>
// 本程序實現matrix類
// 對matrix類的定義
#include "matrix.h"

matrix::matrix(buffer * b): // 缺省構造函數，產生0行0列空矩陣
	rownum(0),colnum(0),istrans(0),isneg(0),issym(0)
{
		if(b){ // 如用戶給出b，則使用b作為數據緩存
			buf=b;
			buf->alloc(0);
			}
		else // 否則，產生一個新的緩存，類別是當前缺省的緩存類
			buf = getnewbuffer(0);
}

matrix::matrix(size_t n, buffer * b): // 產生n階單位方陣
		rownum(n),colnum(1),istrans(0),isneg(0),issym(0)
{
		if(b){ // 如用戶給出b，則使用b作為數據緩存
			buf=b;
			buf->alloc(n);
			}
		else // 否則，產生一個新的緩存，類別是當前缺省的緩存類
			buf = getnewbuffer(n);
}

matrix unit(size_t n)	// 產生n階單位矩陣
{
	if(n==0) throw TMESSAGE("n must larger than 0\n");
	matrix m(n,n);
	for(size_t i=0; i<n; i++)
	for(size_t j=0; j<n; j++)
		if(i==j) m.set(i,j,1.0);
		else m.set(i,j,0.0);
	m.issym = 1;
	return m;
}


matrix::matrix(size_t r, size_t c, buffer * b):
		rownum(r),colnum(c),istrans(0),isneg(0),issym(0)
{
		if(b){ // 如用戶給出b，則使用b作為數據緩存
			buf=b;
			buf->alloc(r*c);
			}
		else // 否則，產生一個新的緩存，類別是當前缺省的緩存類
			buf = getnewbuffer(r*c);
}

matrix::matrix(matrix& m) // 拷貝構造函數
{
	rownum = m.rownum;	// 拷貝行數與列數
	colnum = m.colnum;
	istrans = m.istrans;	// 拷貝轉置標記
	isneg = m.isneg;		// 拷貝負值標記
	issym = m.issym;		// 拷貝對稱標記
	buf = m.buf;		// 設置指向相同緩存的指針
	buf->refnum++;	// 緩存的引用數增1
}

matrix::matrix(const char * filename, buffer * b): // 從數據文件構造矩陣
	istrans(0), isneg(0), issym(0)
{
	char label[10];
	ifstream in(filename); // 打開文件流
	in >> label; // 文件開始必須有matrix關鍵詞
	if(strcmp(label, "matrix")!=0) throw TMESSAGE("format error!");
	in >> rownum >> colnum; // 讀取行數和列數
	if(!in.good()) throw TMESSAGE("read file error!");
	// 申請或產生緩存
	if(b) { buf=b;
		buf->alloc(rownum*colnum);
	}
	else buf = getnewbuffer(rownum*colnum);
	size_t line;
	for(size_t i=0; i<rownum; i++) { // 依次按行讀取
		in >> line; // 讀行號
		if(line != i+1) throw TMESSAGE("format error!");
		in.width(sizeof(label));
		in >> label; // 行號后跟冒號
		if(label[0] != ':') throw TMESSAGE("format error!");
		DOUBLE a;
		for(size_t j=0; j<colnum; j++) { // 讀取一行數據
			in >> a;
			set(i,j,a);
		}
		if(!in.good()) throw TMESSAGE("read file error!");
	}
	checksym(); // 檢查是否為對稱陣
}

matrix::matrix(void * data, size_t r, size_t c, buffer * b):
		rownum(r),colnum(c),istrans(0),isneg(0),issym(0) // 數據構造函數
{
	if(b){
		buf=b;
		buf->alloc(r*c);
		}
	else
		buf = getnewbuffer(r*c);
	DOUBLE * d = (DOUBLE *)data;
	for(size_t i=0; i<r*c; i++)  // 這里進行數據拷貝，因此原數據的內存是可以釋放的
		buf->set(i,d[i]);
	checksym();	// 檢查是否為對稱陣
}

DOUBLE matrix::operator()(size_t r, size_t c)
{
	if(r>= rownum || c>= colnum)
		throw TMESSAGE("Out range!");
	return value(r,c);
}

matrix& matrix::operator=(matrix& m)  // 賦值重載
{
	rownum = m.rownum; //  行數和列數的拷貝
	colnum = m.colnum;
	istrans = m.istrans; // 轉置標志的拷貝
	isneg = m.isneg;	// 取負標志的拷貝
	issym = m.issym;	// 對稱標志的拷貝
	if(buf == m.buf)  // 如果原緩存與m的緩存一樣，則返回
		return (*this);
	buf->refnum--; // 原緩存不同，則原緩存的引用數減1
	if(!buf->refnum)delete buf; // 減1后的原緩存如果引用數為0，則刪除原緩存
	buf = m.buf; // 將原緩存指針指向m的緩存
	buf->refnum++; // 新緩存的引用數增1
	checksym();	// 檢查是否為對稱陣
	return (*this); // 返回自己的引用
}

matrix& matrix::operator=(DOUBLE a) // 通過賦值運算符將矩陣所有元素設為a
{
	if(rownum == 0 || colnum == 0) return (*this);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,a);
	if(rownum == colnum)issym = 1;
	return (*this);
}

matrix matrix::operator-() // 矩陣求負，產生負矩陣
{
	matrix mm(*this);
	mm.neg();
	return mm;
}

ostream& operator<<(ostream& o, matrix& m) // 流輸出運算符
{
	// 先輸出關鍵字matrix,然后是行號，制表符，列號，換行
	o << "matrix " << m.rownum << '\t' << m.colnum << endl;
	for(size_t i=0; i<m.rownum; i++) { // 依次輸出各行
		o<< (i+1) <<':'; // 行號后跟冒號。注意輸出的行號是從1開始
				// 是內部編程的行號加1
		size_t k=8; // 后面輸入一行的內容，每次一行數字超過八個時
			// 就換一行
		for(size_t j=0; j<m.colnum; j++) {
			o<<'\t'<<m.value(i,j);
			if(--k==0) {
				k=8;
				o<<endl;
			}
		}
		o<<endl; // 每行結束時也換行
	}
	return o;
}

istream& operator>>(istream& in, matrix& m) // 流輸入運算符
{
	char label[10];
	in.width(sizeof(label));
	in >> label; // 輸入matrix關鍵字
	if(strcmp(label, "matrix")!=0)
		throw TMESSAGE("format error!\n");
	in >> m.rownum >> m.colnum; // 輸入行號和列號
	if(!in.good()) throw TMESSAGE("read file error!");
	m.buf->refnum--; // 原緩存引用數減1
	if(!m.buf->refnum) delete m.buf; // 如原緩存引用數為0，則刪去原緩存
	m.isneg = m.istrans = 0; // 轉置和取負標志清零
	m.buf = getnewbuffer(m.rownum*m.colnum); // 按缺省子類產生新緩存
	size_t line; // 按矩陣行輸入
	for(size_t i=0; i<m.rownum; i++) {
		in >> line; // 先輸入行號
		if(line != i+1) throw TMESSAGE("format error!\n");
		in.width(sizeof(label)); // 行號后應跟冒號
		in >> label;
		if(label[0] != ':') throw TMESSAGE("format error!\n");
		DOUBLE a; // 隨后是本行的各個數值
		for(size_t j=0; j<m.colnum; j++) {
			in >> a;
			m.set(i,j,a);
		}
		if(!in.good()) throw TMESSAGE("read file error!");
	}
	m.checksym();	// 檢查是否為對稱陣
	return in;
}

matrix& matrix::operator*=(DOUBLE a) // 矩陣數乘常數a，結果放在原矩陣
{
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,a*value(i,j));
	return (*this);
}

matrix matrix::operator*(DOUBLE a) // 矩陣數乘常數a，原矩陣內容不變，返回一新矩陣
{
	matrix m(rownum, colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		m.set(i,j,a*value(i,j));
	return m;
}

matrix matrix::operator+(matrix& m) // 矩陣相加，產生一新的矩陣并返回它
{
	if(rownum != m.rownum || colnum != m.colnum) // 對應行列必須相同
		throw TMESSAGE("can not do add of matrix\n");
	matrix mm(rownum, colnum); // 產生一同自己同形的矩陣
	DOUBLE a;
	for(size_t i=0; i<rownum; i++) // 求和
	for(size_t j=0; j<colnum; j++)
	{
		a = value(i,j)+m.value(i,j);
		mm.set(i,j,a);
	}
	mm.checksym();	// 檢查是否為對稱陣
	return mm;
}

matrix& matrix::operator+=(matrix &m) // 矩陣求和，自己內容改變為和
{
	DOUBLE a;
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
	{
		a = value(i,j)+m.value(i,j);
		set(i,j,a);
	}
	checksym();	// 檢查是否為對稱陣
	return (*this);
}

matrix matrix::operator+(DOUBLE a)	// 矩陣加常數，指每一元素加一固定的常數，產生
											// 新矩陣，原矩陣不變
{
	matrix m(rownum, colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		m.set(i,j,a+value(i,j));
	return m;
}

matrix& matrix::operator+=(DOUBLE a)	// 矩陣自身加常數，自身內容改變
{
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,a+value(i,j));
	return (*this);
}

matrix operator-(DOUBLE a, matrix& m) { // 常數減矩陣，產生新的矩陣
	return (-m)+a;
};


matrix matrix::operator-(matrix& m) // 矩陣相減，產生新的矩陣
{
	matrix mm(*this); // 產生一同自己同形的矩陣
	mm += (-m);	// 加上相應的負矩陣
	return mm;
}

matrix& matrix::operator-=(matrix& m) // 矩陣相減，結果修改原矩陣
{
	(*this) += (-m);
	return (*this);
}

matrix matrix::operator*(matrix& m) // 矩陣相乘，原矩陣內容不變，產生一新矩陣
{
	if(colnum != m.rownum) // 必須滿足相乘條件
		throw TMESSAGE("can not multiply!");
	matrix mm(rownum,m.colnum); // 計算并產生一合要求的矩陣放乘積
	DOUBLE a;
	for(size_t i=0; i<rownum; i++) // 計算乘積
	for(size_t j=0; j<m.colnum; j++){
		a = 0.0;
		for(size_t k=0; k<colnum; k++)
			a += value(i,k)*m.value(k,j);
		mm.set(i,j,a);
	}
	mm.checksym();	// 檢查是否為對稱陣
	return mm; // 返回乘積
}

matrix& matrix::operator*=(matrix& m) // 矩陣相乘，自己修改成積矩陣
{
	(*this) = (*this)*m;
	return (*this);
}

matrix matrix::t()  // 矩陣轉置，產生新的矩陣
{
	matrix mm(*this);
	mm.trans();
	return mm;
}

int matrix::isnear(matrix& m, double e) // 檢查二矩陣是否近似相等
{
	if(rownum != m.rownum || colnum != m.colnum) return 0;
	for(size_t i=0; i< rownum; i++)
	for(size_t j=0; j< colnum; j++)
		if(fabs(value(i,j)-m.value(i,j)) > e) return 0;
	return 1;
}

int matrix::isnearunit(double e) // 檢查矩陣是否近似為單位矩陣
{
	if(rownum != colnum) return 0;
	return isnear(unit(rownum), e);
}

matrix matrix::row(size_t r) // 提取第r行行向量
{
	matrix mm(1, colnum);
	for(int i=0; i< colnum; i++)
		mm.set(0, i, value(r,i));
	return mm;
}

matrix matrix::col(size_t c) // 提取第c列列向量
{
	matrix mm(rownum, 1);
	for(int i=0; i< rownum; i++)
		mm.set(i, value(i, c));
	return mm;
}

void matrix::swapr(size_t r1, size_t r2, size_t k) // 交換矩陣r1和r2兩行
{
	DOUBLE a;
	for(size_t i=k; i<colnum; i++) {
		a = value(r1, i);
		set(r1, i, value(r2, i));
		set(r2, i, a);
	}
}

void matrix::swapc(size_t c1, size_t c2, size_t k) // 交換c1和c2兩列
{
	DOUBLE a;
	for(size_t i=k; i<colnum; i++) {
		a = value(i, c1);
		set(i, c1, value(i, c2));
		set(i, c2, a);
	}
}

DOUBLE matrix::maxabs(size_t &r, size_t &c, size_t k) // 求第k行和第k列后的主元及位置
{
	DOUBLE a=0.0;
	for(size_t i=k;i<rownum;i++)
	for(size_t j=k;j<colnum;j++)
		if(a < fabs(value(i,j))) {
			r=i;c=j;a=fabs(value(i,j));
		}
	return a;
}

size_t matrix::zgsxy(matrix & m, int fn) // 進行主高斯消元運算，fn為參數，缺省為0
 /* 本矩陣其實是常數陣，而矩陣m必須是方陣