操作矩阵的类 cmatrix.txt

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

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=d+mtxU.m_pData[ix]*mtxU.m_pData[iy]/mtxU.m_pData[iz]; 
} 

d=-d; 
for (i=kk; i<=m; i++) 
{ 
ix=(i-1)*m+j-1; 
iy=(i-1)*m+kk-1; 

mtxU.m_pData[ix]=mtxU.m_pData[ix]+d*mtxU.m_pData[iy]; 
} 
} 
} 

for (i=kk; i<=m; i++) 
{ 
ix=(i-1)*m+kk-1; 
mtxU.m_pData[ix]=-mtxU.m_pData[ix]; 
} 

mtxU.m_pData[iz]=1.0+mtxU.m_pData[iz]; 
if (kk-1>=1) 
{ 
for (i=1; i<=kk-1; i++) 

mtxU.m_pData[(i-1)*m+kk-1]=0.0; 
} 
} 
else 
{ 
for (i=1; i<=m; i++) 
mtxU.m_pData[(i-1)*m+kk-1]=0.0; 
mtxU.m_pData[(kk-1)*m+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+mtxV.m_pData[ix]*mtxV.m_pData[iy]/mtxV.m_pData[iz]; 
} 

d=-d; 
for (i=kk+1; i<=n; i++) 
{ 
ix=(i-1)*n+j-1; 
iy=(i-1)*n+kk-1; 
mtxV.m_pData[ix]=mtxV.m_pData[ix]+d*mtxV.m_pData[iy]; 
} 
} 
} 

for (i=1; i<=n; i++) 
mtxV.m_pData[(i-1)*n+kk-1]=0.0; 

mtxV.m_pData[iz-n]=1.0; 
} 

for (i=1; i<=m; i++) 
for (j=1; j<=n; j++) 
m_pData[(i-1)*n+j-1]=0.0; 

m1=mm; 
it=60; 
while (TRUE) 
{ 
if (mm==0) 
{ 
ppp(m_pData,e,s,mtxV.m_pData,m,n); 
return TRUE; 
} 
if (it==0) 
{ 
ppp(m_pData,e,s,mtxV.m_pData,m,n); 
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>eps*d) 
kk=kk-1; 
else 
e[kk-1]=0.0; 
} 

if (kk==mm-1) 
{ 
kk=kk+1; 
if (s[kk-1]<0.0) 
{ 
s[kk-1]=-s[kk-1]; 
for (i=1; i<=n; i++) 
{ 
ix=(i-1)*n+kk-1; 
mtxV.m_pData[ix]=-mtxV.m_pData[ix];} 
} 

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=mtxV.m_pData[ix]; 

mtxV.m_pData[ix]=mtxV.m_pData[iy]; 
mtxV.m_pData[iy]=d; 
} 
} 

if (kk<m) 
{ 
for (i=1; i<=m; i++) 
{ 
ix=(i-1)*m+kk-1; 
iy=(i-1)*m+kk; 
d=mtxU.m_pData[ix]; 

mtxU.m_pData[ix]=mtxU.m_pData[iy]; 
mtxU.m_pData[iy]=d; 
} 
} 

kk=kk+1; 
} 

it=60; 
mm=mm-1; 
} 
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>eps*d) 
ks=ks-1; 
else 
s[ks-1]=0.0; 
} 

if (ks==kk) 
{ 
kk=kk+1; 
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++) 
{ 
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]*mtxV.m_pData[ix]+cs[1]*mtxV.m_pData[iy]; 

mtxV.m_pData[iy]=-cs[1]*mtxV.m_pData[ix]+cs[0]*mtxV.m_pData[iy]; 
mtxV.m_pData[ix]=d; 
} 
} 

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]*mtxU.m_pData[ix]+cs[1]*mtxU.m_pData[iy]; 

mtxU.m_pData[iy]=-cs[1]*mtxU.m_pData[ix]+cs[0]*mtxU.m_pData[iy]; 

mtxU.m_pData[ix]=d; 
} 
} 
} 
} 

e[mm-2]=fg[0]; 
it=it-1; 
} 
else 
{ 
if (ks==mm) 
{ 
kk=kk+1; 
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]; 
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]*mtxV.m_pData[ix]+cs[1]*mtxV.m_pData[iy]; 

mtxV.m_pData[iy]=-cs[1]*mtxV.m_pData[ix]+cs[0]*mtxV.m_pData[iy]; 

mtxV.m_pData[ix]=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]; 
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]*mtxU.m_pData[ix]+cs[1]*mtxU.m_pData[iy]; 

mtxU.m_pData[iy]=-cs[1]*mtxU.m_pData[ix]+cs[0]*mtxU.m_pData[iy]; 

mtxU.m_pData[ix]=d; 
} 
} 
} 
} 
} 
} 
} 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 内部函数，由SplitUV函数调用 
////////////////////////////////////////////////////////////////////// 
void CMatrix::ppp(double a[], double e[], double s[], double 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-1; j++) 
{ 
a[(j-1)*n+j-1]=s[j-1]; 
a[(j-1)*n+j]=e[j-1]; 
} 

a[(i-1)*n+i-1]=s[i-1]; 
if (m<n) 
a[(i-1)*n+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[p]; 
v[p]=v[q]; 
v[q]=d; 
} 
} 
} 

////////////////////////////////////////////////////////////////////// 
// 内部函数，由SplitUV函数调用 
////////////////////////////////////////////////////////////////////// 
void CMatrix::sss(double fg[2], double cs[2]) 
{ 
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; 
} 

////////////////////////////////////////////////////////////////////// 
// 求广义逆的奇异值分解法，分解成功后，原矩阵对角线元素就是矩阵的奇异值 
// 
// 参数： 
// 1. CMatrix& mtxAP - 返回原矩阵的广义逆矩阵 
// 2. CMatrix& mtxU - 返回分解后的U矩阵 
// 3. CMatrix& mtxV - 返回分解后的V矩阵 
// 4. double eps - 计算精度，默认值为0.000001 
// 
// 返回值：BOOL型，求解是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::GInvertUV(CMatrix& mtxAP, CMatrix& mtxU, CMatrix& mtxV, double 
eps /*= 0.000001*/) 
{ 
int i,j,k,l,t,p,q,f; 

// 调用奇异值分解 
if (! SplitUV(mtxU, mtxV, eps)) 
return FALSE; 

int m = m_nNumRows; 
int n = m_nNumColumns; 

// 初始化广义逆矩阵 
if (! mtxAP.Init(n, m)) 
return FALSE; 

// 计算广义逆矩阵 

j=n; 
if (m<n) 
j=m; 
j=j-1; 
k=0; 
while ((k<=j)&&(m_pData[k*n+k]!=0.0)) 
k=k+1; 

k=k-1; 
for (i=0; i<=n-1; i++) 
{ 
for (j=0; j<=m-1; j++) 
{ 
t=i*m+j; 
mtxAP.m_pData[t]=0.0; 
for (l=0; l<=k; l++) 
{ 
f=l*n+i; 
p=j*m+l; 
q=l*n+l; 

mtxAP.m_pData[t]=mtxAP.m_pData[t]+mtxV.m_pData[f]*mtxU.m_pData[p]/m_pData[q]; 
} 
} 
} 

return TRUE; 
} 

////////////////////////////////////////////////////////////////////// 
// 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法 
// 
// 参数： 
// 1. CMatrix& mtxQ - 返回豪斯荷尔德变换的乘积矩阵Q 
// 2. CMatrix& mtxT - 返回求得的对称三对角阵 
// 3. double dblB[] - 一维数组，长度为矩阵的阶数，返回对称三对角阵的主对角线元 
素 
// 4. double dblC[] - 一维数组，长度为矩阵的阶数，前n-1个元素返回对称三对角阵的 
次对角线元素 
// 
// 返回值：BOOL型，求解是否成功 
////////////////////////////////////////////////////////////////////// 
BOOL CMatrix::MakeSymTri(CMatrix& mtxQ, CMatrix& mtxT, double dblB[], double 
dblC[]) 
{ 
int i,j,k,u; 
double h,f,g,h2; 

// 初始化矩阵Q和T 
if (! mtxQ.Init(m_nNumColumns, m_nNumColumns) || 
! mtxT.Init(m_nNumColumns, m_nNumColumns)) 
return FALSE; 

if (dblB == NULL || dblC == NULL) 
return FALSE; 

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

for (i=m_nNumColumns-1; i>=1; i--) 
{ 
h=0.0; 
if (i>1) 
{ 
for (k=0; k<=i-1; k++) 
{ 
u=i*m_nNumColumns+k; 
h=h+mtxQ.m_pData[u]*mtxQ.m_pData[u]; 
} 
} 

if (h == 0.0) 
{ 
dblC[i]=0.0; 
if (i==1) 
dblC[i]=mtxQ.m_pData[i*m_nNumColumns+i-1]; 
dblB[i]=0.0; 
} 
else 
{ 
dblC[i]=sqrt(h); 
u=i*m_nNumColumns+i-1; 
if (mtxQ.m_pData[u]>0.0) 
dblC[i]=-dblC[i]; 

h=h-mtxQ.m_pData[u]*dblC[i]; 
mtxQ.m_pData[u]=mtxQ.m_pData[u]-dblC[i]; 
f=0.0; 
for (j=0; j<=i-1; j++) 
{ 

mtxQ.m_pData[j*m_nNumColumns+i]=mtxQ.m_pData[i*m_nNumColumns+j]/h; 
g=0.0; 
for (k=0; k<=j; k++) 

g=g+mtxQ.m_pData[j*m_nNumColumns+k]*mtxQ.m_pData[i*m_nNumColumns+k]; 

if (j+1<=i-1) 
for (k=j+1; k<=i-1; k++) 

g=g+mtxQ.m_pData[k*m_nNumColumns+j]*mtxQ.m_pData[i*m_nNumColumns+k]; 

dblC[j]=g/h; 
f=f+g*mtxQ.m_pData[j*m_nNumColumns+i]; 
} 

h2=f/(h+h); 
for (j=0; j<=i-1; j++) 
{ 
f=mtxQ.m_pData[i*m_nNumColumns+j]; 
g=dblC[j]-h2*f; 
dblC[j]=g; 
for (k=0; k<=j; k++) 
{ 
u=j*m_nNumColumns+k; 

mtxQ.m_pData[u]=mtxQ.m_pData[u]-f*dblC[k]-g*mtxQ.m_pData[i*m_nNumColumns+k]; 
} 
} 

dblB[i]=h; 
} 
} 

for (i=0; i<=m_nNumColumns-2; i++) 
dblC[i]=dblC[i+1]; 

dblC[m_nNumColumns-1]=0.0; 
dblB[0]=0.0; 
for (i=0; i<=m_nNumColumns-1; i++) 
{ 
if ((dblB[i]!=0.0)&&(i-1>=0)) 
{ 
for (j=0; j<=i-1; j++) 
{ 
g=0.0; 
for (k=0; k<=i-1; k++) 

g=g+mtxQ.m_pData[i*m_nNumColumns+k]*mtxQ.m_pData[k*m_nNumColumns+j]; 

for (k=0; k<=i-1; k++) 
{ 
u=k*m_nNumColumns+j; 

mtxQ.m_pData[u]=mtxQ.m_pData[u]-g*mtxQ.m_pData[k*m_nNumColumns+i]; 
} 
} 
} 

u=i*m_nNumColumns+i; 
dblB[i]=mtxQ.m_pData[u]; mtxQ.m_pData[u]=1.0; 
if (i-1>=0) 
{ 
for (j=0; j<=i-1; j++) 
{ 
mtxQ.m_pData[i*m_nNumColumns+j]=0.0; 
mtxQ.m_pData[j*m_nNumColumns+i]=0.0; 
} 
} 
} 

// 构造对称三对角矩阵 
for (i=0; i<m_nNumColumns; ++i)