- //////////////////////////////////////////////////////////////////////
- // 求赫申伯格矩阵全部特征值的QR方法
- //
- // 参数:
- // 1. double dblU[] - 一维数组,长度为矩阵的阶数,返回时存放特征值的实部
- // 2. double dblV[] - 一维数组,长度为矩阵的阶数,返回时存放特征值的虚部
- // 3. int nMaxIt - 迭代次数,默认值为60
- // 4. double eps - 计算精度,默认值为0.000001
- //
- // 返回值:BOOL型,求解是否成功
- //////////////////////////////////////////////////////////////////////
- BOOL CMatrix::HBergEigenv(double dblU[], double dblV[], int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
- {
- int m,it,i,j,k,l,ii,jj,kk,ll;
- double b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
- int n = m_nNumColumns;
- it=0;
- m=n;
- while (m!=0)
- {
- l=m-1;
- while ((l>0)&&(fabs(m_pData[l*n+l-1]) >
- eps*(fabs(m_pData[(l-1)*n+l-1])+fabs(m_pData[l*n+l]))))
- l=l-1;
- ii=(m-1)*n+m-1;
- jj=(m-1)*n+m-2;
- kk=(m-2)*n+m-1;
- ll=(m-2)*n+m-2;
- if (l==m-1)
- {
- dblU[m-1]=m_pData[(m-1)*n+m-1];
- dblV[m-1]=0.0;
- m=m-1;
- it=0;
- }
- else if (l==m-2)
- {
- b=-(m_pData[ii]+m_pData[ll]);
- c=m_pData[ii]*m_pData[ll]-m_pData[jj]*m_pData[kk];
- w=b*b-4.0*c;
- y=sqrt(fabs(w));
- if (w>0.0)
- {
- xy=1.0;
- if (b<0.0)
- xy=-1.0;
- dblU[m-1]=(-b-xy*y)/2.0;
- dblU[m-2]=c/dblU[m-1];
- dblV[m-1]=0.0; dblV[m-2]=0.0;
- }
- else
- {
- dblU[m-1]=-b/2.0;
- dblU[m-2]=dblU[m-1];
- dblV[m-1]=y/2.0;
- dblV[m-2]=-dblV[m-1];
- }
- m=m-2;
- it=0;
- }
- else
- {
- if (it>=nMaxIt)
- return FALSE;
- it=it+1;
- for (j=l+2; j<=m-1; j++)
- m_pData[j*n+j-2]=0.0;
- for (j=l+3; j<=m-1; j++)
- m_pData[j*n+j-3]=0.0;
- for (k=l; k<=m-2; k++)
- {
- if (k!=l)
- {
- p=m_pData[k*n+k-1];
- q=m_pData[(k+1)*n+k-1];
- r=0.0;
- if (k!=m-2)
- r=m_pData[(k+2)*n+k-1];
- }
- else
- {
- x=m_pData[ii]+m_pData[ll];
- y=m_pData[ll]*m_pData[ii]-m_pData[kk]*m_pData[jj];
- ii=l*n+l;
- jj=l*n+l+1;
- kk=(l+1)*n+l;
- ll=(l+1)*n+l+1;
- p=m_pData[ii]*(m_pData[ii]-x)+m_pData[jj]*m_pData[kk]+y;
- q=m_pData[kk]*(m_pData[ii]+m_pData[ll]-x);
- r=m_pData[kk]*m_pData[(l+2)*n+l+1];
- }
- if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
- {
- xy=1.0;
- if (p<0.0)
- xy=-1.0;
- s=xy*sqrt(p*p+q*q+r*r);
- if (k!=l)
- m_pData[k*n+k-1]=-s;
- e=-q/s;
- f=-r/s;
- x=-p/s;
- y=-x-f*r/(p+s);
- g=e*r/(p+s);
- z=-x-e*q/(p+s);
- for (j=k; j<=m-1; j++)
- {
- ii=k*n+j;
- jj=(k+1)*n+j;
- p=x*m_pData[ii]+e*m_pData[jj];
- q=e*m_pData[ii]+y*m_pData[jj];
- r=f*m_pData[ii]+g*m_pData[jj];
- if (k!=m-2)
- {
- kk=(k+2)*n+j;
- p=p+f*m_pData[kk];
- q=q+g*m_pData[kk];
- r=r+z*m_pData[kk];
- m_pData[kk]=r;
- }
- m_pData[jj]=q; m_pData[ii]=p;
- }
- j=k+3;
- if (j>=m-1)
- j=m-1;
- for (i=l; i<=j; i++)
- {
- ii=i*n+k;
- jj=i*n+k+1;
- p=x*m_pData[ii]+e*m_pData[jj];
- q=e*m_pData[ii]+y*m_pData[jj];
- r=f*m_pData[ii]+g*m_pData[jj];
- if (k!=m-2)
- {
- kk=i*n+k+2;
- p=p+f*m_pData[kk];
- q=q+g*m_pData[kk];
- r=r+z*m_pData[kk];
- m_pData[kk]=r;
- }
- m_pData[jj]=q;
- m_pData[ii]=p;
- }
- }
- }
- }
- }
- return TRUE;
- }
- //////////////////////////////////////////////////////////////////////
- // 求实对称矩阵特征值与特征向量的雅可比法
- //
- // 参数:
- // 1. double dblEigenValue[] - 一维数组,长度为矩阵的阶数,返回时存放特征值
- // 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵,其中第i列为与
- // 数组dblEigenValue中第j个特征值对应的特征向量
- // 3. int nMaxIt - 迭代次数,默认值为60
- // 4. double eps - 计算精度,默认值为0.000001
- //
- // 返回值:BOOL型,求解是否成功
- //////////////////////////////////////////////////////////////////////
- BOOL CMatrix::JacobiEigenv(double dblEigenValue[], CMatrix& mtxEigenVector, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
- {
- int i,j,p,q,u,w,t,s,l;
- double fm,cn,sn,omega,x,y,d;
- if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns))
- return FALSE;
- l=1;
- for (i=0; i<=m_nNumColumns-1; i++)
- {
- mtxEigenVector.m_pData[i*m_nNumColumns+i]=1.0;
- for (j=0; j<=m_nNumColumns-1; j++)
- if (i!=j)
- mtxEigenVector.m_pData[i*m_nNumColumns+j]=0.0;
- }
- while (TRUE)
- {
- fm=0.0;
- for (i=1; i<=m_nNumColumns-1; i++)
- {
- for (j=0; j<=i-1; j++)
- {
- d=fabs(m_pData[i*m_nNumColumns+j]);
- if ((i!=j)&&(d>fm))
- {
- fm=d;
- p=i;
- q=j;
- }
- }
- }
- if (fm<eps)
- {
- for (i=0; i<m_nNumColumns; ++i)
- dblEigenValue = GetElement(i,i);
- return TRUE;
- }
- if (l>nMaxIt)
- return FALSE;
- l=l+1;
- u=p*m_nNumColumns+q;
- w=p*m_nNumColumns+p;
- t=q*m_nNumColumns+p;
- s=q*m_nNumColumns+q;
- x=-m_pData;
- y=(m_pData[s]-m_pData[w])/2.0;
- omega=x/sqrt(x*x+y*y);
- if (y<0.0)
- omega=-omega;
- sn=1.0+sqrt(1.0-omega*omega);
- sn=omega/sqrt(2.0*sn);
- cn=sqrt(1.0-sn*sn);
- fm=m_pData[w];
- m_pData[w]=fm*cn*cn+m_pData[s]*sn*sn+m_pData*omega;
- m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData*omega;
- m_pData=0.0;
- m_pData[t]=0.0;
- for (j=0; j<=m_nNumColumns-1; j++)
- {
- if ((j!=p)&&(j!=q))
- {
- u=p*m_nNumColumns+j; w=q*m_nNumColumns+j;
- fm=m_pData;
- m_pData=fm*cn+m_pData[w]*sn;
- m_pData[w]=-fm*sn+m_pData[w]*cn;
- }
- }
- for (i=0; i<=m_nNumColumns-1; i++)
- {
- if ((i!=p)&&(i!=q))
- {
- u=i*m_nNumColumns+p;
- w=i*m_nNumColumns+q;
- fm=m_pData;
- m_pData=fm*cn+m_pData[w]*sn;
- m_pData[w]=-fm*sn+m_pData[w]*cn;
- }
- }
- for (i=0; i<=m_nNumColumns-1; i++)
- {
- u=i*m_nNumColumns+p;
- w=i*m_nNumColumns+q;
- fm=mtxEigenVector.m_pData;
- mtxEigenVector.m_pData=fm*cn+mtxEigenVector.m_pData[w]*sn;
- mtxEigenVector.m_pData[w]=-fm*sn+mtxEigenVector.m_pData[w]*cn;
- }
- }
- for (i=0; i<m_nNumColumns; ++i)
- dblEigenValue = GetElement(i,i);
- return TRUE;
- }
复制代码
[ 本帖最后由 风花雪月 于 2006-11-14 19:23 编辑 ] |