d8r7.cpp

Visual C++ 常用数值算法集 源代码

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