sbqr.cpp

这个是数值分析中双步QR算法的实现程序

#include "stdio.h"
#include "math.h"
#define n 3
#define eps 1e-12
#define L 100
void print(char arr[2],double A[n][n])
{
	printf("%s is:\n",arr);
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
		{	
			
			printf("%12.7f,",A[i][j]);
			if(j==n-1) printf("\n");
		}
}
//打印函数
void print1(char ar[2],double *arr)
{
	printf("%s is:",ar);
	for(int j=0;j<n;j++)				
	{	
		if(j==0) printf("\t[");
		printf("%14.9f",*(arr+j));
		if(j!=n-1) printf(",");
		else printf("\t]");
	}
	printf("\n");
}
void chenfa(double arr1[n][n],double arr2[n][n],double re[n][n])
{
	double temp;
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
		{
			temp=0;
			for(int k=0;k<n;k++)
				temp=temp+arr1[i][k]*arr2[k][j];
			re[i][j]=temp;
		}
}
int sign(double x)
{
	if(x>0) return 1;
	else return -1;
}
double neiji(double arr1[n],double arr2[n])
{
	double value=0;
	for(int i=0;i<n;i++)
		value=value+arr1[i]*arr2[i];
	return value;
}
void Sanjiaohua(double A[n][n])
{
	int r,i,j;
	double dr,cr,hr,tr;
	double I[n][n],H[n][n];
	double u[n],w[n],p[n],q[n];
	int pd;
	double temp;
	for(i=0;i<n;i++)//Q is I
		for(j=0;j<n;j++)
		{
			if(i==j) I[i][j]=1;
			else I[i][j]=0;
		}
	for(r=1;r<=n-2;r++)
	{
		for(i=r+2;i<=n;i++)
		{
			if(A[i-1][r-1]==0) 
			{
				if(i==n) pd=1;
				else continue;
			}
			else {pd=0;break;}
		}
		if(pd==1) continue;
		else
		{
			for(temp=0,i=r+1;i<=n;i++)
				temp=temp+A[i-1][r-1]*A[i-1][r-1];
			dr=sqrt(temp);
			cr=-sign(A[r][r-1])*dr;
			hr=cr*cr-cr*A[r][r-1];
			for(i=1;i<=n;i++)//ur
			{
				if(i<r+1) u[i-1]=0;
				else if(i==r+1) u[i-1]=A[r][r-1]-cr;
				else u[i-1]=A[i-1][r-1];
			}
			for(i=0;i<n;i++)//pr
			{
				temp=0;
				for(j=0;j<n;j++)	
					temp=temp+A[j][i]*u[j];
				p[i]=temp/hr;
			}
			for(i=0;i<n;i++)//qr
			{
				temp=0;
				for(j=0;j<n;j++)	
					temp=temp+A[i][j]*u[j];
				q[i]=temp/hr;
			}
			tr=neiji(p,u)/hr;//tr
			for(i=0;i<n;i++)//wr
			{
				w[i]=q[i]-tr*u[i];			
			}
			for(i=0;i<n;i++)//A(r+1)
			{
				for(j=0;j<n;j++)	
					A[i][j]=A[i][j]-w[i]*u[j]-u[i]*p[j];
			}
			/*for(i=0;i<n;i++)
				for(j=0;j<n;j++)
				{
					H[i][j]=I[i][j]-2*u[i]*u[j]/(neiji(u,u));
				}
			print("H",H);*/
		}//else
	}//for(r)
	print("A",A);
}
void QR(double A[n][n],int m)
{
	int r,i,j;
	double Q[n][n];
	double dr,cr,hr;
	double u[n],w[n],p[n],H[n][n];
	int pd;
	double temp;
	for(i=0;i<n;i++)//Q is I
		for(j=0;j<n;j++)
		{
			if(i==j) Q[i][j]=1;
			else Q[i][j]=0;
		}
	for(r=1;r<=m-1;r++)
	{
		for(i=r+1;i<=m;i++)
		{
			if(A[i-1][r-1]==0) 
			{
				if(i==n) pd=1;
				else continue;
			}
			else {pd=0;break;}
		}
		if(pd==1) continue;
		else
		{
			for(temp=0,i=r;i<=m;i++)
				temp=temp+A[i-1][r-1]*A[i-1][r-1];
			dr=sqrt(temp);
			cr=-sign(A[r-1][r-1])*dr;
			hr=cr*cr-cr*A[r-1][r-1];
			for(i=1;i<=m;i++)//ur
			{
				if(i<r) u[i-1]=0;
				else if(i==r) u[i-1]=A[r-1][r-1]-cr;
				else u[i-1]=A[i-1][r-1];
			}
			/*for(i=0;i<n;i++)//wr
			{
				temp=0;
				for(j=0;j<n;j++)	
					temp=temp+Q[i][j]*u[j];
				w[i]=temp;
			}
			for(i=0;i<n;i++)//Q(r+1)
			{
				for(j=0;j<n;j++)	
					Q[i][j]=Q[i][j]-(w[i]*u[j])/hr;
			}
			for(i=0;i<n;i++)//pr
			{
				temp=0;
				for(j=0;j<n;j++)	
					temp=temp+A[j][i]*u[j];
				p[i]=temp/hr;
			}
			for(i=0;i<n;i++)//A(r+1)
			{
				for(j=0;j<n;j++)	
					A[i][j]=A[i][j]-u[i]*p[j];
			}*/
			for(i=0;i<m;i++)//v
			{
				temp=0;
				for(j=0;j<m;j++)	
					temp=temp+A[j][i]*u[j];
				v[i]=temp/hr;
			}
			for(i=0;i<m;i++)//M
				for(j=0;j<m;j++)	
					A[i][j]=A[i][j]-(u[i]*v[j]);
			for(i=0;i<n;i++)//pr
			{
				temp=0;
				for(j=0;j<n;j++)	
					temp=temp+A[j][i]*u[j];
				p[i]=temp/hr;
			}

			//printf("dr is %10.5f,cr is %10.5f,hr is %10.5f\n",dr,cr,hr);
			//print1("u",u);
			//print1("w",w);
			//print("Q",Q);		
			//print1("p",p);
			//print("A",A);
			print("H",H);

		}//else
	}//for(r)
	print("Q",Q);
	print("R",A);
	//chenfa(Q,A);

}
void fanchen(double a11,double a12,double a21,double a22,double result[3])
{
	double delta,det;
	det=a11*a22-a12*a21;
	delta=(a11+a22)*(a11+a22)-4*det;
	if(delta>=0) 
	{
		result[0]=(a11+a22+sqrt(detlta))/2;
		result[1]=(a11+a22-sqrt(detlta))/2;
		result[2]=0;//实数根
		printf("x1=%10.5f,x2=%10.5f\n",result[0],result[1]);
	}
	else
	{
		result[0]=(a11+a22)/2;
		result[1]=sqrt(detlta)/2;
		result[2]=1;//复数根
		printf("x1=%10.5f+j%10.5f,x2=%10.5f+j%10.5f\n",result[0],result[1],result[0],result[1]);
	}
}
int main()
{
	int i,j,m,k;
	double A[n][n]={3,1,0,1,2,1,0,1,1};
	double nmd[n],jieguo[3],s,t,I[n][n];
	for(i=0;i<n;i++)//Q is I
		for(j=0;j<n;j++)
		{
			if(i==j) I[i][j]=1;
			else I[i][j]=0;
		}
	Sanjiaohua(A);
	while(1)
	{
		k=1
		m=n;
		while(fabs(A[m-1][m-2])<=eps) 
		{
			nmd[m-1]=A[m-1][m-1];m=m-1;
			if(m==1) {nmd[0]=A[0][0];break;}
			if(m>1) continue;
		}
		if(m==1) {nmd[0]=A[0][0];break;}
		fanchen(A[m-2][m-2],A[m-2][m-1],A[m-1][m-2],A[m-1][m-1],jieguo)
		if(m==2) break;
		while(fabs(A[m-1][m-2])<=eps) 
		{
			m=m-2;
			if(m==1) {break;}
		}
		if(m==1) {nmd[0]=A[0][0];break;}
		if(k==L) break;
		s=A[m-2][m-2]+A[m-1][m-1];
		t=A[m-2][m-2]*A[m-1][m-1]+A[m-1][m-2]*A[m-2][m-1];
		chenfa(A,A,M);//A*A
		for(i=0;i<m;i++)
			for(j=0;j<m;j++)
				M[i][j]=M[i][j]-s*A[i][j]+t*I[i][j];
		QR(M);
		
		
		
		
	}
	QR(B);
}