C语言实现矩阵QR分解（利用householder变换）

不多说，直接上代码。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 100
struct vec
{
	float v[N];
	float beta;
};
//struct vec save;
float A[N][N];
struct vec house(float [],int n);
float two_norm(float x[],int begin,int end);
float max(float x[],int n);

int main(void)
{
	int i,j,t,k,m,n;
	float d[N],sum[N][N],temp[N],tem[N][N],temp_v[N][N];
	struct vec left;
	
	printf("Please enter the row and col of matrix A:\n");
	scanf("%d %d",&m,&n);
	printf("Please enter the elements of the matrix A:\n");
	for (i=1;i<=m;i++)
	for (j=1;j<=n;j++)
	scanf("%f",&A[i][j]);
	
	for (j=1;j<=n;j++)
	{
		for (i=0;i<=N-1;i++)
		{
			for (k=0;k<=N-1;k++)
			{
				sum[i][k]=0.0;
			}
		}
		if (j<m)
		{
			for (i=j;i<=m;i++)
			temp[i-j+1]=A[i][j];
			
			left=house(temp,m-j+1);
			d[j]=left.beta;
			for (i=1;i<=m-j+1;i++)
			{
				for (k=1;k<=m-j+1;k++)
				{
					temp_v[i][k]=left.v[i]*left.v[k];
					if (i==k)
					{
						tem[i][k]=1.0-d[j]*temp_v[i][k];
					}
					else
					tem[i][k]=-d[j]*temp_v[i][k];
				}
			}
			//存储矩阵I-beta*vv^T
			for (i=j;i<=m;i++)
			{
				for (k=j;k<=n;k++)
				{
					for (t=j;t<=m;t++)
					{
						sum[i-j+1][k-j+1]+=tem[i-j+1][t-j+1]*A[t][k];
					}
				}
			} 
			
			for (i=j;i<=m;i++)
			{   
				for (k=j;k<=n;k++)
				{
					A[i][k]=sum[i-j+1][k-j+1];
				}
			}

			for (i=j+1;i<=m;i++)
			A[i][j]=left.v[i-j+1];
		}
	}
	
	for (i=1;i<=m;i++)
	{
		for (k=1;k<=n;k++)
		{
			printf("%f ",A[i][k]);
		}
		printf("\n");
	}
	
	return 0;
}


//householder变换的实现
struct vec house(float x[],int n)
{
	int i,j,k;
	float m,sigma,alpha;
	struct vec save;
	m=max(x,n);
	printf(" %f\n",m);
	for (i=1;i<=n;i++)
	{
		x[i]=x[i]/m;
	}
	sigma=two_norm(x,2,n);
	sigma=pow(sigma,2);
	for (j=2;j<=n;j++)
	{
		save.v[j]=x[j];
	}
	if (sigma==0)
	{
		save.beta=0.0;
	}
	else
	{
		alpha=sqrt(pow(x[1],2)+sigma);
		if (x[1]<=0)
		{
			save.v[1]=x[1]-sigma;
		}
		else 
		{
			save.v[1]=-(sigma/(x[1]+alpha));
		}
		save.beta=2*pow(save.v[1],2)/(sigma+pow(save.v[1],2));
	}
	for (k=n;k>=1;k--)
    {
		save.v[k]=save.v[k]/save.v[1];
    }
	
	return save;
} 

float two_norm(float x[],int begin,int end)
{
	int i;
	float val=0.0;
	
	for (i=begin;i<=end;i++)
	{
		val+=pow(x[i],2);
	}
	val=sqrt(val);
	
	return val;
}

float max(float x[],int n)
{
	int i,j;
	float max_val;
	max_val=x[1];
	
	for (i=1;i<=n;i++)
	{
		if (max_val<x[i])
		{
			max_val=x[i];
		}
	}
	
	return max_val;
}