或 http://topic.csdn.net/t/20031205/13/2531097.html
http://www.koders.com/c/fidAFD2E5385C63D0DEB3D455AA065A7F0D9A2064AC.aspx
  
#include   "math.h"   
  #include   "stdio.h"   
    
  #define   SIGN(a,b)   ((b)   >   0   ?   fabs(a)   :   -fabs(a))   
  #define   SWAP(g,h)   {y=(g);(g)=(h);(h)=y;}   
  #define   RADIX   2.0   
    
  void   balanc(double   **a,   int   n)   
  {   
  int   last,j,i;   
  double   s,r,g,f,c,sqrdx;   
  sqrdx=RADIX*RADIX;   
  last=0;   
  while   (last   ==   0)     
  {   
  last=1;   
  for   (i=1;i<=n;i++)     
  {   
  r=c=0.0;   
  for   (j=1;j<=n;j++)   
  if   (j   !=   i)     
  {   
  c   +=   fabs(a[j][i]);   
  r   +=   fabs(a[i][j]);   
  }   
  if   (c   &&   r)     
  {   
  g=r/RADIX;   
  f=1.0;   
  s=c+r;   
  while   (c<g)     
  {   
  f   *=   RADIX;   
  c   *=   sqrdx;   
  }   
  g=r*RADIX;   
  while   (c>g)   
  {   
  f   /=   RADIX;   
  c   /=   sqrdx;   
  }   
  if   ((c+r)/f   <   0.95*s)     
  {   
  last=0;   
  g=1.0/f;   
  for   (j=1;j<=n;j++)     
  a[i][j]   *=   g;   
  for   (j=1;j<=n;j++)   
  a[j][i]   *=   f;   
  }   
  }   
  }   
  }   
  }   
    
  void   elmhes(double   **a,int   n)   
  {   
  int   m,j,i;   
  double   y,x;   
    
  for   (m=2;m<n;m++)   {   
  x=0.0;   
  i=m;   
  for   (j=m;j<=n;j++)   {   
  if   (fabs(a[j][m-1])   >   fabs(x))   {   
  x=a[j][m-1];   
  i=j;   
  }   
  }   
  if   (i   !=   m)   {   
  for   (j=m-1;j<=n;j++)   SWAP(a[i][j],a[m][j])   
  for   (j=1;j<=n;j++)   SWAP(a[j][i],a[j][m])   
  }   
  if   (x)   {   
  for   (i=m+1;i<=n;i++)   {   
  if   (y=a[i][m-1])   {   
  y   /=   x;   
  a[i][m-1]=y;   
  for   (j=m;j<=n;j++)   
  a[i][j]   -=   y*a[m][j];   
  for   (j=1;j<=n;j++)   
  a[j][m]   +=   y*a[j][i];   
  }   
  }   
  }   
  }   
  }   
    
  void   hqr(double   **a,int   n,double   *wr,double   *wi)   
  {   
  int   nn,m,l,k,j,its,i,mmin;   
  double   z,y,x,w,v,u,t,s,r,q,p,anorm;   
    
  anorm=fabs(a[1][1]);   
  for   (i=2;i<=n;i++)   
  for   (j=(i-1);j<=n;j++)   
  anorm   +=   fabs(a[i][j]);   
  nn=n;   
  t=0.0;   
  while   (nn   >=   1)   {   
  its=0;   
  do   {   
  for   (l=nn;l>=2;l--)   {   
  s=fabs(a[l-1][l-1])+fabs(a[l][l]);   
  if   (s   ==   0.0)   s=anorm;   
  if   (fabs(a[l][l-1])   +   s   ==   s)   break;   
  }   
  x=a[nn][nn];   
  if   (l   ==   nn)   {   
  wr[nn]=x+t;   
  wi[nn--]=0.0;   
  }   else   {   
  y=a[nn-1][nn-1];   
  w=a[nn][nn-1]*a[nn-1][nn];   
  if   (l   ==   (nn-1))   {   
  p=0.5*(y-x);   
  q=p*p+w;   
  z=sqrt(fabs(q));   
  x   +=   t;   
  if   (q   >=   0.0)   {   
  z=p+SIGN(z,p);   
  wr[nn-1]=wr[nn]=x+z;   
  if   (z)   wr[nn]=x-w/z;   
  wi[nn-1]=wi[nn]=0.0;   
  }   else   {   
  wr[nn-1]=wr[nn]=x+p;   
  wi[nn-1]=   -(wi[nn]=z);   
  }   
  nn   -=   2;   
  }   else   {   
  if   (its   ==   30)   printf("Too   many   iterations   in   HQR");   
  if   (its   ==   10   ||   its   ==   20)   {   
  t   +=   x;   
  for   (i=1;i<=nn;i++)   a[i][i]   -=   x;   
  s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);   
  y=x=0.75*s;   
  w   =   -0.4375*s*s;   
  }   
  ++its;   
  for   (m=(nn-2);m>=l;m--)   {   
  z=a[m][m];   
  r=x-z;   
  s=y-z;   
  p=(r*s-w)/a[m+1][m]+a[m][m+1];   
  q=a[m+1][m+1]-z-r-s;   
  r=a[m+2][m+1];   
  s=fabs(p)+fabs(q)+fabs(r);   
  p   /=   s;   
  q   /=   s;   
  r   /=   s;   
  if   (m   ==   l)   break;   
  u=fabs(a[m][m-1])*(fabs(q)+fabs(r));   
  v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));   
  if   (u+v   ==   v)   break;   
  }   
  for   (i=m+2;i<=nn;i++)   {   
  a[i][i-2]=0.0;   
  if     (i   !=   (m+2))   a[i][i-3]=0.0;   
  }   
  for   (k=m;k<=nn-1;k++)   {   
  if   (k   !=   m)   {   
  p=a[k][k-1];   
  q=a[k+1][k-1];   
  r=0.0;   
  if   (k   !=   (nn-1))   r=a[k+2][k-1];   
  if   (x=fabs(p)+fabs(q)+fabs(r))   {   
  p   /=   x;   
  q   /=   x;   
  r   /=   x;   
  }   
  }   
  if   (s=SIGN(sqrt(p*p+q*q+r*r),p))   {   
  if   (k   ==   m)   {   
  if   (l   !=   m)   
  a[k][k-1]   =   -a[k][k-1];   
  }   else   
  a[k][k-1]   =   -s*x;   
  p   +=   s;   
  x=p/s;   
  y=q/s;   
  z=r/s;   
  q   /=   p;   
  r   /=   p;   
  for   (j=k;j<=nn;j++)   {   
  p=a[k][j]+q*a[k+1][j];   
  if   (k   !=   (nn-1))   {   
  p   +=   r*a[k+2][j];   
  a[k+2][j]   -=   p*z;   
  }   
  a[k+1][j]   -=   p*y;   
  a[k][j]   -=   p*x;   
  }   
  mmin   =   nn<k+3   ?   nn   :   k+3;   
  for   (i=l;i<=mmin;i++)   {   
  p=x*a[i][k]+y*a[i][k+1];   
  if   (k   !=   (nn-1))   {   
  p   +=   z*a[i][k+2];   
  a[i][k+2]   -=   p*r;   
  }   
  a[i][k+1]   -=   p*q;   
  a[i][k]   -=   p;   
  }   
  }   
  }   
  }   
  }   
  }   while   (l   <   nn-1);   
  }   
  }   
    
  void   eig(double   **a,double   *wr,double   *wi,int   n)   
  {   
  balanc(a,n);   
  elmhes(a,n);   
  hqr(a,n,wr,wi);   
  }   
    
  void   solve(double   *aa,double   *wr,double   *wi,int   n)   
  //aa:1*(n+1),x:1:n   
  {   
  int   i,j;   
  double   **a;   
  a=new   double*[n+1];   
  for(i=0;i<n+1;i++)   
  a[i]=new   double[n+1];   
  for(i=0;i<n+1;i++)   
  for(j=0;j<n+1;j++)   
  a[i][j]=0;   
  for(i=1;i<n+1;i++)   
  a[1][i]=-aa[i]/aa[0];   
  for(i=1;i<n;i++)   
  a[i+1][i]=1;   
  eig(a,wr,wi,n);   
  }   
    
  int   main1(int   argc,   char*   argv[])   
  {   
  int   i,j,n;   
  double   **a,*wr,*wi;   
 FILE   *fp=fopen("d:\\hqr.txt","r");   
  fscanf(fp,"%d",&n);   
  a=new   double*[n+1];   
  wr=new   double[n+1];   
  wi=new   double[n+1];   
  for(i=0;i<n+1;i++)   
  a[i]=new   double[n+1];   
  for(i=1;i<n+1;i++)   
  for(j=1;j<n+1;j++)   
  fscanf(fp,"%lf",a[i]+j);   
  fclose(fp);   
    
  printf("原矩陣\n");   
  for(i=1;i<n+1;i++)   
  {   
  for(j=1;j<n+1;j++)   
  printf("%f\t",a[i][j]);   
  printf("\n");   
  }   
  printf("\n");   
    
  eig(a,wr,wi,n);   
    
  printf("特征值\n");   
  for(i=1;i<n+1;i++)   
  printf("%f\t+i*%f\n",wr[i],wi[i]);   
  return   0;   
  }   
    
  void   main()   
  {   
  int   i;   
  const   int   n=6;   
  double   A=1e-2,B=1e2,R=1e-2;   
  double   aa[n+1],wr[n+1],wi[n+1];   
  //aa:方程系數，從高次到低次   
  aa[0]=5.814e47;   
  aa[1]=4.41037e45+3.42e43*B+6.12e42*A*R;   
  aa[2]=2.8072e38-4.64249e40*A+2.59434e41*B+4.6425e40*A*R;   
  aa[3]=-4.3949e31-5.88422e33*A+1.63691e34*B+2.95501e33*A*R;   
  aa[4]=-2.80726e24-3.25852e24*A+7.41539e21*B+1.62862e24*A*R;   
  aa[5]=1.54719e15-2.13108e12*A+7.10361e11*A*R;   
  aa[6]=-674.843;   
  solve(aa,wr,wi,n);   
  printf("方程的解\n");   
  for(i=1;i<n+1;i++)   
  printf("%18.16f\t+i*%18.16f\n",wr[i],wi[i]);   
  }