cmatrix.cpp

陈必红 编写的《用C++语言编写的数学常用算法》光盘内容

#include "stdio.h"
#include "string.h"
#include "complex.h"
#include "cmatrix.h"
#include "fstream.h"

// 缺省的构造函数，产生0行0列空复矩阵
cmatrix::cmatrix(cbuffer * b): rownum(0),colnum(0),
	isneg(0),istrans(0),isconj(0)
{
		if(b){
			buf=b;
			buf->alloc(0);
			}
		else
			buf = getnewcbuffer(0); // 如果缺省的b没有给出，则产生一
				// 长度为行乘列个复数和缓存
}
// 构造一r行c列的复矩阵，可选的复缓存指向b
cmatrix::cmatrix(size_t r, size_t c, cbuffer * b):
		rownum(r),colnum(c),istrans(0),isneg(0),isconj(0)
{
		if(b){
			buf=b;
			buf->alloc(r*c);
			}
		else
			buf = getnewcbuffer(r*c); // 如果缺省的b没有给出，则产生一
				// 长度为行乘列个复数和缓存
}

cmatrix::cmatrix(size_t n, cbuffer * b):
		rownum(n),colnum(1),istrans(0),isneg(0),isconj(0)
{
		if(b){
			buf=b;
			buf->alloc(n);
			}
		else
			buf = getnewcbuffer(n); // 如果缺省的b没有给出，则产生一
				// 长度为行乘列个复数和缓存
}

cmatrix cunit(size_t n) // 产生n阶单位矩阵
{
	cmatrix 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,i,1.0);
			else m.set(i,j,0.0);
	return m;
}

// 拷贝构造函数
cmatrix::cmatrix(cmatrix& m)
{
	rownum = m.rownum; // 复制行数和列数
	colnum = m.colnum;
	istrans = m.istrans; // 复制标志
	isneg = m.isneg;
	isconj = m.isconj;
	buf = m.buf; // 指向相同的缓存
	buf->refnum++;	// 缓存的引用数增1
}

// 用数据文件构造复矩阵
cmatrix::cmatrix(const char * filename, cbuffer * b) :
	istrans(0),isneg(0),isconj(0)
{
	char label[10];
	ifstream in(filename); // 打开文件输入
	in >> label; // 文件以cmatrix关键词打头
	if(strcmp(label, "cmatrix")!=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 = getnewcbuffer(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!");
		COMPLEX a;
		for(size_t j=0; j<colnum; j++) { // 随后是一行的内容
			in >> a;
			set(i,j,a);
		}
		if(!in.good()) throw TMESSAGE("read file error!");
	}
}

// 从内存数据块构造复矩阵。data指向内存数据块，构造完毕后data指向的内存可以释放
cmatrix::cmatrix(void * data, size_t r, size_t c, cbuffer * b):
		rownum(r),colnum(c),istrans(0),isneg(0),isconj(0)
{
	if(b){
		buf=b;
		buf->alloc(r*c);
		}
	else
		buf = getnewcbuffer(r*c);
	COMPLEX * d = (COMPLEX *)data;
	for(size_t i=0; i<r*c; i++)
		buf->set(i,d[i]);
}

cmatrix::cmatrix(matrix& m, cbuffer* b) : // 由实矩阵产生一个内容相同的复矩阵
	istrans(0),isneg(0),isconj(0)
{
	rownum = m.rownum;
	colnum = m.colnum;
	if(b) { buf=b;	// 产生所需长度的缓存
		buf->alloc(rownum*colnum);
	}
	else buf = getnewcbuffer(rownum*colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,m.value(i,j));
}

cmatrix::cmatrix(matrix& mr, matrix& mi, cbuffer * b): // 两个实矩阵成为实部和虚部
	istrans(0),isneg(0),isconj(0)
{
	if(mr.rownum != mi.rownum || mr.colnum != mi.colnum)
		throw TMESSAGE("dimension not same!");
	rownum = mi.rownum;
	colnum = mr.colnum;
	if(b) { buf=b;	// 产生所需长度的缓存
		buf->alloc(rownum*colnum);
	}
	else buf = getnewcbuffer(rownum*colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,COMPLEX(mr.value(i,j),mi.value(i,j)));
}

matrix cmatrix::real() // 由实部生成实矩阵
{
	matrix m(rownum, colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++) {
		m.set(i,j,::real(value(i,j)));
	}
	return m;
}

matrix cmatrix::image() // 由虚部生成实矩阵
{
	matrix m(rownum, colnum);
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		m.set(i,j,imag(value(i,j)));
	return m;
}

// 矩阵赋值
cmatrix& cmatrix::operator=(cmatrix& m)
{
	rownum = m.rownum; // 复制行数和列数
	colnum = m.colnum;
	if(buf == m.buf) // 如果本来就指向相同的缓存，则返回
		return (*this);
	buf->refnum--;	// 否则放弃原来的缓存，先将原缓存引用数减1
						// 如果不为0说明还有其它的矩阵引用数据
	if(!buf->refnum)delete buf; // 否则删除缓存
	istrans = m.istrans;	// 复制标志
	isneg = m.isneg;
	isconj = m.isconj;
	buf = m.buf; // 复制缓存指针
	buf->refnum++; // 缓存引用数增1
	return (*this);
}

cmatrix& cmatrix::operator=(matrix& m)
{
	rownum = m.rownum; // 复制行数和列数
	colnum = m.colnum;
	buf->refnum--;
	if(!buf->refnum) delete buf; // 放弃原来缓存
	istrans = isneg = isconj = 0; // 标志清零
	buf = getnewcbuffer(rownum*colnum); // 申请新的缓存
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,m.value(i,j));
	return (*this);
}

cmatrix& cmatrix::operator=(COMPLEX a) // 将矩阵所有元素设为常数a
{
	if(rownum == 0 || colnum==0 ) return (*this);
	istrans = isneg = isconj = 0; // 标志清零
	for(size_t i=0; i<rownum; i++)
	for(size_t j=0; j<colnum; j++)
		set(i,j,a);
	return (*this);
}

COMPLEX cmatrix::operator()(size_t r, size_t c)
 // 重载函数运算符，返回第r行c列的值
{
	if(r>= rownum || c >= colnum)
		throw TMESSAGE("range overflow");
	return value(r,c);
}

COMPLEX cmatrix::operator()(size_t r)
	// 重载函数运算符，返回第r行0列的值，用于向量
{
	return (*this)(r,0);
}


// 流输出
ostream& operator<<(ostream& o, cmatrix& m)
{	// 先输出关键词cmatrix,后跟行数和列数
	o << "cmatrix " << m.rownum << '\t' << m.colnum << endl;
	// 每行由行号，冒号，后跟一行的复数构成
	for(size_t i=0; i<m.rownum; i++) {
		o<< (i+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, cmatrix& m)
{
	char label[10];
	in.width(sizeof(label));
	in >> label;
	if(strcmp(label, "cmatrix")!=0) // 以关键词cmatrix打到
		throw TMESSAGE("format error!");
	in >> m.rownum >> m.colnum; // 后跟行数和列数
	if(!in.good()) throw TMESSAGE("read file error!");
	m.buf->refnum--; // 放弃原来的缓存
	if(!m.buf->refnum) delete m.buf;
	m.buf = getnewcbuffer(m.rownum*m.colnum); // 申请新的缓存
	m.istrans = m.isneg = m.isconj = 0; // 标志清零
	size_t line;
	for(size_t i=0; i<m.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!"); // 后跟冒号
		COMPLEX 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!");
	}
	return in;
}

// 数乘矩阵，原矩阵内容改变为结果
cmatrix& cmatrix::operator*=(COMPLEX 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);
}

// 数乘矩阵，原矩阵内容不变，返回结果矩阵
cmatrix cmatrix::operator*(COMPLEX a)
{
	cmatrix 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;
}

cmatrix cmatrix::operator+(COMPLEX a) // 矩阵加常数，产生新矩阵
{
	cmatrix 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;

}

cmatrix& cmatrix::operator+=(COMPLEX 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);
}

// 矩阵相加
cmatrix cmatrix::operator+(cmatrix& m)
{
	if(rownum != m.rownum || colnum != m.colnum) // 行列必须对应相等
		throw TMESSAGE("can not do add of cmatrix");
	cmatrix mm(rownum, colnum);
	COMPLEX 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);
	}
	return mm;
}

// 矩阵相加，不产生新的矩阵，修改原矩阵的内容
cmatrix& cmatrix::operator+=(cmatrix &m)
{
	COMPLEX 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);
	}
	return (*this);
}

cmatrix cmatrix::operator-(cmatrix& a)	// 矩阵相减
{
	return (*this)+(-a);
}

cmatrix& cmatrix::operator-=(cmatrix& a) // 矩阵自身减矩阵
{
	(*this)+=(-a);
	return (*this);
}

// 矩阵相乘，返回新产生的结果矩阵
cmatrix cmatrix::operator*(cmatrix& m)
{
	if(colnum != m.rownum) // 前矩阵的列数必须等于后矩阵的行数
		throw TMESSAGE("can not multiply!");
	cmatrix mm(rownum,m.colnum);
	COMPLEX 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);
	}
	return mm;
}

// 矩阵相乘，改变原矩阵的内容
cmatrix& cmatrix::operator*=(cmatrix& m)
{
	(*this) = (*this)*m;
	return (*this);
}

cmatrix cmatrix::operator-() // 矩阵求负，产生新的矩阵
{
	cmatrix m(*this);
	m.neg();
	return m;
}