1005ac.txt

数值分析的c语言几种算法实现

#include <cmath>
#include <iostream>
using namespace std;

int EigenV(int n, double a[][MAX_SIZE], double *lambda, double v[], 

double TOL, int MAXN)
{
	double b[MAX_SIZE][MAX_SIZE];
	int i,j,r;
	double q;
	int p;
	int k;
	double max=0.0,mm=0.0;
	double x[MAX_SIZE];
	double temp1=0.0,temp2=0.0,value=0.0,tmp;
	double u[MAX_SIZE];
	
	
	
	int m=0;
	for(i=0;i<n;i++) x[i]=0.0;
	for(i=0;i<n;i++) 
	{
		for(j=0;j<n;j++)
			x[i]+=(v[j]*a[j][i]);
	}
	tmp=0.0;
	for(i=0;i<n;i++) tmp+=x[i]*v[i];
	for(i=0;i<n;i++) m+=v[i]*v[i];
	q=*lambda;
	//q=tmp/m;
    //cout << tmp << endl << m << endl;
	//cout << q << endl;
	max=0.0;
	p=0;
   	for(k=0;k<n;k++)
		if(fabs(v[k])>fabs(max))
		{
				max=v[k];
				p=k;
		}
    //cout << max << endl << p << endl;
	for(i=0;i<n;i++){ v[i]=v[i]/max;x[i]=v[i];}
    //   cout << x[p] << endl;
	for(i=0;i<MAXN;i++)
	{
		for(r=0;r<n;r++)
			for(j=0;j<n;j++)
				b[r][j]=a[r][j];

		for(k=0;k<n;k++)
			b[k][k]-=q;

		for(r=0;r<n;r++){ u[r]=x[r];}

		for(j=0;j<n-1;j++)
		{
			for(k=j+1;k<n;k++)
			{
				tmp=b[k][j]/b[j][j];
				for(r=j+1;r<n;r++)
					b[k][r]-=(tmp*b[j][r]);
				u[k]-=(tmp*u[j]);
			}
		}
		//if((b[n-1][n-1]==0.0&&u[n-1]==0.0)||b[n-1][n-1]==0.0) return -1;

		v[n-1]=u[n-1]/b[n-1][n-1];
		for(j=n-2;j>=0;j--)
		{
			tmp=0;		
			for(k=j+1;k<n;k++)
			{
				tmp+=v[k]*b[j][k];
			}
			v[j]=(u[j]-tmp)/b[j][j];
		}
		//for(j=0;j<n;j++) cout << v[j] << " ";cout << endl;
		 
		 max=v[p];
		 mm=0.0;
		 for(k=0;k<n;k++)
			if(fabs(v[k])>fabs(mm))
			{
				mm=v[k];
				p=k;
			}
		tmp=0;
		for(k=0;k<n;k++)
		if(fabs(x[k]-v[k]/mm)>tmp)
			tmp=fabs(x[k]-v[k]/mm);
		for(k=0;k<n;k++)
			x[k]=v[k]/mm;
		if(tmp<=TOL) {
			if(1/max+q==*lambda) return -1;
			*lambda=1/max+q;
			for(k=0;k<n;k++)
				v[k]=x[k];
			return 1;
		}
		if(i==MAXN-1) 
			return 0;
		
	}
}