tqli.cpp

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

void tqli(double d[4],double  e[4],int  n, double z[4][4])
{
	int i,l,m,iter,k,t;
	double dd,g,r,s,c,p,f,b;
    if (n > 1)
	{
        for (i = 2; i<=n; i++)
		{
            e[i - 1] = e[i];
        }
        e[n] = 0.0;
        for (l = 1; l<=n; l++)
		{
            iter = 0;
loop1:      for (m = l; m<=n - 1; m++)
			{
				dd = fabs(d[m]) + fabs(d[m + 1]);
				if ((fabs(e[m]) + dd) == dd ) goto  loop2;
			}
			m = n;
loop2:		if (m != l) 
			{
				if (iter == 30 ) cout<<" too many iterations ";
				iter = iter + 1;
				g = (d[l + 1] - d[l]) / (2.0 * e[l]);
				r = sqrt(g * g + 1.0);
				if (g>0) t=1;
				if (g==0) t=0;
				if (g<0) t=-1;
				g = d[m] - d[l] + e[l] / (g + fabs(r) * t);
				s = 1.0;
				c = 1.0;
				p = 0.0;
				for (i = m - 1 ; i>=l; i--)
				{
					f = s * e[i];
					b = c * e[i];
					if (fabs(f) >= fabs(g))
					{
						c = g / f;
						r = sqrt(c *c + 1.0);
						e[i + 1] = f * r;
						s = 1.0 / r;
						c = c * s;
					}
					else
					{
						s = f / g;
						r =sqrt(s *s + 1.0);
						e[i + 1] = g * r;
						c = 1.0 / r;
						s=s*c;
					}
					g = d[i + 1] - p;
					r = (d[i] - g) * s + 2.0 * c * b;
					p = s * r;
					d[i + 1] = g + p;
					g = c * r - b;
					//omit lines from here ...
					for (k = 1; k<= n; k++)
					{
						f = z[k][i + 1];
						z[k][i + 1] = s * z[k][i] + c * f;
						z[k][i] = c * z[k][i] - s * f;
					}
					//to here when finding only eigenvalues.
				}
				d[l] = d[l] - p;
				e[l] = g;
				e[m] = 0.0;
				goto loop1;
			}
		}
    }
}