📄 d8r7.cpp
字号:
#include <math.h>
#include <iomanip.h>
#include <iostream.h>
void balanc(double a[6][6],int n)
{
double radix,sqrdx,c,r,f,g,s;
int last,i,j;
radix = 2.0;
sqrdx = radix *radix;
loop1: last = 1;
for (i = 1;i<=n;i++)
{
c = 0.0;
r = 0.0;
for (j = 1;j<=n;j++)
{
if (j != i)
{
c = c + fabs(a[j][i]);
r = r + fabs(a[i][j]);
}
}
if ((c != 0.0) && (r != 0.0))
{
g = r / radix;
f = 1.0;
s = c + r;
loop2: if (c < g)
{
f = f * radix;
c = c * sqrdx;
goto loop2;
}
g = r * radix;
loop3: if (c > g)
{
f = f / radix;
c = c / sqrdx;
goto loop3;
}
if (((c + r) / f) < (0.95 * s))
{
last = 0;
g = 1.0 / f;
for (j = 1 ;j<=n;j++)
{
a[i][j] = a[i][j] * g;
}
for (j = 1 ;j<=n;j++)
{
a[j][i] = a[j][i] * f;
}
}
}
}
if (last == 0) goto loop1;
}
void elmhes(double a[6][6],int n)
{
int m,j,i;
double x,y;
if (n > 2)
{
for (m = 2 ;m<=n - 1;m++)
{
x = 0.0;
i = m;
for (j = m; j<= n;j++)
{
if ((fabs(a[j][m - 1])) > fabs(x))
{
x = a[j][m - 1];
i = j;
}
}
if (i != m)
{
for (j = m - 1;j<=n;j++)
{
y = a[i][j];
a[i][j] = a[m][j];
a[m][j] = y;
}
for (j = 1;j<= n;j++)
{
y = a[j][i];
a[j][i] = a[j][m];
a[j][m] = y;
}
}
if (x != 0.0)
{
for (i = m + 1;i<=n;i++)
{
y = a[i][m - 1];
if (y != 0.0)
{
y = y / x;
a[i][m - 1] = y;
for (j = m;j<=n;j++)
{
a[i][j] = a[i][j] - y * a[m][j];
}
for (j = 1;j<=n;j++)
{
a[j][m] = a[j][m] + y * a[j][i];
}
}
}
}
}
}
}
void hqr(double a[6][6],int n,double wr[6],double wi[6])
{
int i,j,nn,its,l,k,iii,m,temp;
double anorm,t,s,x,y,w,p,q,z,r,v,u,bbb,aaa;
anorm = fabs(a[1][1]);
for (i = 2; i<=n; i++)
{
for (j = i - 1; j<=n; j++)
{
anorm = anorm + fabs(a[i][j]);
}
}
nn = n;
t = 0.0;
loop1: if (nn >= 1)
{
its = 0;
loop2: for (l = nn ; l >=2; l--)
{
s = fabs(a[l - 1][l - 1]) + fabs(a[l][l]);
if (s == 0.0) s = anorm;
if ((fabs(a[l][l - 1]) + s) ==s) goto loop3;
}
l = 1;
loop3: x = a[nn][nn];
if (l == nn)
{
wr[nn] = x + t;
wi[nn] = 0.0;
nn = nn - 1;
}
else
{
y = a[nn - 1][nn - 1];
w = a[nn][nn - 1] * a[nn - 1][nn];
if (l == (nn - 1))
{
p = 0.5 * (y - x);
q = p *p + w;
z = sqrt(fabs(q));
x = x + t;
if (q >= 0.0)
{
if (p>0) t=1;
if (p==0) t=0;
if (p<0) t=-1;
z = p + fabs(z) * t;
wr[nn] = z + x;
wr[nn - 1] = wr[nn];
if (z != 0.0) wr[nn] = x - w / z;
wi[nn] = 0.0;
wi[nn - 1] = 0.0;
}
else
{
wr[nn] = x + p;
wr[nn - 1] = wr[nn];
wi[nn] = z;
wi[nn - 1] = -z;
}
nn = nn - 2;
}
else
{
if (its == 30) cout<<"too many iterations"<<endl;
if ((its == 10) || (its == 20))
{
t = t + x;
for (i = 1;i<=nn;i++)
{
a[i][i] = a[i][i] - x;
}
s = fabs(a[nn][nn - 1]) + fabs(a[nn - 1][nn - 2]);
x = 0.75 * s;
y = x;
w = -0.4375 * s *s;
}
its = its + 1;
for (m = nn - 2 ; m>=l; m--)
{
z = a[m][m];
r = x - z;
s = y - z;
p = (r * s - w) / a[m + 1][m] + a[m][m + 1];
q = a[m + 1][m + 1] - z - r - s;
r = a[m + 2][m + 1];
s = fabs(p) + fabs(q) + fabs(r);
p = p / s;
q = q / s;
r = r / s;
if (m == l) goto loop4;
u = fabs(a[m][m - 1]) * (fabs(q) + fabs(r));
bbb = fabs(a[m + 1][m + 1]);
aaa = fabs(a[m - 1][m - 1]) + fabs(z) + bbb;
v = fabs(p) * aaa;
if ((u + v) == v) goto loop4;
}
loop4 : for (i = m + 2; i<=nn; i++)
{
a[i][i - 2] = 0.0;
if (i !=( m + 2)) a[i][i - 3] = 0.0;
}
for (k = m ; k<=nn - 1; k++)
{
if (k != m)
{
p = a[k][k - 1];
q = a[k + 1][k - 1];
r = 0.0;
if (k != (nn - 1)) r = a[k + 2][k - 1];
x = fabs(p) + fabs(q) + fabs(r);
if (x != 0.0)
{
p = p / x;
q = q / x;
r = r / x;
}
}
if (p>0) temp=1;
if (p==0) temp=0;
if (p<0) temp=-1;
s = sqrt(p*p + q*q + r*r) * temp;
if (s != 0.0)
{
if (k == m)
{
if (l != m ) a[k][k - 1] = -a[k][k - 1];
}
else
{
a[k][k - 1] = -s * x;
}
p = p + s;
x = p / s;
y = q / s;
z = r / s;
q = q / p;
r = r / p;
for (j = k; j<= nn; j++)
{
p = a[k][j] + q * a[k + 1][j];
if (k != (nn - 1 ))
{
p = p + r * a[k + 2][j];
a[k + 2][j] = a[k + 2][j] - p * z;
}
a[k + 1][j] = a[k + 1][j] - p * y;
a[k][j] = a[k][j] - p * x;
}
if (nn > (k + 3))
{
iii = k + 3;
}
else
{
iii = nn;
}
for (i = l; i<=iii; i++)
{
p = x * a[i][k] + y * a[i][k + 1];
if (k != (nn - 1))
{
p = p + z * a[i][k + 2];
a[i][k + 2] = a[i][k + 2] - p * r;
}
a[i][k + 1] = a[i][k + 1] - p * q;
a[i][k] = a[i][k] - p;
}
}
}
goto loop2;
}
}
goto loop1;
}
}
void main()
{
//program d8r7
//driver for routine hqr
int i,j,np = 5;
double a[6][6],wr[6],wi[6];
a[1][1]=1.0; a[1][2]=2.0; a[1][3]=0.0; a[1][4]=0.0; a[1][5]=0.0;
a[2][1]=-2.0; a[2][2]=3.0; a[2][3]=0.0; a[2][4]=0.0; a[2][5]=0.0;
a[3][1]=3.0; a[3][2]=4.0; a[3][3]=50.0; a[3][4]=0.0; a[3][5]=0.0;
a[4][1]=-4.0; a[4][2]=5.0; a[4][3]=-60.0; a[4][4]=7.0; a[4][5]=0.0;
a[5][1]=-5.0; a[5][2]=6.0; a[5][3]=-70.0; a[5][4]=8.0; a[5][5]=-9.0;
cout<< "Matrix:"<<endl;
cout<<endl;
for (i = 1;i<=np;i++)
{
cout<<setprecision(3)<<setiosflags(ios::fixed);
for (j = 1; j<= np; j++)
cout<<setw(12)<<a[i][j];
cout<<endl;
}
balanc(a,np);
elmhes(a,np);
hqr(a,np,wr,wi);
cout<<endl;
cout<<"Eigenvalues:"<<endl;
cout<<endl;
cout<<" Real Imag."<<endl;
cout<<setprecision(6)<<setiosflags(ios::fixed);
for (i = 1; i<=np; i++)
{
cout<<setw(14)<<wr[i];
cout<<setw(14)<<wi[i];
cout<<endl;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -