#ifndef CMATRIX_H

#define CMATRIX_H

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

#include "cbuffer.h"
#include "matrix.h"

// 定義復矩陣類
class cmatrix {
 public:
	cbuffer * buf;	// 指向復矩陣緩存類的指針
	size_t rownum, colnum;	// 矩陣的行數與列數
	unsigned istrans:1;	// 轉置標記
	unsigned isneg:1;		// 取負標記
	unsigned isconj:1;	// 共軛標記
	cmatrix(cbuffer * b=0); // 缺省構造函數，產生0行0列空矩陣
	cmatrix(size_t n, cbuffer * b=0); // 產生n維列向量
	cmatrix(size_t r, size_t c, cbuffer * b=0); // 構造函數，構造一個空的復矩陣
	cmatrix(cmatrix& m);	// 拷貝構造函數
	cmatrix(const char * filename, cbuffer * b=0); // 從數據文件構造一個復矩陣
	cmatrix(void * data, size_t r, size_t c=1, cbuffer * b=0);
		 // 從內存數據產生復矩陣
	cmatrix(matrix& m, cbuffer * b=0); // 由實矩陣產生一個內容相同的復矩陣
	cmatrix(matrix& mr, matrix& mi, cbuffer * b=0); // 兩個實矩陣成為實部和虛部
												// 構成一復矩陣
	matrix real();	// 由矩陣的實部生成實矩陣
	matrix image(); // 由矩陣的虛部生成實矩陣
	COMPLEX operator()(size_t r, size_t c); // 重載函數運算符，返回第r行c列的值
	COMPLEX operator()(size_t r);	// 重載函數運算符，返回第r行0列的值，用于向量
	virtual void set(size_t r, size_t c, COMPLEX v) { // 將第r行c列的值設為v
		size_t l;
		v = isneg ? -v : v;
		v = isconj ? ::conj(v) : v;
		l = istrans ? r+c*rownum : r*colnum+c;
		buf = buf->set(l, v); };
	void set(size_t r, COMPLEX v) {
		set(r,0,v);};	// 設置向量中第r個元素的值
	virtual ~cmatrix(){ // 析構函數
		buf->refnum--;	// 緩存的引用數減1
		if(!buf->refnum) delete buf;}; // 如果緩存的引用數為0，則釋放緩存
	cmatrix& operator=(cmatrix& m);	// 重載賦值運算
	cmatrix& operator=(matrix& m); // 將實矩陣賦給復矩陣
	cmatrix& operator=(COMPLEX a); // 將矩陣所有元素設為常數a
	cmatrix operator*(COMPLEX a);	// 數乘矩陣，產生新的矩陣
	cmatrix& operator*=(COMPLEX a); // 數乘矩陣，只是改動原矩陣
	friend cmatrix operator*(COMPLEX a, cmatrix& m); // 數乘矩陣
	cmatrix operator+(COMPLEX a); // 矩陣加常數，產生新矩陣
	cmatrix& operator+=(COMPLEX a); // 矩陣加常數，更動原矩陣
	friend cmatrix operator+(COMPLEX a, cmatrix& m); // 矩陣加常數，常數在前
	cmatrix operator+(cmatrix& m); // 矩陣加法，產生新的矩陣
	cmatrix& operator+=(cmatrix &m); // 矩陣加法，結果取代原矩陣

	cmatrix operator*(cmatrix& m);	// 矩陣乘法，產生并返回新的矩陣
	cmatrix& operator*=(cmatrix& m);	// 矩陣乘法，只是更動原矩陣

	cmatrix operator-(); // 矩陣求負，產生新的矩陣
	cmatrix& neg();	// 將自己求負，不產生新的矩陣

	cmatrix	t();		// 矩陣轉置，產生新的矩陣
	cmatrix& trans()	// 矩陣自身轉置
		{  size_t r = rownum; rownum = colnum; colnum = r; // 交換行數和列數
		istrans = !istrans; return (*this);};
	cmatrix&	conj()	// 矩陣共軛，原矩陣更動
		{isconj = !isconj; return (*this);};
	cmatrix	operator!();		// 矩陣共軛，產生新的矩陣
	cmatrix& transconj()		// 矩陣自身共軛轉置
		{ isconj = !isconj; trans(); return (*this);};
	cmatrix	tc();		// 矩陣共軛轉置，產生新的矩陣

	cmatrix operator-(cmatrix& a);	// 矩陣相減
	cmatrix& operator-=(cmatrix& a); // 矩陣自身減矩陣
	cmatrix& operator-=(COMPLEX a); // 矩陣自身減常數
	cmatrix operator-(COMPLEX a);	// 矩陣減常數，產生新的矩陣
	friend cmatrix operator-(COMPLEX a, cmatrix m); // 常數減矩陣，產生新的矩陣

	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){ // 交換兩個元素
		COMPLEX 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(cmatrix & m, int fn=0);	// 主高斯消元運算
	cmatrix& operator/=(cmatrix m); // 用主高斯消元法求解線性方程的解
		// 矩陣本身為系數矩陣，m為常數向量，返回解向量

	cmatrix& inv();		// 用全選主元高斯-約當法求逆矩陣
	cmatrix operator~();	// 求逆矩陣，但產生新矩陣
	friend cmatrix operator/(COMPLEX a, cmatrix& m); // 求逆矩陣再乘常數
	cmatrix& operator/=(COMPLEX a); // 所有元素乘a的倒數，自身改變
	cmatrix operator/(COMPLEX a); // 所有元素乘a的倒數，產生新的矩陣

	cmatrix& fft(int l=0);


	friend ostream& operator<<(ostream& o, cmatrix& m); // 矩陣的流輸出
	friend istream& operator>>(istream& in, cmatrix& m); // 矩陣的流輸入

	virtual COMPLEX value(size_t r, size_t c) { // 返回第r行c列的復數值
		COMPLEX v;
		v = istrans ? (*buf)[r+c*rownum] : (*buf)[r*colnum+c];
		v = isneg ? -v : v;
		v = isconj ? ::conj(v) : v;
		return v; };
};

inline cmatrix operator*(COMPLEX a, cmatrix& m) {
	return m*a;
};

inline cmatrix operator+(COMPLEX a, cmatrix& m){ // 矩陣加常數，常數在前
	return m+a;
};

cmatrix operator/(COMPLEX a, cmatrix& m); // 求逆矩陣再乘常數


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

cmatrix cunit(size_t n);	// 產生n階單位復矩陣

cmatrix fft(matrix& m);	// 對實向量作fft變換
matrix convolute(matrix&a, matrix& b);	// 求矩陣a與b的離散卷積

#endif // CMATRIX_H