用householder 变换化一般矩阵为hesssenberg矩阵的C程序
#include "stdio.h"#include "math.h"
#define n 10
int sign(double x)
{
int y;
if(x>0) y=1;
else if(x<0) y=-1;
else y=0;
return(y);
}
void matrix_print(double A[n][n])
{
int i,j;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
printf("%8.4lf\t",A[i][j]);
printf("\n");
}
}
double matrix_norm(double a[n])
{
int i;
double d=0;
for(i=0;i<n;i++)
d+=a[i]*a[i];
return(sqrt(d));
}
void matrix_multiply(double A[n][n],double B[n][n],double C[n][n])
{
int i,j,t;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
C[i][j]=0;
for(t=0;t<n;t++)
C[i][j]+=A[i][t]*B[t][j];
}
}
void matrix_copy(double A[n][n],double B[n][n])
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
A[i][j]=B[i][j];
}
void householder_trans(double A[n][n],int k)
{
double a[n];
double H[n][n];
double d=0;
int i,j;
for(i=0;i<=k;i++)
a[i]=0;
for(i=k+1;i<n;i++)
a[i]=A[i][k];
d=matrix_norm(a);
if(fabs(d)!=fabs(A[k+1][k]))
a[k+1]+=sign(A[k+1][k])*d;
d=matrix_norm(a);
for(i=0;i<n;i++)
a[i]=a[i]/d;
printf("%8.4lf\t",d);
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
H[i][j]=-2*a[i]*a[j];
H[i][i]++;
}
double temp[n][n];
printf("H:\n");
matrix_print(H);
matrix_multiply(H,A,temp);
matrix_copy(A,temp);
matrix_multiply(A,H,temp);
matrix_copy(A,temp);
}
void matrix_input(double A[n][n])
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
printf("A[%d][%d]=",i,j);
scanf("%lf",&A[i][j]);
}
}
void main()
{
int i;
double A[n][n];
matrix_input(A);
printf("A:\n");
matrix_print(A);
for(i=0;i<n-2;i++)
householder_trans(A,i);
printf("hessenberg is:\n");
matrix_print(A);
}
