d8r1.cpp

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

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

void main()
{
    //program d8r1
    //driver for routine jacobi
	int np,ii,jj,j,k,l,kk,ll,nrot;
    double  d[4],v[4][4],r[4];
    double  a[4][4],e[4][4],ratio;
	np = 3;
    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 (ii = 1; ii<=3; ii++)
	{
        for (jj = 1; jj<=3; jj++)
		{
            e[ii][ jj] = a[ii][ jj];
        }
    }
    jacobi(e, np,d ,v, nrot);
    cout<< "Number of jacobi rotations;  "<< nrot<<endl;
    cout<< "Eigenvalues;"<<endl;
    for (j = 1; j<=np; j++)
	{
        cout<<setw(13)<<d[j];
    }
    cout<<endl;
    cout<< "Eigenvectors;"<<endl;
    for (j = 1; j<=np; j++)
	{
        cout<< setw(5)<< "Number"<<j;
        for (k = 1; k<=np; k++)
		{
            cout<< setw(11)<<v[k][j];
        }
        cout<<endl;
    }
    cout<<endl;
    //eigenvector test
    cout<< setw(5)<< "Eigenvector test";
    for (j = 1; j<=np; j++)
	{
        for (l = 1; l<=np; l++)
		{
            r[l] = 0.0;
            for (k = 1; k<=np; k++)
			{
               if (k > l)
			   {
                    kk = l;
                    ll = k;
			   }
                else
				{
                    kk = k;
                    ll = l;
                }
                r[l] = r[l] + a[ll][kk] * v[k][j];
            }
        }
        cout<<endl;
        cout<< setw(5)<<"Vector number"<<j<<endl;
        cout<<endl;
        cout<< setw(15)<<"  vector     mtrx*vec.     ratio"<<endl;
        for (l = 1; l<=np; l++)
		{
            ratio = double(r[l]) / v[l][j];
            cout<< setw(15)<<v[l][j];
            cout<< setw(15)<<r[l];
            cout<< setw(15)<<ratio<<endl;
        }
    }
}