hqr.cpp

VC++常用数据算法集,里面包含了大量的常用算法程序,很有用的哟!

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