svector.txt

计算矩阵特征值的小程序

代码较长，如果不能执行，就是要建立结构体，大家试试吧，希望能用。

//////////////////////////////////////////////////////////////////////
// 实对称三对角阵的全部特征值与特征向量的计算
//
// 参数：
// 1. double dblB[] - 一维数组，长度为矩阵的阶数，传入对称三对角阵的主对角线元素；
// 返回时存放全部特征值。
// 2. double dblC[] - 一维数组，长度为矩阵的阶数，前n-1个元素传入对称三对角阵的次对角线元素
// 3. CMatrix& mtxQ - 如果传入单位矩阵，则返回实对称三对角阵的特征值向量矩阵；
// 如果传入MakeSymTri函数求得的矩阵A的豪斯荷尔德变换的乘积矩阵Q，则返回矩阵A的
// 特征值向量矩阵。其中第i列为与数组dblB中第j个特征值对应的特征向量。
// 4. int nMaxIt - 迭代次数，默认值为60
// 5. double eps - 计算精度，默认值为0.000001
//
// 返回值：BOOL型，求解是否成功
//////////////////////////////////////////////////////////////////////
BOOL CMatrix::SymTriEigenv(double dblB[], double dblC[], CMatrix& mtxQ, int nMaxIt /*= 60*/, double eps /*= 0.000001*/)
{
int i,j,k,m,it,u,v;
double d,f,h,g,p,r,e,s;

// 初值
int n = mtxQ.GetNumColumns();
dblC[n-1]=0.0; 
d=0.0; 
f=0.0;

// 迭代计算

for (j=0; j<=n-1; j++)
{ 
it=0;
h=eps*(fabs(dblB[j])+fabs(dblC[j]));
if (h>d) 
d=h;

m=j;
while ((m<=n-1)&&(fabs(dblC[m])>d)) 
m=m+1;

if (m!=j)
{ 
do
{ 
if (it==nMaxIt)
return FALSE;

it=it+1;
g=dblB[j];
p=(dblB[j+1]-g)/(2.0*dblC[j]);
r=sqrt(p*p+1.0);
if (p>=0.0) 
dblB[j]=dblC[j]/(p+r);
else 
dblB[j]=dblC[j]/(p-r);

h=g-dblB[j];
for (i=j+1; i<=n-1; i++)
dblB=dblB-h;

f=f+h; 
p=dblB[m]; 
e=1.0; 
s=0.0;
for (i=m-1; i>=j; i--)
{ 
g=e*dblC; 
h=e*p;
if (fabs(p)>=fabs(dblC))
{ 
e=dblC/p; 
r=sqrt(e*e+1.0);
dblC[i+1]=s*p*r; 
s=e/r; 
e=1.0/r;
}
else
{ 
e=p/dblC; 
r=sqrt(e*e+1.0);
dblC[i+1]=s*dblC*r;
s=1.0/r; 
e=e/r;
}

p=e*dblB-s*g;
dblB[i+1]=h+s*(e*g+s*dblB);
for (k=0; k<=n-1; k++)
{ 
u=k*n+i+1; 
v=u-1;
h=mtxQ.m_pData; 
mtxQ.m_pData=s*mtxQ.m_pData[v]+e*h;
mtxQ.m_pData[v]=e*mtxQ.m_pData[v]-s*h;
}
}

dblC[j]=s*p; 
dblB[j]=e*p;

} while (fabs(dblC[j])>d);
}

dblB[j]=dblB[j]+f;
}

for (i=0; i<=n-1; i++)
{ 
k=i; 
p=dblB;
if (i+1<=n-1)
{ 
j=i+1;
while ((j<=n-1)&&(dblB[j]<=p))
{ 
k=j; 
p=dblB[j]; 
j=j+1;
}
}

if (k!=i)
{ 
dblB[k]=dblB; 
dblB=p;
for (j=0; j<=n-1; j++)
{ 
u=j*n+i; 
v=j*n+k;
p=mtxQ.m_pData; 
mtxQ.m_pData=mtxQ.m_pData[v]; 
mtxQ.m_pData[v]=p;
}
}
}

return TRUE;
}

//////////////////////////////////////////////////////////////////////
// 约化一般实矩阵为赫申伯格矩阵的初等相似变换法
//
// 参数：无
//
// 返回值：无
//////////////////////////////////////////////////////////////////////
void CMatrix::MakeHberg()
{ 
int i,j,k,u,v;
double d,t;

for (k=1; k<=m_nNumColumns-2; k++)
{ 
d=0.0;
for (j=k; j<=m_nNumColumns-1; j++)
{ 
u=j*m_nNumColumns+k-1; 
t=m_pData;
if (fabs(t)>fabs(d))
{ 
d=t; 
i=j;
}
}

if (d != 0.0)
{ 
if (i!=k)
{ 
for (j=k-1; j<=m_nNumColumns-1; j++)
{ 
u=i*m_nNumColumns+j; 
v=k*m_nNumColumns+j;
t=m_pData; 
m_pData=m_pData[v]; 
m_pData[v]=t;
}

for (j=0; j<=m_nNumColumns-1; j++)
{ 
u=j*m_nNumColumns+i; 
v=j*m_nNumColumns+k;
t=m_pData; 
m_pData=m_pData[v]; 
m_pData[v]=t;
}
}

for (i=k+1; i<=m_nNumColumns-1; i++)
{ 
u=i*m_nNumColumns+k-1; 
t=m_pData/d; 
m_pData=0.0;
for (j=k; j<=m_nNumColumns-1; j++)
{ 
v=i*m_nNumColumns+j;
m_pData[v]=m_pData[v]-t*m_pData[k*m_nNumColumns+j];
}

for (j=0; j<=m_nNumColumns-1; j++)
{ 
v=j*m_nNumColumns+k;
m_pData[v]=m_pData[v]+t*m_pData[j*m_nNumColumns+i];
}
}
}
}
}


//////////////////////////////////////////////////////////////////////
// 求赫申伯格矩阵全部特征值的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;
}