不多说,直接上代码。
#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;
}