d8r2.cpp

C语言的编程合集！是我自己编的！很有用的！大家可以下载！有意见麻烦大家提出来以便本人及时修改！

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

void eigsrt(double d[], double v[4][4], int n)
{
	double p;
	int i,j,k;
    for (i = 1; i<=n - 1; i++)
	{
        k = i;
        p = d[i];
        for (j = i + 1;j<=n;j++)
		{
            if (d[j] >= p)
			{
				k = j;
                p = d[j];
			}
        }
        if (k != i )
		{
            d[k] = d[i];
            d[i] = p;
            for (j = 1;j<=n;j++)
			{
                p = v[j][ i];
                v[j][ i] = v[j][ k];
                v[j][ k] = p;
            }
		}
    }
}

void jacobi(double a[4][4], int n, double d[], double v[4][4], int & nrot)
{
    double b[100], z[100],g,s,c,h,tau,sss,ddd,t,tresh,sm,theta;
	int ip,iq,i,j;
    for (ip = 1;ip<=n;ip++)
	{
        for (iq = 1;iq<=n;iq++)
		{
            v[ip][ iq] = 0.0;
        }
        v[ip][ ip] = 1.0;
    }
    for (ip = 1;ip<=n;ip++)
	{
        b[ip] = a[ip][ ip];
        d[ip] = b[ip];
        z[ip] = 0.0;
    }
    nrot = 0;
    for (i = 1;i<=50;i++)
	{
        sm = 0.0;
        for (ip = 1;ip<=n - 1;ip++)
		{
            for (iq = ip + 1;iq<=n;iq++)
			{
                sm = sm + fabs(a[ip][iq]);
            }
        }
        if (sm == 0.0)
		{
			return;
		}
        if (i < 4)
		{
            tresh = 0.2 * sm /double( n *n);
		}
        else
		{
            tresh = 0.0;
        }
        for (ip = 1;ip<=n - 1;ip++)
		{
            for (iq = ip + 1;iq<=n;iq++)
			{
                g = 100.0 * fabs(a[ip][ iq]);
                sss = fabs(d[ip]) + g;
                ddd = fabs(d[iq]) + g;
                if ((i > 4) && (sss == fabs(d[ip])) && (ddd == fabs(d[iq])))
				{
                    a[ip][ iq] = 0.0;
				}
                else 
				{
					if (fabs(a[ip][ iq]) > tresh)
					{
						h = d[iq] - d[ip];
						if ((fabs(h) + g) == fabs(h))
						{
                        t = a[ip][ iq] / h;
						}
					
                        else
						{
                        theta = 0.5 * h / a[ip][ iq];
                        t = 1.0 / (fabs(theta) + sqrt(1.0 + pow(theta , 2)));
						if (theta < 0.0 ) 
						t = -t;
						}
							c = 1.0 / sqrt(1.0 + t *t);
							s = t * c;
							tau = s / (1.0 + c);
							h = t * a[ip][ iq];
							z[ip] = z[ip] - h;
							z[iq] = z[iq] + h;
							d[ip] = d[ip] - h;
							d[iq] = d[iq] + h;
							a[ip][ iq] = 0.0;
							for (j = 1;j<=ip - 1;j++)
								{
								g = a[j][ ip];
								h = a[j][ iq];
								a[j][ ip] = g - s * (h + g * tau);
								a[j][ iq] = h + s * (g - h * tau);
								}
							for (j = ip + 1 ;j<=iq - 1;j++)
								{
								g = a[ip][ j];
								h = a[j][ iq];
								a[ip][ j] = g - s * (h + g * tau);
								a[j][ iq] = h + s * (g - h * tau);
								}
							for (j = iq + 1 ;j<=n;j++)
								{
								g = a[ip][ j];
								h = a[iq][ j];
								a[ip][ j] = g - s * (h + g * tau);
								a[iq][ j] = h + s * (g - h * tau);
								}
							 for (j = 1 ;j<=n;j++)
								{
								g = v[j][ ip];
								h = v[j][ iq];
								v[j][ ip] = g - s * (h + g * tau);
								v[j][ iq] = h + s * (g - h * tau);
								}
							nrot = nrot + 1;
						}
				}
			}
        }
        for (ip = 1;ip<=n;ip++)
		{
            b[ip] = b[ip] + z[ip];
            d[ip] = b[ip];
            z[ip] = 0.0;
        }
    }
cout<< " 50 iterations should never happen"<<endl;
}

void main()
{
    //program d8r2
    //driver for routine eigsrt
    int i,nrot,j,np = 3;
    double d[4], v[4][4], a[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;
    jacobi(a, np, d, v, nrot);
    cout<< "Unsorted eigenvectors:"<<endl;
    for (i = 1; i<=np; i++)
	{
        cout<< setw(5)<<endl;
        cout<< setw(5)<<"Eigenvalue "<< i<<" = ";
        cout<< d[i];
        cout<<endl;
        cout<< setw(5)<< "Eigenvector:  "<<endl;
        for (j = 1; j<=np; j++)
		{
            cout<<setw(15)<<v[j][ i];
        }
		cout<<endl;
    }
    cout<< endl;
    cout<< setw(5)<< "********** sorting **********"<<endl;
    eigsrt(d, v, np);
	cout<<endl;
    cout<< setw(5)<< "Sorted eigenvectors:"<<endl;
    for (i = 1; i<=np; i++)
	{
        cout<< endl;
        cout<< setw(5)<<"Eigenvalue"<< i<< " = ";
        cout<< d[i]<<endl;
        cout<< setw(5)<< "Eigenvector:"<<endl;
        for (j = 1; j<=np; j++)
		{
            cout<< setw(15)<<v[j][ i];
        }
         cout<< endl;
    }
}