约化对称矩阵为三对角对称矩阵

//约化对称矩阵为三对角对称矩阵

//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵

//a-长度为n*n的数组，存放n阶实对称矩阵

//n-矩阵的阶数

//q-长度为n*n的数组，返回时存放Householder变换矩阵

//b-长度为n的数组，返回时存放三对角阵的主对角线元素

//c-长度为n的数组，返回时前n-1个元素存放次对角线元素

void eastrq(double a[],int n,double q[],double b[],double c[]);

//////////////////////////////////////////////////////////////

//求实对称三对角对称矩阵的全部特征值及特征向量

//利用变型QR方法计算实对称三对角矩阵全部特征值及特征向量

//n-矩阵的阶数

//b-长度为n的数组，返回时存放三对角阵的主对角线元素

//c-长度为n的数组，返回时前n-1个元素存放次对角线元素

//q-长度为n*n的数组，若存放单位矩阵，则返回实对称三对角矩阵的特征向量组

// 若存放Householder变换矩阵，则返回实对称矩阵A的特征向量组

//a-长度为n*n的数组，存放n阶实对称矩阵

int ebstq(int n,double b[],double c[],double q[],double eps,int l);

//////////////////////////////////////////////////////////////

//约化实矩阵为赫申伯格(Hessen berg)矩阵

//利用初等相似变换将n阶实矩阵约化为上H矩阵

//a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵

//n-矩阵的阶数

void echbg(double a[],int n);

//////////////////////////////////////////////////////////////

//求赫申伯格(Hessen berg)矩阵的全部特征值

//返回值小于0表示超过迭代jt次仍未达到精度要求

//返回值大于0表示正常返回

//利用带原点位移的双重步QR方法求上H矩阵的全部特征值

//a-长度为n*n的数组，存放上H矩阵

//n-矩阵的阶数

//u-长度为n的数组，返回n个特征值的实部

//v-长度为n的数组，返回n个特征值的虚部

//eps-控制精度要求

//jt-整型变量，控制最大迭代次数

int edqr(double a[],int n,double u[],double v[],double eps,int jt);

//////////////////////////////////////////////////////////////

//求实对称矩阵的特征值及特征向量的雅格比法

//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量

//返回值小于0表示超过迭代jt次仍未达到精度要求

//返回值大于0表示正常返回

//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值

//n-矩阵的阶数

//u-长度为n*n的数组，返回特征向量(按列存储)

//eps-控制精度要求

//jt-整型变量，控制最大迭代次数

int eejcb(double a[],int n,double v[],double eps,int jt);

//////////////////////////////////////////////////////////////

选自<<徐世良数值计算程序集(C)>>

每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。今天都给贴出来了

#include "stdio.h"

#include "math.h"

//约化对称矩阵为三对角对称矩阵

//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵

//a-长度为n*n的数组，存放n阶实对称矩阵

//n-矩阵的阶数

//q-长度为n*n的数组，返回时存放Householder变换矩阵

//b-长度为n的数组，返回时存放三对角阵的主对角线元素

//c-长度为n的数组，返回时前n-1个元素存放次对角线元素

void eastrq(double a[],int n,double q[],double b[],double c[])

{

int i,j,k,u,v;

double h,f,g,h2;

for (i=0; i<=n-1; i++)

{

for (j=0; j<=n-1; j++)

{

u=i*n+j; q[u]=a[u];

}

}

for (i=n-1; i>=1; i--)

{

h=0.0;

if (i>1)

{

for (k=0; k<=i-1; k++)

{

u=i*n+k;

h=h+q[u]*q[u];

}

}

if (h+1.0==1.0)

{

c[i-1]=0.0;

if (i==1)

{

c[i-1]=q[i*n+i-1];

}

b[i]=0.0;

}

else

{

c[i-1]=sqrt(h);

u=i*n+i-1;

if (q[u]>0.0)

{

c[i-1]=-c[i-1];

}

h=h-q[u]*c[i-1];

q[u]=q[u]-c[i-1];

f=0.0;

for (j=0; j<=i-1; j++) {

q[j*n+i]=q[i*n+j]/h;

g=0.0;

for (k=0; k<=j; k++) {

g=g+q[j*n+k]*q[i*n+k];

}

if (j+1<=i-1)

{

for (k=j+1; k<=i-1; k++)

{

g=g+q[k*n+j]*q[i*n+k];

}

}

c[j-1]=g/h;

f=f+g*q[j*n+i];

}

h2=f/(h+h);

for (j=0; j<=i-1; j++) {

f=q[i*n+j];