📄 matrix.cpp
字号:
}
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]*m_pData[(i-1)*n+j-1];
for (j=kk+1; j<=n; j++)
{
for (i=kk+1; i<=m; i++)
{
ix=(i-1)*n+j-1;
m_pData[ix]=m_pData[ix]-w[i-1]*e[j-1]/e[kk];
}
}
}
for (i=kk+1; i<=n; i++)
mtxV.m_pData[(i-1)*n+kk-1]=e[i-1];
}
}
}
mm=n;
if (m+1<n)
mm=m+1;
if (k<n)
s[k]=m_pData[k*n+k];
if (m<mm)
s[mm-1]=0.0;
if (l+1<mm)
e[l]=m_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++)
mtxU.m_pData[(i-1)*m+j-1]=0.0;
mtxU.m_pData[(j-1)*m+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=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;
}
}
}
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;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -