qr_householder.cpp

用householder变换求QR分解

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef EPSLION
#define EPSLION 1.0e-10
#define M 84
#define N 84
#endif



//输入1-1
void Input_1(double A[M][N])
{
	int i,j;
	for (i=0;i<M;i++){
	for(j=0;j<N;j++){
			A[i][j]=0;}}
	for (i=0;i<M;i++){
		A[i][i]=6;}
	for (i=0;i<N-1;i++){
		A[i][i+1]=1;}
	for (i=1;i<N;i++){
		A[i][i-1]=8;}
}

//输入1-2
void Input_2(double A[M][N])
{
	int i,j;
	for (i=0;i<M;i++){
	for(j=0;j<N;j++){
			A[i][j]=0;}}
	for (i=0;i<M;i++){
		A[i][i]=10;}
	for (i=0;i<N-1;i++){
		A[i][i+1]=1;}
	for (i=1;i<N;i++){
		A[i][i-1]=1;}
}
//输入1-3
void Input_3(double A[M][N])
{
	int i,j;
	double temp;
	double a=1.0;
	for (i=0;i<M;i++){
		for (j=0;j<N;j++){
			temp=i+j+1;
			A[i][j]=a/temp;
		}
	}
}

//householder transformation
void householder (double A[M][N],double y[M])
{
	int i,j,l;
	static double aa,ki,u[M],alpha;
	for (i=0;i<N;i++)
	{
		aa=0.0;
		for (l=i;l<M;l++) aa=aa+A[l][i]*A[l][i];
		aa=sqrt(aa);
		if(aa<EPSLION)
		{
			printf("The rank of overdetermined matrix A is less than COLUMN\n");
			printf("Please use singluar value decomposition if matrix A for\n");
			printf("minimal least square solution.\n ");
			printf("Press any key to exit!\n");
			exit(1);
		}
		if (A[i][i]>0.0) ki=-aa;
		else ki=aa;
		alpha=1.0/(ki*(ki-A[i][i]));
		u[i]=A[i][i]-ki;
		A[i][i]=ki;
		for(l=i+1;l<M;l++) {u[l]=A[l][i];A[l][i]=0.0;}
		for(j=i+1;j<N;j++)
		{
			aa=0.0;
			for(l=i;l<M;l++) aa=aa+u[l]*A[l][j];
			aa=alpha*aa;
			for(l=i;l<M;l++) A[l][j]=A[l][j]-aa*u[l];
		}
		aa=0.0;
		for(l=i;l<M;l++) aa=aa+u[l]*y[l];
		aa=alpha*aa;
		for(l=i;l<M;l++) y[l]=y[l]-aa*u[l];
	}
	printf("Matrix QA,and Qy:\n");
	for(i=0;i<M;i++)
	{
		for(j=0;j<N;j++){ printf("%14.8f",A[i][j]);}
		printf("   |%14.8f\n",y[i]);	
	}
}



void pofitho (double y[M],double b[N])
{
	int i,j;
	double A[M][N];
	printf("Overdetermined equations:\n");
	Input_1(A);
	for (i=0;i<M;i++)
	{
		for (j=0;j<N;j++) printf("%10.5f",A[i][j]);
		printf("  |%10.5f \n",y[i]);
	}
	householder(A,y);
	for (i=N-1;i>=0;i--)
	{
		for (j=i+1;j<N;j++) y[i]=y[i]-A[i][j]*b[j];
		b[i]=y[i]/A[i][i];
	}
}


void main()
{
	int i;
	double b[N];
	double y[M];
	printf("Polynomial fit, Householder transformation:\n");
	pofitho(y,b);
	for (i=0;i<N;i++)
		printf("a(%2d)=%14.8e\n",i,b[i]);
}