cmatrix.cpp

用c++编写数学常用算法的源代码

cmatrix& cmatrix::neg() // 矩阵自身求负
{
	isneg = !isneg;	// 取负标志求反
	return (*this);
}

cmatrix cmatrix::t()		// 矩阵转置，产生新的矩阵
{
	cmatrix m(*this);
	m.trans();
	return m;
}

cmatrix cmatrix::operator!()		// 矩阵共轭，产生新的矩阵
{
	cmatrix m(*this);
	return m.conj();
}

cmatrix cmatrix::tc()		// 矩阵共轭转置，产生新的矩阵
{
	cmatrix m(*this);
	return m.transconj();
};

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

cmatrix cmatrix::operator-(COMPLEX a)	// 矩阵减常数，产生新的矩阵
{
	return (*this)+(-a);
}

cmatrix operator-(COMPLEX a, cmatrix m) // 常数减矩阵，产生新的矩阵
{
	return (-m)+a;
}

void cmatrix::swapr(size_t r1, size_t r2, size_t k) // 交换矩阵r1和r2两行
{
	COMPLEX a;
	for(size_t i=k; i<colnum; i++) {
		a = value(r1, i);
		set(r1, i, value(r2, i));
		set(r2, i, a);
	}
}

void cmatrix::swapc(size_t c1, size_t c2, size_t k) // 交换c1和c2两列
{
	COMPLEX a;
	for(size_t i=k; i<colnum; i++) {
		a = value(i, c1);
		set(i, c1, value(i, c2));
		set(i, c2, a);
	}
}

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

size_t cmatrix::zgsxy(cmatrix & m, int fn) // 进行主高斯消元运算，fn为参数，缺省为0
 /* 本矩阵其实是常数阵，而矩阵m必须是方阵
	运算过程其实是对本矩阵和m同时作行初等变换，
	运算结果m的对角线相乘将是行列式，而本矩阵则变换成
	自己的原矩阵被m的逆阵左乘，m的秩被返回，如果秩等于阶数
	则本矩阵中的内容已经是唯一解
 */
{
	if(rownum != m.rownum || m.rownum != m.colnum) // 本矩阵行数必须与m相等
						// 且m必须是方阵
		throw TMESSAGE("can not divide!");
	lbuffer * bb = getnewlbuffer(rownum); // 产生一维数为行数的长整数缓存区
	lbuffer & b = (*bb); // 用引用的办法使下面的程序容易懂
	size_t is;
	COMPLEX a;
	size_t i,j,rank=0;
	for(size_t k=0; k<rownum; k++) { // 从第0行到第k行进行主高斯消元
		if(m.maxabs(is, i, k)==0) // 求m中第k级主元，主元所在的行，列存在is,i中
				// 返回的主元值放在a中
			break; // 如果主元为零，则m不可逆，运算失败
		rank = k+1; // rank存放当前的阶数
		b.retrieve(k) = i;  // 将长整数缓存区的第k个值设为i
		if(i != k)
			m.swapc(k, i); // 交换m中i,k两列
		if(is != k) {
			m.swapr(k, is, k); // 交换m中i,k两行,从k列以后交换
			swapr(k, is); // 交换本矩阵中i,k两行
		}
		a = m.value(k,k); // 取出主元值
		for (j=k+1;j<rownum;j++) // 本意是将m的第k行除以主元
				// 但只需把第k行的第k+1列以上除以主元即可
				// 这样还保留了主元作行列式运算用
			m.set(k,j,m.value(k,j)/a);
		for (j=0;j<colnum;j++) // 将本矩阵的第k行除以主元
			set(k,j,value(k,j)/a);
		// 上面两步相当于将m和本矩阵构成的增广矩阵第k行除以主元
		// 下面对增广矩阵作行基本初等变换使第k行的其余列均为零
		// 但0值无必要计算，因此从第k+1列开始计算
		for(j=k+1;j<rownum;j++) // j代表列，本矩阵的行数就是m的列数
		for(i=0;i<rownum;i++) //i代表行，依次对各行计算，k行除外
			if(i!=k)
				m.set(i,j,m.value(i,j)-m.value(i,k)*m.value(k,j));
		// 对本矩阵亦作同样的计算
		for(j=0;j<colnum;j++)
		for(i=0;i<rownum;i++)
			if(i!=k)
				set(i,j,value(i,j)-m.value(i,k)*value(k,j));
	} // 主高斯消元循环k结束
	if(fn == 1) {
		for(j=0; j<rank; j++)
		for(i=0; i<rownum; i++)
			if(i==j) m.set(i,i,1.0);
			else
				m.set(i,j,0.0);
		for(k = rank; k>0; k--)
			m.swapc(k-1,(size_t)b[k-1]);
	}
	for(k = rank; k>0; k--) // 将本矩阵中的各行按b中内容进行调节
		if(b[k-1] != k-1)
			swapr(k-1,(size_t)b[k-1]); // 行交换
	delete bb; // 释放长整数缓存
	return rank; // 返回mm的秩
}

cmatrix& cmatrix::operator/=(cmatrix m) // 利用重载的除法符号/=来解方程
	// 本矩阵设为d,则方程为mx=d,考虑解写成x=d/m的形式，
	// 而方程的解也存放在d中，则实际编程时写d/=m
{
	if(zgsxy(m)!=rownum) // 如秩不等于m的阶数，则方程无解
		throw TMESSAGE("can not divide!");
	return *this;
}

cmatrix& cmatrix::fft(int l)
{
	if(colnum != 1) throw TMESSAGE("fft need vector!");
	size_t n,i;
	int k;
	n = 1;
	for(k=0;k<10000;k++) {
		if(n>=rownum) break;
		n*=2;
	}
	if(rownum < n) {	// 将数组扩大为2的乘幂维
		cmatrix m(n);
		for(i=0; i<rownum; i++)
			m.set(i,value(i,0));
		for(i=rownum; i<n; i++)
			m.set(i,0.0);
		(*this)=m;
	}
	size_t it,m,is,j,nv;
	int l0;
	COMPLEX q,s,v,podd;
	DOUBLE p;
	cmatrix f(n);
	for (it=0; it<n; it++)
		{ m=it; is=0;
		  for (i=0; i<k; i++)
			 { j=m/2; is=2*is+(m-2*j); m=j;}
		  f.set(it,value(is,0));
		}
	set(0,1.0);
	cbuffer & pp = *buf;
	cbuffer & ff = *(f.buf);
	p = 2.0*M_PI/n;
	pp[1] = COMPLEX(cos(p),-sin(p));
	if (l!=0) pp[1] = ::conj(pp[1]);
	for (i=2; i<n; i++)
		pp[i] = pp[1]*pp[i-1];
	for (it=0; it<n-1; it+=2)
		{  v = ff[it];
			ff[it] += ff[it+1];
			ff[it+1] = v - ff[it+1];
		}
	m=n/2; nv=2;
	for (l0=k-2; l0>=0; l0--)
		{ m/=2; nv*=2;
		  for (it=0; it<=(m-1)*nv; it+=nv)
			 for (j=0; j<=(nv/2)-1; j++)
				{
					podd = pp[m*j]*ff[it+j+nv/2];
					ff[it+j+nv/2] = ff[it+j]-podd;
					ff[it+j] += podd;
				}
		}
	 if (l!=0)
		for (i=0; i<n; i++)
		  ff[i]/=n;
	(*this) = f;
	return (*this);
}

cmatrix fft(matrix& m)
{
	cmatrix mm(m);
	return mm.fft();
}

matrix convolute(matrix&a, matrix& b)	// 求矩阵a与b的离散卷积
{
	if(a.colnum != 1 || b.colnum != 1)throw TMESSAGE("must be a vector");
	size_t nn = a.rownum+b.rownum;
	size_t n;
	int k;
	n = 1;
	for(k=0;k<10000;k++) {
		if(n>=nn) break;
		n*=2;
	}
	matrix aa(n),bb(n);
	size_t i;
	aa = 0.0;
	bb = 0.0;
	for(i=0; i<a.rownum; i++)
		aa.set(i,a(i));
	for(i=0; i<b.rownum; i++)
		bb.set(i,b(i));
	cmatrix aafft(fft(aa)),bbfft(fft(bb));
	for(i=0; i<n; i++)
		aafft.set(i,aafft(i)*bbfft(i));
	nn--;
	matrix c(nn),d;
	d = (aafft.fft(1)).real();
	for(i=0;i<nn;i++)
		c.set(i,0,d(i));
	return c;
}

cmatrix& cmatrix::inv()		// 用全选主元高斯-约当法求逆矩阵
{
	if(rownum != colnum) throw TMESSAGE("must be a square matrix");
	size_t n = rownum;
	lbuffer * isp = getnewlbuffer(n);
	lbuffer * jsp = getnewlbuffer(n);
	lbuffer& is = (*isp);
	lbuffer& js = (*jsp);
	size_t i,j,k,u,v;
	DOUBLE d;
	for (k=0; k<n; k++)
		{
			if((d=maxabs(u,v,k))== 0.0) {
				delete isp;
				delete jsp;
				throw TMESSAGE("can not inverse");
			}
			is[k]=u;
			js[k]=v;
/*
		 d=0.0;
		  for (i=k; i<=n-1; i++)
		  for (j=k; j<=n-1; j++)
			 { u=i*n+j;
				p=ar[u]*ar[u]+ai[u]*ai[u];
				if (p>d) { d=p; is[k]=i; js[k]=j;}
			 }
		  if (d+1.0==1.0)
			 { free(is); free(js); printf("err**not inv\n");
				return(0);
			 }
			 */
		  if (is[k]!=k) swapr(k,(size_t)is[k]);
/*
			 for (j=0; j<=n-1; j++)
				{ u=k*n+j; v=is[k]*n+j;
				  t=ar[u]; ar[u]=ar[v]; ar[v]=t;
				  t=ai[u]; ai[u]=ai[v]; ai[v]=t;
				} */
		  if (js[k]!=k) swapc(k,(size_t)js[k]);
/*			 for (i=0; i<=n-1; i++)
				{ u=i*n+k; v=i*n+js[k];
				  t=ar[u]; ar[u]=ar[v]; ar[v]=t;
				  t=ai[u]; ai[u]=ai[v]; ai[v]=t;
				} */
//		  l=k*n+k;
		  set(k,k,::conj(value(k,k))/d);
//		  ar[l]=ar[l]/d; ai[l]=-ai[l]/d;
		  for (j=0; j<n; j++)
			 if (j!=k)
				set(k,j,value(k,j)*value(k,k));
/*
				{ u=k*n+j;
				  p=ar[u]*ar[l]; q=ai[u]*ai[l];
				  s=(ar[u]+ai[u])*(ar[l]+ai[l]);
				  ar[u]=p-q; ai[u]=s-p-q;
				} */
		  for (i=0; i<n; i++)
			 if (i!=k)
				{ // v=i*n+k;
				  for (j=0; j<n; j++)
					 if (j!=k)
						set(i,j,value(i,j)-value(k,j)*value(i,k));
/*
						{ u=k*n+j;  w=i*n+j;
						  p=ar[u]*ar[v]; q=ai[u]*ai[v];
						  s=(ar[u]+ai[u])*(ar[v]+ai[v]);
						  t=p-q; b=s-p-q;
						  ar[w]=ar[w]-t;
						  ai[w]=ai[w]-b;
						} */
				}
		  for (i=0; i<n; i++)
			 if (i!=k)
				set(i,k,-value(i,k)*value(k,k));
/*				{ u=i*n+k;
				  p=ar[u]*ar[l]; q=ai[u]*ai[l];
				  s=(ar[u]+ai[u])*(ar[l]+ai[l]);
				  ar[u]=q-p; ai[u]=p+q-s;
				} */
		}
	 for (k=n; k>0; k--)
		{ if (js[k-1]!=k-1) swapr((size_t)js[k-1],k-1);
/*			 for (j=0; j<n; j++) swapr((size_t)js[k-1],k-1);
				{ u=k*n+j; v=js[k]*n+j;
				  t=ar[u]; ar[u]=ar[v]; ar[v]=t;
				  t=ai[u]; ai[u]=ai[v]; ai[v]=t;
				} */
		  if (is[k-1]!=k-1) swapc((size_t)is[k-1],k-1);
/*			 for (i=0; i<=n-1; i++)
				{ u=i*n+k; v=i*n+is[k];
				  t=ar[u]; ar[u]=ar[v]; ar[v]=t;
				  t=ai[u]; ai[u]=ai[v]; ai[v]=t;
				} */
		}
	delete isp;
	delete jsp;
	return (*this);
}

cmatrix cmatrix::operator~()	// 求逆矩阵，但产生新矩阵
{
	cmatrix m(*this);
	m.inv();
	return m;
}

cmatrix operator/(COMPLEX a, cmatrix m) // 求逆矩阵再乘常数
{
	return (~m)*a;
}

cmatrix& cmatrix::operator/=(COMPLEX a) // 所有元素乘a的倒数，自身改变
{
	(*this)*=1.0/a;
	return (*this);
}

cmatrix cmatrix::operator/(COMPLEX a) // 所有元素乘a的倒数，产生新的矩阵
{
	cmatrix m(*this);
	m/=a;
	return m;
}