matrix.cpp

快速fft变换的c实现

      return false;
   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 if (_m->Val[i][j] != _m->Val[j][i])
	    return false;
   return true;
}
	   
// Determine if the matrix is skew-symmetric
 bool CMatrix::IsSkewSymmetric () 
{
   if (_m->Row != _m->Col)
      return false;
   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 if (_m->Val[i][j] != -_m->Val[j][i])
	    return false;
   return true;
}
   
// Determine if the matrix is upper triangular
 bool CMatrix::IsUpperTriangular () 
{
   if (_m->Row != _m->Col)
      return false;
   for (size_t i=1; i < _m->Row; i++)
      for (size_t j=0; j < i-1; j++)
	 if (_m->Val[i][j] != double(0))
	    return false;
   return true;
}

// Determine if the matrix is lower triangular
bool CMatrix::IsLowerTriangular () 
{
   if (_m->Row != _m->Col)
      return false;

   for (size_t j=1; j < _m->Col; j++)
      for (size_t i=0; i < j-1; i++)
	 if (_m->Val[i][j] != double(0))
	    return false;
   return true;
}

CMatrix CMatrix::Inversion ()
{
	CMatrix *temp = new CMatrix(_m->Col, _m->Row);
	for (size_t i=0; i<_m->Col; i++)
	{
		for (size_t j=0; j<_m->Row; j++)
		{
			(*temp)(i, j) = _m->Val[j][i];
		}
	}
	return *temp;
}

//用豪斯荷尔德(Householder)变换及变形QR算法对一般实矩阵进行奇异值分解
//设A为m*n实矩阵，则存在一个m*m的列正交矩阵U和n*n的列正交矩阵V，使A=U*S*V'
//A : 双精度实型二维数组,体积为m*n
//U : 双精度实型二维数组,体积为m*m,返回时存放左奇异向量U
//S : 双精度实型二维数组,体积为m*n,返回时其对角线给出奇异值(以降次序排列),其余元素均为0
//V : 双精度实型二维数组,体积为n*n,返回时存放右奇异向量V (注意:已转置)   A=U*S*V'
//dbEps : 双精度实型变量,给定的精度要求
//nIt : 整型变量,指明迭代的最大次数,缺省值为60
//返回: 若工作正常，返回TRUE；若迭代了nIt次还未求得某个奇异值的情况，返回FALSE
BOOL CMatrix::USVtDecompose(const CMatrix &A, CMatrix &U, CMatrix &S, CMatrix &V, double dbEps, int nIt/*=60*/)
{
	int i,j,k,l,ll,kk,ix,iy,mm,nn,iz,m1,ks,m,n,ka;
	double d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];
	double *s,*e,*w;

	//initial:
	m=A.RowNo();	//行数
	n=A.ColNo();	//列数

	if (m>=n) ka=m+1;
	else ka=n+1;	//ka=max(m,n)+1

	s=new double[ka];
	e=new double[ka];
	w=new double[ka];

//	U.Initialization(m,m);
	U.SetSize(m,m);
	S=A;
//	V.Initialization(n,n);
	V.SetSize(n,n);
	k=n;
	if (m-1<n) k=m-1;
	l=m;
	if (n-2<m) l=n-2;
	if (l<0) l=0;
	ll=k;
	if (l>k) ll=l;
	if (ll>=1)
	{
		for (kk=1;kk<=ll;kk++)		//kk相当于列
		{
			if (kk<=k)
			{
				d=0.0;
				for (i=kk;i<=m;i++)  //i相当于行
				{
					ix=(i-1)*n+kk-1;	//ix相当于i行kk列元素位置
					d+=S(i-1,kk-1)*S(i-1,kk-1);
					//d=d+S.pData[ix]*S.pData[ix];
				}
				s[kk-1]=sqrt(d);
				if (s[kk-1]!=0.0)
				{
					ix=(kk-1)*n+kk-1;
					if ( S(kk-1,kk-1) !=0.0       )// S.pData[ix]!=0.0)
					{
						s[kk-1]=fabs(s[kk-1]);
					//	if (S.pData[ix]<0.0) 
						if(S(kk-1,kk-1)<0.0)		
							s[kk-1]=-s[kk-1];
					}
					for (i=kk;i<=m;i++)
					{
						iy=(i-1)*n+kk-1;
						S(i-1,kk-1)=S(i-1,kk-1)/s[kk-1];
						//S.pData[iy]=S.pData[iy]/s[kk-1];
					}

					//S.pData[ix]=1.0+S.pData[ix];
					S(kk-1,kk-1)+=1.0;
				}
				s[kk-1]=-s[kk-1];
			}

			if (n>=kk+1)
			{
				for (j=kk+1;j<=n;j++)
				{
					if ((kk<=k)&&(s[kk-1]!=0.0))
					{
						d=0.0;
						for (i=kk;i<=m;i++)
						{
							ix=(i-1)*n+kk-1;
							iy=(i-1)*n+j-1;
							d+=S(i-1,kk-1)*S(i-1,j-1);
							//d=d+S.pData[ix]*S.pData[iy];
						}

					//	d=-d/S.pData[(kk-1)*n+kk-1];
						d=-d/S(kk-1,kk-1);
						for (i=kk;i<=m;i++)
						{
							ix=(i-1)*n+j-1;
							iy=(i-1)*n+kk-1;
							S(i-1,j-1)+=	d*S(i-1,kk-1);
							//S.pData[ix]=S.pData[ix]+d*S.pData[iy];
						}
					}
					//e[j-1]=S.pData[(kk-1)*n+j-1];
					e[j-1]=S(kk-1,j-1);
				}
			}
			if (kk<=k)
			{
				for (i=kk;i<=m;i++)
				{
					ix=(i-1)*m+kk-1;
					iy=(i-1)*n+kk-1;
					//U.pData[ix]=S.pData[iy];
					U(i-1,kk-1)=S(i-1,kk-1);
				}
			}
			if (kk<=l)
			{
				d=0.0;
				for (i=kk+1;i<=n;i++)
					d=d+e[i-1]*e[i-1];

				e[kk-1]=sqrt(d);
				if (e[kk-1]!=0.0)
				{
					if (e[kk]!=0.0)
					{
						e[kk-1]=fabs(e[kk-1]);
						if (e[kk]<0.0)
							e[kk-1]=-e[kk-1];
					}
					for (i=kk+1;i<=n;i++)
						e[i-1]=e[i-1]/e[kk-1];

					e[kk]=1.0+e[kk];
				}
				e[kk-1]=-e[kk-1];
				if ((kk+1<=m)&&(e[kk-1]!=0.0))
				{
					for (i=kk+1;i<=m;i++)
						w[i-1]=0.0;

					for (j=kk+1;j<=n;j++)
						for (i=kk+1;i<=m;i++)

							//w[i-1]=w[i-1]+e[j-1]*S.pData[(i-1)*n+j-1];
							w[i-1]+=e[j-1]*S(i-1,j-1);

					for (j=kk+1;j<=n;j++)
						for (i=kk+1;i<=m;i++)
						{
							ix=(i-1)*n+j-1;
							S(i-1,j-1)-=w[i-1]*e[j-1]/e[kk];
							//S.pData[ix]=S.pData[ix]-w[i-1]*e[j-1]/e[kk];
						}
				}

				for (i=kk+1;i<=n;i++)
					//V.pData[(i-1)*n+kk-1]=e[i-1];
					V(i-1,kk-1)=e[i-1];
			}
		}
	}

	mm=n;
	if (m+1<n) mm=m+1;
	if (k<n) //s[k]=S.pData[k*n+k];
		s[k]=S(k,k);

	if (m<mm) s[mm-1]=0.0;
	if (l+1<mm)  e[l]=S(l,mm-1);
		//e[l]=S.pData[l*n+mm-1];

	e[mm-1]=0.0;
	nn=m;
	if (m>n) nn=n;
	if (nn>=k+1)
	{
		for (j=k+1;j<=nn;j++)
		{
			for (i=1;i<=m;i++)
					U(i-1,j-1)=0.0;

				//U.pData[(i-1)*m+j-1]=0.0;
			//U.pData[(j-1)*m+j-1]=1.0;
			U(j-1,j-1)=1.0;
		}
	}
	if (k>=1)
	{
		for (ll=1;ll<=k;ll++)
		{
			kk=k-ll+1;
			iz=(kk-1)*m+kk-1;
			if (s[kk-1]!=0.0)
			{
				if (nn>=kk+1)
				{
					for (j=kk+1;j<=nn;j++)
					{
						d=0.0;
						for (i=kk;i<=m;i++)
						{
							ix=(i-1)*m+kk-1;
							iy=(i-1)*m+j-1;
							d+=U(i-1,kk-1)*U(i-1,j-1)/U(kk-1,kk-1);
							//d=d+U.pData[ix]*U.pData[iy]/U.pData[iz];
						}
						d=-d;
						for (i=kk;i<=m;i++)
						{
							ix=(i-1)*m+j-1;
							iy=(i-1)*m+kk-1;
							U(i-1,j-1)+=d*U(i-1,kk-1);
							//U.pData[ix]=U.pData[ix]+d*U.pData[iy];
						}
					}
				}
				for (i=kk;i<=m;i++)
				{
					ix=(i-1)*m+kk-1;
					U(i-1,kk-1)=-U(i-1,kk-1);
					//U.pData[ix]=-U.pData[ix];
				}
//				U.pData[iz]=1.0+U.pData[iz];
				U(kk-1,kk-1)+=1.0;
				if (kk-1>=1)
					for (i=1;i<=kk-1;i++)
						//U.pData[(i-1)*m+kk-1]=0.0;
						U(i-1,kk-1)=0.0;
			}
			else
			{
				for (i=1;i<=m;i++)
					//U.pData[(i-1)*m+kk-1]=0.0;
					U(i-1,kk-1)=0.0;
				//U.pData[(kk-1)*m+kk-1]=1.0;
				U(kk-1,kk-1)=1.0;
			}
		}
	}

	for (ll=1;ll<=n;ll++)
	{
		kk=n-ll+1; iz=kk*n+kk-1;
		if ((kk<=l)&&(e[kk-1]!=0.0))
		{
			for (j=kk+1;j<=n;j++)
			{
				d=0.0;
				for (i=kk+1;i<=n;i++)
				{
					ix=(i-1)*n+kk-1; iy=(i-1)*n+j-1;
				//	d=d+V.pData[ix]*V.pData[iy]/V.pData[iz];
					d+=V(i-1,kk-1)*V(i-1,j-1)/V(kk,kk-1);
				}
				d=-d;
				for (i=kk+1;i<=n;i++)
				{
					ix=(i-1)*n+j-1; iy=(i-1)*n+kk-1;
					//V.pData[ix]=V.pData[ix]+d*V.pData[iy];
					V(i-1,j-1)+=d*V(i-1,kk-1);
				}
			}
		}

		for (i=1;i<=n;i++)
			//V.pData[(i-1)*n+kk-1]=0.0;
			V(i-1,kk-1)=0.0;

	//	V.pData[iz-n]=1.0;
		V(kk-1,kk-1)=1.0;
	}

	for (i=1;i<=m;i++)
		for (j=1;j<=n;j++)
			//S.pData[(i-1)*n+j-1]=0.0;
			S(i-1,j-1)=0.0;

	m1=mm;
	while (1==1)	//
	{
		if (mm==0)
		{
			USVtDecompose_ppp(S,e,s,V,m,n);
			V=~V;	//转置
			delete [] s;
			delete [] e;
			delete [] w;
			return TRUE;
		}
		if (nIt==0)
		{
			USVtDecompose_ppp(S,e,s,V,m,n);
			V=~V;	//转置
			delete [] s;
			delete [] e;
			delete [] w;
			return FALSE;
		}
		kk=mm-1;
		while ((kk!=0)&&(fabs(e[kk-1])!=0.0))
		{
			d=fabs(s[kk-1])+fabs(s[kk]);
			dd=fabs(e[kk-1]);
			if (dd>dbEps*d) 
				kk=kk-1;
			else
				e[kk-1]=0.0;
		}

		if (kk==mm-1)
		{
			kk++;
			if (s[kk-1]<0.0)
			{
				s[kk-1]=-s[kk-1];
				for (i=1;i<=n;i++)
				{
					ix=(i-1)*n+kk-1;
					//V.pData[ix]=-V.pData[ix];
					V(i-1,kk-1)=-V(i-1,kk-1);
				}
			}
			while ((kk!=m1)&&(s[kk-1]<s[kk]))
			{
				d=s[kk-1]; s[kk-1]=s[kk]; s[kk]=d;
				if (kk<n)
					for (i=1;i<=n;i++)
					{
						ix=(i-1)*n+kk-1; iy=(i-1)*n+kk;
						//d=V.pData[ix]; V.pData[ix]=V.pData[iy]; V.pData[iy]=d;
						d=V(i-1,kk-1);
						V(i-1,kk-1)=V(i-1,kk);
						V(i-1,kk)=d;
					}
				if (kk<m)
					for (i=1;i<=m;i++)
					{
						ix=(i-1)*m+kk-1; iy=(i-1)*m+kk;
						//d=U.pData[ix]; U.pData[ix]=U.pData[iy]; U.pData[iy]=d;
						d=U(i-1,kk-1);
						U(i-1,kk-1)=U(i-1,kk);
						U(i-1,kk)=d;
					}
					kk++;
			}
			//it=60;
			mm--;
		}
		else
		{
			ks=mm;
			while ((ks>kk)&&(fabs(s[ks-1])!=0.0))
			{
				d=0.0;
				if (ks!=mm) d=d+fabs(e[ks-1]);
				if (ks!=kk+1) d=d+fabs(e[ks-2]);
				dd=fabs(s[ks-1]);
				if (dd>dbEps*d) 
					ks=ks-1;
				else
					s[ks-1]=0.0;
			}
			if (ks==kk)
			{
				kk++;
				d=fabs(s[mm-1]); 
				t=fabs(s[mm-2]);
				if (t>d) d=t;
				t=fabs(e[mm-2]);
				if (t>d) d=t;
				t=fabs(s[kk-1]);
				if (t>d) d=t;
				t=fabs(e[kk-1]);
				if (t>d) d=t;
				sm=s[mm-1]/d; sm1=s[mm-2]/d;
				em1=e[mm-2]/d;
				sk=s[kk-1]/d; ek=e[kk-1]/d;
				b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0;
				c=sm*em1; c=c*c; shh=0.0;
				if ((b!=0.0)||(c!=0.0))
				{
					shh=sqrt(b*b+c);
					if (b<0.0) shh=-shh;
					shh=c/(b+shh);
				}
				fg[0]=(sk+sm)*(sk-sm)-shh;
				fg[1]=sk*ek;
				for (i=kk;i<=mm-1;i++)
				{
					USVtDecompose_sss(fg,cs);
					if (i!=kk) e[i-2]=fg[0];
					fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1];
					e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1];
					fg[1]=cs[1]*s[i];
					s[i]=cs[0]*s[i];
					if ((cs[0]!=1.0)||(cs[1]!=0.0))
						for (j=1;j<=n;j++)
						{
							ix=(j-1)*n+i-1;
							iy=(j-1)*n+i;
					/*		d=cs[0]*V.pData[ix]+cs[1]*V.pData[iy];
							V.pData[iy]=-cs[1]*V.pData[ix]+cs[0]*V.pData[iy];
							V.pData[ix]=d;  */
							d=cs[0]*V(j-1,i-1)+cs[1]*V(j-1,i);
							V(j-1,i)=-cs[1]*V(j-1,i-1)+cs[0]*V(j-1,i);
							V(j-1,i-1)=d;
						}
					USVtDecompose_sss(fg,cs);
					s[i-1]=fg[0];
					fg[0]=cs[0]*e[i-1]+cs[1]*s[i];
					s[i]=-cs[1]*e[i-1]+cs[0]*s[i];
					fg[1]=cs[1]*e[i];
					e[i]=cs[0]*e[i];
					if (i<m)
						if ((cs[0]!=1.0)||(cs[1]!=0.0))
							for (j=1;j<=m;j++)
							{
								ix=(j-1)*m+i-1;
								iy=(j-1)*m+i;
						/*		d=cs[0]*U.pData[ix]+cs[1]*U.pData[iy];
								U.pData[iy]=-cs[1]*U.pData[ix]+cs[0]*U.pData[iy];
								U.pData[ix]=d;*/

								d=cs[0]*U(j-1,i-1)+cs[1]*U(j-1,i);
								U(j-1,i)=-cs[1]*U(j-1,i-1)+cs[0]*U(j-1,i);
								U(j-1,i-1)=d;
							}
				}
				e[mm-2]=fg[0];
				nIt--;
			}
			else
			{
				if (ks==mm)
				{
					kk++;
					fg[1]=e[mm-2]; e[mm-2]=0.0;
					for (ll=kk;ll<=mm-1;ll++)
					{
						i=mm+kk-ll-1;
						fg[0]=s[i-1];
						USVtDecompose_sss(fg,cs);
						s[i-1]=fg[0];
						if (i!=kk)
						{
							fg[1]=-cs[1]*e[i-2];
							e[i-2]=cs[0]*e[i-2];
						}
						if ((cs[0]!=1.0)||(cs[1]!=0.0))
							for (j=1;j<=n;j++)
							{
								ix=(j-1)*n+i-1;
								iy=(j-1)*n+mm-1;
						/*		d=cs[0]*V.pData[ix]+cs[1]*V.pData[iy];
								V.pData[iy]=-cs[1]*V.pData[ix]+cs[0]*V.pData[iy];
								V.pData[ix]=d;*/
								d=cs[0]*V(j-1,i-1)+cs[1]*V(j-1,mm-1);
								V(j-1,mm-1)=-cs[1]*V(j-1,i-1)+cs[0]*V(j-1,mm-1);
								V(j-1,i-1)=d;
							}
					}
				}
				else
				{
					kk=ks+1;	///////
					fg[1]=e[kk-2];
					e[kk-2]=0.0;
					for (i=kk;i<=mm;i++)
					{
						fg[0]=s[i-1];
						USVtDecompose_sss(fg,cs);
						s[i-1]=fg[0];
						fg[1]=-cs[1]*e[i-1];
						e[i-1]=cs[0]*e[i-1];
						if ((cs[0]!=1.0)||(cs[1]!=0.0))
							for (j=1;j<=m;j++)
							{
								ix=(j-1)*m+i-1;
								iy=(j-1)*m+kk-2;

							/*	d=cs[0]*U.pData[ix]+cs[1]*U.pData[iy];

								U.pData[iy]=-cs[1]*U.pData[ix]+cs[0]*U.pData[iy];

								U.pData[ix]=d; */

								d=cs[0]*U(j-1,i-1)+cs[1]*U(j-1,kk-2);
								U(j-1,kk-2)=-cs[1]*U(j-1,i-1)+cs[0]*U(j-1,kk-2);
								U(j-1,i-1)=d;
							}
					}
				}
			}
		}
	}

	return TRUE;
}



//USVtDecompose的子程序
void CMatrix::USVtDecompose_ppp(CMatrix &S, double *e, double *s, CMatrix &V, int m, int n)
{
	int i,j,p,q;
	double d;
	if (m>=n) i=n;
	else i=m;
	for (j=1;j<i;j++)
	{
	/*	S.pData[(j-1)*n+j-1]=s[j-1];
		S.pData[(j-1)*n+j]=e[j-1];  */
		S(j-1,j-1)=s[j-1];
		S(j-1,j)=e[j-1];
	}
	//S.pData[(i-1)*n+i-1]=s[i-1];
	S(i-1,i-1)=s[i-1];
	if (m<n) //S.pData[(i-1)*n+i]=e[i-1];
		S(i-1,i)=e[i-1];
	for (i=1;i<=n-1;i++)
		for (j=i+1;j<=n;j++)
		{
			p=(i-1)*n+j-1; q=(j-1)*n+i-1;
			//d=V.pData[p]; V.pData[p]=V.pData[q]; V.pData[q]=d;
			d=V(i-1,j-1);V(i-1,j-1)=V(j-1,i-1);V(j-1,i-1)=d;
		}
}

//USVtDecompose的子程序
void CMatrix::USVtDecompose_sss(double *fg, double *cs)
{
	double r,d;
	if ((fabs(fg[0])+fabs(fg[1]))==0.0)
	{
		cs[0]=1.0; cs[1]=0.0; d=0.0;
	}
	else
	{
		d=sqrt(fg[0]*fg[0]+fg[1]*fg[1]);
		if (fabs(fg[0])>fabs(fg[1]))
		{
			d=fabs(d);
			if (fg[0]<0.0) d=-d;
		}
		if (fabs(fg[1])>=fabs(fg[0]))
		{
			d=fabs(d);
			if (fg[1]<0.0) d=-d;
		}
		cs[0]=fg[0]/d;
		cs[1]=fg[1]/d;
	}
	r=1.0;
	if (fabs(fg[0])>fabs(fg[1])) r=cs[1];
	else
		if (cs[0]!=0.0) r=1.0/cs[0];

	fg[0]=d;
	fg[1]=r;
}