📄 matrix.cs
字号:
{
ix=(i-1)*m+kk-1;
mtxU.elements[ix]=-mtxU.elements[ix];
}
mtxU.elements[iz]=1.0+mtxU.elements[iz];
if (kk-1>=1)
{
for (i=1; i<=kk-1; i++)
mtxU.elements[(i-1)*m+kk-1]=0.0;
}
}
else
{
for (i=1; i<=m; i++)
mtxU.elements[(i-1)*m+kk-1]=0.0;
mtxU.elements[(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.elements[ix]*mtxV.elements[iy]/mtxV.elements[iz];
}
d=-d;
for (i=kk+1; i<=n; i++)
{
ix=(i-1)*n+j-1;
iy=(i-1)*n+kk-1;
mtxV.elements[ix]=mtxV.elements[ix]+d*mtxV.elements[iy];
}
}
}
for (i=1; i<=n; i++)
mtxV.elements[(i-1)*n+kk-1]=0.0;
mtxV.elements[iz-n]=1.0;
}
for (i=1; i<=m; i++)
for (j=1; j<=n; j++)
elements[(i-1)*n+j-1]=0.0;
m1=mm;
it=60;
while (true)
{
if (mm==0)
{
ppp(elements,e,s,mtxV.elements,m,n);
return true;
}
if (it==0)
{
ppp(elements,e,s,mtxV.elements,m,n);
return false;
}
kk=mm-1;
while ((kk!=0) && (Math.Abs(e[kk-1])!=0.0))
{
d=Math.Abs(s[kk-1])+Math.Abs(s[kk]);
dd=Math.Abs(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.elements[ix]=-mtxV.elements[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.elements[ix];
mtxV.elements[ix]=mtxV.elements[iy];
mtxV.elements[iy]=d;
}
}
if (kk<m)
{
for (i=1; i<=m; i++)
{
ix=(i-1)*m+kk-1;
iy=(i-1)*m+kk;
d=mtxU.elements[ix];
mtxU.elements[ix]=mtxU.elements[iy];
mtxU.elements[iy]=d;
}
}
kk=kk+1;
}
it=60;
mm=mm-1;
}
else
{
ks=mm;
while ((ks>kk) && (Math.Abs(s[ks-1])!=0.0))
{
d=0.0;
if (ks!=mm)
d=d+Math.Abs(e[ks-1]);
if (ks!=kk+1)
d=d+Math.Abs(e[ks-2]);
dd=Math.Abs(s[ks-1]);
if (dd>eps*d)
ks=ks-1;
else
s[ks-1]=0.0;
}
if (ks==kk)
{
kk=kk+1;
d=Math.Abs(s[mm-1]);
t=Math.Abs(s[mm-2]);
if (t>d)
d=t;
t=Math.Abs(e[mm-2]);
if (t>d)
d=t;
t=Math.Abs(s[kk-1]);
if (t>d)
d=t;
t=Math.Abs(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=Math.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.elements[ix]+cs[1]*mtxV.elements[iy];
mtxV.elements[iy]=-cs[1]*mtxV.elements[ix]+cs[0]*mtxV.elements[iy];
mtxV.elements[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.elements[ix]+cs[1]*mtxU.elements[iy];
mtxU.elements[iy]=-cs[1]*mtxU.elements[ix]+cs[0]*mtxU.elements[iy];
mtxU.elements[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.elements[ix]+cs[1]*mtxV.elements[iy];
mtxV.elements[iy]=-cs[1]*mtxV.elements[ix]+cs[0]*mtxV.elements[iy];
mtxV.elements[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.elements[ix]+cs[1]*mtxU.elements[iy];
mtxU.elements[iy]=-cs[1]*mtxU.elements[ix]+cs[0]*mtxU.elements[iy];
mtxU.elements[ix]=d;
}
}
}
}
}
}
}
}
/**
* 内部函数，由SplitUV函数调用
*/
private void 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函数调用
*/
private void sss(double[] fg, double[] cs)
{
double r,d;
if ((Math.Abs(fg[0])+Math.Abs(fg[1]))==0.0)
{
cs[0]=1.0;
cs[1]=0.0;
d=0.0;
}
else
{
d=Math.Sqrt(fg[0]*fg[0]+fg[1]*fg[1]);
if (Math.Abs(fg[0])>Math.Abs(fg[1]))
{
d=Math.Abs(d);
if (fg[0]<0.0)
d=-d;
}
if (Math.Abs(fg[1])>=Math.Abs(fg[0]))
{
d=Math.Abs(d);
if (fg[1]<0.0)
d=-d;
}
cs[0]=fg[0]/d;
cs[1]=fg[1]/d;
}
r=1.0;
if (Math.Abs(fg[0])>Math.Abs(fg[1]))
r=cs[1];
else if (cs[0]!=0.0)
r=1.0/cs[0];
fg[0]=d;
fg[1]=r;
}
/**
* 求广义逆的奇异值分解法，分解成功后，原矩阵对角线元素就是矩阵的奇异值
*
* @param mtxAP - 返回原矩阵的广义逆矩阵
* @param mtxU - 返回分解后的U矩阵
* @param mtxV - 返回分解后的V矩阵
* @param eps - 计算精度
* @return bool型，求解是否成功
*/
public bool InvertUV(Matrix mtxAP, Matrix mtxU, Matrix mtxV, double eps)
{
int i,j,k,l,t,p,q,f;
// 调用奇异值分解
if (! SplitUV(mtxU, mtxV, eps))
return false;
int m = numRows;
int n = numColumns;
// 初始化广义逆矩阵
if (! mtxAP.Init(n, m))
return false;
// 计算广义逆矩阵
j=n;
if (m<n)
j=m;
j=j-1;
k=0;
while ((k<=j) && (elements[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.elements[t]=0.0;
for (l=0; l<=k; l++)
{
f=l*n+i;
p=j*m+l;
q=l*n+l;
mtxAP.elements[t]=mtxAP.elements[t]+mtxV.elements[f]*mtxU.elements[p]/elements[q];
}
}
}
return true;
}
/**
* 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法
*
* @param mtxQ - 返回豪斯荷尔德变换的乘积矩阵Q
* @param mtxT - 返回求得的对称三对角阵
* @param dblB - 一维数组，长度为矩阵的阶数，返回对称三对角阵的主对角线元素
* @param dblC - 一维数组，长度为矩阵的阶数，前n-1个元素返回对称三对角阵的
* 次对角线元素
* @return bool型，求解是否成功
*/
public bool MakeSymTri(Matrix mtxQ, Matrix mtxT, double[] dblB, double[] dblC)
{
int i,j,k,u;
double h,f,g,h2;
// 初始化矩阵Q和T
if (! mtxQ.Init(numColumns, numColumns) ||
! mtxT.Init(numColumns, numColumns))
return false;
if (dblB == null || dblC == null)
return false;
for (i=0; i<=numColumns-1; i++)
{
for (j=0; j<=numColumns-1; j++)
{
u=i*numColumns+j;
mtxQ.elements[u]=elements[u];
}
}
for (i=numColumns-1; i>=1; i--)
{
h=0.0;
if (i>1)
{
for (k=0; k<=i-1; k++)
{
u=i*numColumns+k;
h=h+mtxQ.elements[u]*mtxQ.elements[u];
}
}
if (h == 0.0)
{
dblC[i]=0.0;
if (i==1)
dblC[i]=mtxQ.elements[i*numColumns+i-1];
dblB[i]=0.0;
}
else
{
dblC[i]=Math.Sqrt(h);
u=i*numColumns+i-1;
if (mtxQ.elements[u]>0.0)
dblC[i]=-dblC[i];
h=h-mtxQ.elements[u]*dblC[i];
mtxQ.elements[u]=mtxQ.elements[u]-dblC[i];
f=0.0;
for (j=0; j<=i-1; j++)
{
mtxQ.elements[j*numColumns+i]=mtxQ.elements[i*numColumns+j]/h;
g=0.0;
for (k=0; k<=j; k++)
g=g+mtxQ.elements[j*numColumns+k]*mtxQ.elements[i*numColumns+k];
if (j+1<=i-1)
for (k=j+1; k<=i-1; k++)
g=g+mtxQ.elements[k*numColumns+j]*mtxQ.elements[i*numColumns+k];
dblC[j]=g/h;
f=f+g*mtxQ.elements[j*numColumns+i];
}
h2=f/(h+h);
for (j=0; j<=i-1; j++)
{
f=mtxQ.elements[i*numColumns+j];
g=dblC[j]-h2*f;
dblC[j]=g;
for (k=0; k<=j; k++)
{
u=j*numColumns+k;
mtxQ.elements[u]=mtxQ.elements[u]-f*dblC[k]-g*mtxQ.elements[i*numColumns+k];
}
}
dblB[i]=h;
}
}
for (i=0; i<=numColumns-2; i++)
dblC[i]=dblC[i+1];
dblC[numColumns-1]=0.0;
dblB[0]=0.0;
for (i=0; i<=numColumns-1; i++)
{
if ((dblB[i]!=(double)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.elements[i*numColumns+k]*mtxQ.elements[k*numColumns+j];
for (k=0; k<=i-1; k++)
{
u=k*numColumns+j;
mtxQ.elements[u]=mtxQ.elements[u]-g*mtxQ.elements[k*numColumns+i];
}
}
}
u=i*numColumns+i;
dblB[i]=mtxQ.elements[u]; mtxQ.elements[u]=1.0;
if (i-1>=0)
{
for (j=0; j<=i-1; j++)
{
mtxQ.elements[i*numColumns+j]=0.0;
mtxQ.elements[j*numColumns+i]=0.0;
}
}
}
// 构造对称三对角矩阵
for (i=0; i<numColumns; ++i)
{
for (j=0; j<numColumns; ++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;
}
/**
* 实对称三对角阵的全部特征值与特征向量的计算
*
* @param dblB - 一维数组，长度为矩阵的阶数，传入对称三对角阵的主对角线元素；
* 返回时存放全部特征值。
* @param dblC - 一维数组，长度为矩阵的阶数，前n-1个元素传入对称三对角阵的
* 次对角线元素
* @param mtxQ - 如果传入单位矩阵，则返回实对称三对角阵的特征值向量矩阵；
* 如果传入MakeSymTri函数求得的矩阵A的豪斯荷尔德变换的乘积
* 矩阵Q，则返回矩阵A的特征值向量矩阵。其中第i列为与数组dblB
* 中第j个特征值对应的特征向量。
* @param nMaxIt - 迭代次数
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -