tred2.cpp

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

void tred2(double a[4][4], int n, double d[4], double e[4])
{
	int i,j,k,l;
	double h,scale1,f,g,hh,t;
    if (n > 1)
	{
        for (i = n; i>=2; i--)
		{
            l = i - 1;
            h = 0.0;
            scale1 = 0.0;
            if (l > 1)
			{
                for (k = 1; k<=l; k++)
					scale1 = scale1 + fabs(a[i][k]);
                if (scale1 == 0.0)
					e[i] = a[i][l];
			    else
				{
                    for (k = 1; k<=l; k++)
					{
                        a[i][k] = a[i][ k] / scale1;
                        h = h + a[i][k]*a[i][k];
                    }
                    f = a[i][ l];
					if (f>0.0) t=1;
					if (f==0.0) t=0;
					if (f<0.0)  t=-1;
                    g = -sqrt(h) * t;
                    e[i] = scale1 * g;
                    h = h - f * g;
                    a[i][l] = f - g;
                    f = 0.0;
                    for (j = 1; j<=l; j++)
					{
                        a[j][i] = a[i][j] / h;
                        g = 0.0;
                        for (k = 1; k<=j; k++)
					    g = g + a[j][ k] * a[i][k];
                        if (l > j )
						{
                           for (k = j + 1; k<=l; k++)
						   g = g + a[k][j] * a[i][k];
                        }
                        e[j] = g / h;
                        f = f + e[j] * a[i][j];
                    }
                    hh = f / (h + h);
                    for (j = 1 ;j<=l;j++)
					{
                        f = a[i][j];
                        g = e[j] - hh * f;
                        e[j] = g;
                        for (k = 1;k<=j;k++)
						a[j][k] = a[j][k] - f * e[k] - g * a[i][k];
                    }
               }
			}
            else
			{
                e[i] = a[i][l];
            }
            d[i] = h;
		}
    }
    //omit following line if finding only eigenvalues.
    d[1] = 0.0;
    e[1] = 0.0;
    for (i = 1; i<=n; i++)
	{
    //delete lines from here ...
        l = i - 1;
        if (d[i] != 0.0)
		{
            for (j = 1 ; j<=l; j++)
			{
                g = 0.0;
                for (k = 1 ; k<=l; k++)
				{
                    g = g + a[i][k] * a[k][j];
                }
                for (k = 1; k<=l; k++)
				{
                    a[k][j] = a[k][j] - g * a[k][i];
                }
            }
		}
    //... to here when finding only eibenvalues.
        d[i] = a[i][ i];
    //also delete lines from here ...
        a[i][i] = 1.0;
        if (l >= 1)
		{
            for (j = 1; j<=l; j++)
			{
                a[i][j] = 0.0;
                a[j][i] = 0.0;
            }
        }
    //... to here when finding only eigenvalues.
    }
}