操作矩阵的类 cmatrix.txt

一个操作矩阵的类CMatrix的算法

{ 
for (j=0; j<m_nNumColumns; ++j) 
{ 
mtxT.SetElement(i, j, 0); 
k = i - j; 
if (k == 0) 
mtxT.SetElement(i, j, dblB[j]); 
else if (k == 1) 
mtxT.SetElement(i, j, dblC[j]); 
else if (k == -1) 
mtxT.SetElement(i, j, dblC[i]); 
} 
} 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 实对称三对角阵的全部特征值与特征向量的计算 
// 
// 参数： 
// 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[i]=dblB[i]-h; 

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

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

mtxQ.m_pData[u]=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[i]; 
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[i]; 
dblB[i]=p; 
for (j=0; j<=n-1; j++) 
{ 
u=j*n+i; 
v=j*n+k; 
p=mtxQ.m_pData[u]; 
mtxQ.m_pData[u]=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[u]; 
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[u]; 
m_pData[u]=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[u]; 
m_pData[u]=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[u]/d; 
m_pData[u]=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; 
} 

////////////////////////////////////////////////////////////////////// 
// 求实对称矩阵特征值与特征向量的雅可比法 
// 
// 参数： 
// 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[i] = 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[u]; 
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[u]*omega; 
m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData[u]*omega; 
m_pData[u]=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[u]; 
m_pData[u]=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[u]; 
m_pData[u]=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[u]; 
mtxEigenVector.m_pData[u]=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[i] = GetElement(i,i); 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 求实对称矩阵特征值与特征向量的雅可比过关法 
// 
// 参数： 
// 1. double dblEigenValue[] - 一维数组，长度为矩阵的阶数，返回时存放特征值 
// 2. CMatrix& mtxEigenVector - 返回时存放特征向量矩阵，其中第i列为与 
// 数组dblEigenValue中第j个特征值对应的特征向量 
// 3. double eps - 计算精度，默认值为0.000001 
// 
// 返回值：BOOL型，求解是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::JacobiEigenv2(double dblEigenValue[], CMatrix& mtxEigenVector, 
double eps /*= 0.000001*/) 
{ 
int i,j,p,q,u,w,t,s; 
double ff,fm,cn,sn,omega,x,y,d; 

if (! mtxEigenVector.Init(m_nNumColumns, m_nNumColumns)) 
return FALSE; 

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; 
} 

ff=0.0; 
for (i=1; i<=m_nNumColumns-1; i++) 
{ 
for (j=0; j<=i-1; j++) 
{ 
d=m_pData[i*m_nNumColumns+j]; 
ff=ff+d*d; 
} 
} 

ff=sqrt(2.0*ff); 

Loop_0: 

ff=ff/(1.0*m_nNumColumns); 

Loop_1: 

for (i=1; i<=m_nNumColumns-1; i++) 
{ 
for (j=0; j<=i-1; j++) 
{ 
d=fabs(m_pData[i*m_nNumColumns+j]); 
if (d>ff) 
{ 
p=i; 
q=j; 
goto Loop_2; 
} 
} 
} 

if (ff<eps) 
{ 
for (i=0; i<m_nNumColumns; ++i) 
dblEigenValue[i] = GetElement(i,i); 
return TRUE; 
} 

goto Loop_0; 

Loop_2: 

u=p*m_nNumColumns+q; 
w=p*m_nNumColumns+p; 
t=q*m_nNumColumns+p; 
s=q*m_nNumColumns+q; 
x=-m_pData[u]; 
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[u]*omega; 
m_pData[s]=fm*sn*sn+m_pData[s]*cn*cn-m_pData[u]*omega; 
m_pData[u]=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[u]; 
m_pData[u]=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[u]; 
m_pData[u]=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[u]; 
mtxEigenVector.m_pData[u]=fm*cn+mtxEigenVector.m_pData[w]*sn; 
mtxEigenVector.m_pData[w]=-fm*sn+mtxEigenVector.m_pData[w]*cn; 
} 

goto Loop_1; 
}