d8r3.cpp

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

# include <math.h>
# include <iomanip.h>
# include <iostream.h>

void main()
{
    //program d8r3
    //driver for routine tred2
    int i,j,l,k,np = 3,m;
    double a[4][4], c[4][4], d[4], e[4], f[4][4];
    a[1][1] = 1.0; a[1][2] = 2.0; a[1][3] = 3.0;
    a[2][1] = 2.0; a[2][2] = 2.0; a[2][3] = 3.0;
    a[3][1] = 3.0; a[3][2] = 3.0; a[3][3] = 3.0;
    for (i = 1 ;i<=np;i++)
	{
        for (j = 1 ; j<=np; j++)
		{
            c[i][j] = a[i][j];
        }
    }
    tred2(c, np, d, e);
    cout<<endl;
    cout<< "Diagonal elements";
    cout<< endl;
	cout<<setprecision(5)<<setiosflags(ios::fixed);
    for (i = 1 ; i<=np; i++)
	{
        cout<<setw(15)<<d[i];
	}
    cout<< endl;
    cout<<"'Off-diagonal elements"<<endl;
    cout<< endl;
    for (i = 2 ; i<=np; i++)
	{
        cout<<setw(15)<<e[i];
    }
	cout<<endl;
    //check transformation matrix
    for (j = 1 ; j<=np; j++)
	{
        for (k = 1 ; k<=np; k++)
		{
            f[j][k] = 0.0;
            for (l = 1 ;l<=np;l++)
			{
                for (m = 1 ; m<=np; m++)
				{
                    f[j][k] = f[j][k] + c[l][j] * a[l][m] * c[m][k];
                }
            }
        }
    }
    //how does it look?
    cout<<endl;
    cout<<setw(5)<< "Tridiagonal matrix";
    cout<<endl;
   for (i = 1 ; i<=np; i++)
	{
         for (j = 1 ; j<=np; j++)
		 {
            cout<< setw(15)<<f[i][j];
		 }
        cout<<endl;
	}
}