d8r6.cpp

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 main()
{
    //program d8r6
    //driver for routine elmhes
    int i,j,np = 5;
    double  a[6][6],r[6],c[6];
    a[1][1]=1.0; a[1][2]=2.0; a[1][3]=300.0; a[1][4]=4.0; a[1][5]=5.0;
    a[2][1]=2.0; a[2][2]=3.0; a[2][3]=400.0; a[2][4]=5.0; a[2][5]=6.0;
    a[3][1]=3.0; a[3][2]=4.0; a[3][3]=5.0;   a[3][4]=6.0; a[3][5]=7.0;
    a[4][1]=4.0; a[4][2]=5.0; a[4][3]=600.0; a[4][4]=7.0; a[4][5]=8.0;
    a[5][1]=5.0; a[5][2]=6.0; a[5][3]=700.0; a[5][4]=8.0; a[5][5]=9.0;
    cout<< endl;
    cout<<"***** Original matrix *****"<< endl;
    cout<< endl;
    for (i = 1; i<=np; i++)
	{
        for (j = 1; j<=np; j++)
        cout<< setw(9)<<a[i][j];
	    cout<<endl;
	}
    cout<< endl;
    cout<<"***** balance matrix *****";
    cout<< endl;
    balanc(a,np);
    for (i = 1; i<=np; i++)
	{
        for (j = 1; j<=np; j++)
        cout<< setw(9)<<a[i][j];
	    cout<<endl;
	}
    cout<< endl;
    cout<< "*****reduce to Hessenherg form *****";
    cout<< endl;
    elmhes(a,np);
    for (j = 1; j<=np - 2; j++)
	{
        for (i = j + 2; i<= np; i++)
		{
        a[i][j] = 0.0;
		}
	}
	cout<<setprecision(4);
    for (i = 1; i<=np; i++)
	{
        for (j = 1; j<=np; j++)
		{
        cout<< setw(13)<<a[i][j];
		}
        cout<<endl;
	}
}