博主系某985信息与计算科学大三,刚调试成功,分享给正在学习数值计算的伙伴们,有问题还请大家批评指正。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 100
struct vec
{
float v[N];
float beta;
};
//struct vec save;
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,k,n;
float x[N];
struct vec left;
//读取需要规范化的向量x
printf("Please enter the length of your vector:\n");
scanf("%d",&n);
printf("Please enter the elements of your vector:\n");
for (i=1;i<=n;i++)
{
scanf("%f",&x[i]);
}
//检查
for (i=1;i<=n;i++)
{
printf("%f",x[i]);
}
//读取数据完毕
left=house(x,n);
for (j=1;j<=n;j++)
{
printf("%f \n",left.v[j]);
}
printf("%f",left.beta);
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=1;k<=n;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;
}