Программа tred2
/* Редукция Хаусхолдера действительной симметричной матрицы a[1...n][1...n].
   На выходе a заменяется ортогональной матрицей трансформации q.
   d[1...n] возвращает диагональ трехдиагональной матрицы.
   e[1...n] возвращает внедиагональные элементы, причем e[1]=0.
   Некоторые инструкции программы могут быть опущены (это указано в комментариях),
   если требуется отыскать только собственные значения, а не вектора. В
   этом случае на выходе массив a будет содержать мусор.
*/
#include <math.h>

void tred2(float **a, int n, float *d, float *e) {
  int l,k,j,i;
  float scale,hh,h,g,f;
  /* Проход по стадиям процесса редукции */
  for(i=n;i>=2;i--) {
    l=i-1; h=scale=0.;
    /* сложный процесс везде, кроме последней стадии */
    if(l>1) {
      /* вычислить шкалу */
      for(k=1;k<=l;k++) scale += fabs(a[i][k]);
      /* малая величина шкалы -> пропустить преобразование */
      if(scale==0.) e[i]=a[i][l];
      else {
       /* отмасштабировать строку и вычислить s2 в h */
        for(k=1;k<=l;k++) {
          a[i][k]/=scale; h += a[i][k]*a[i][k];
        }
       /* вычислить вектор u */
        f=a[i][l];
        g=(f>=0.?-sqrt(h):sqrt(h));
        e[i]=scale*g; h -= f*g;
       /* записать u на место i-го ряда a */
        a[i][l]=f-g;
       /* вычисление u/h, Au, p, K */
        f=0.;
        for(j=1;j<=l;j++) {
       /* следующая инструкция не нужна, если не требуются вектора,
             она содержит загрузку u/h в столбец a */
          a[j][i]=a[i][j]/h;
       /* сформировать элемент Au (в g) */
          g=0.;
          for(k=1;k<=j;k++) g += a[j][k]*a[i][k];
          for(k=j+1;k<=l;k++) g += a[k][j]*a[i][k];
       /* загрузить элемент p во временно неиспользуемую область e */
          e[j]=g/h;
       /* подготовка к формированию K */
          f += e[j]*a[i][j];
        }
       /* Сформировать K */
        hh=f/(h+h);
        for(j=1;j<=l;j++) {
       /* Сформировать q и поместить на место p (в e) */
          f=a[i][j]; e[j]=g=e[j]-hh*f;
       /* Трансформировать матрицу a */
          for(k=1;k<=j;k++) a[j][k] -= (f*e[k]+g*a[i][k]);
        }
      }
    }
    else e[i]=a[i][l];
    d[i]=h;
  }
  /* если не нужны собственные вектора, опустите следующую инструкцию */
  d[1]=0.;
  /* эту опускать не надо */
  e[1]=0.;
  /* Все содержание цикла, кроме одной инструкции, можно опустить, если не
     требуются собственные вектора */
  for(i=1;i<=n;i++) {
    l=i-1;
    /* этот блок будет пропущен при i=1 */
    if(d[i]!=0.) {
      for(j=1;j<=l;j++) {
        g=0.;
       /* формируем PQ, используя u и u/H */
        for(k=1;k<=l;k++) g += a[i][k]*a[k][j];
        for(k=1;k<=l;k++) a[k][j] -= g*a[k][i];
      }
    }
    /* эта инструкция остается */
    d[i]=a[i][i];
    /* ряд и колонка матрицы a преобразуются к единичной, для след. итерации */
    a[i][i]=0.;
    for(j=1;j<=l;j++) a[j][i]=a[i][j]=0.;
  }
} 

Программа tqli
#include <math.h>
#include "nrutil.h"    /* Здесь определяются некоторые утилиты типа выделения памяти */

/* QL-алгоритм с неявными сдвигами для определения собственных значений (и собственных
   векторов) действительной, симметричной, трехдиагональной матрицы. Эта матрица может
   быть предварительно получена с помощью программы tred2. На входе d[1...n] содержит
   диагональ исходной матрицы, на выходе - собственные значения. На входе e[1...n]
   содержит поддиагональные элементы, начиная с e[2]. На выходе массив e разрушается.
   При необходимости поиска только собственных значений в программе следует
   закомментировать или удалить инструкции, необходимые только для поиска собственных
   векторов. Если требуются собственные вектора трехдиагональной матрицы, массив
   z[1...n][1...n] необходимо инициализировать на входе единичной матрицей. Если
   требуются собственные вектора матрицы, сведенной к трехдиагональному виду с помощью
   программы tred2, в массив z требуется загрузить соответствующий выход tred2. В
   обоих случаях на выходе массив z возвращает матрицу собственных векторов, расположенных
   по столбцам.
*/

/* максимальное число итераций */
#define MAXITER 30

void tqli(float *d, float *e, int n, float **z) {
  int m,l,iter,i,k;
  float s,r,p,g,f,dd,c,b;
  /* удобнее будет перенумеровать элементы e */
  for(i=2;i<=n;i++) e[i-1]=e[i];
  e[n]=0.;
  /* главный цикл идет по строкам матрицы */
  for(l=1;l<=n;l++) {
    /* обнуляем счетчик итераций для этой строки */
    iter=0;
    /* цикл проводится, пока минор 2х2 в левом верхнем углу начиная со строки l
       не станет диагональным */
    do {
      /* найти малый поддиагональный элемент, дабы расщепить матрицу */
      for(m=l;m<=n-1;m++) {
        dd=fabs(d[m])+fabs(d[m+1]);
        if((float)(fabs(e[m]+dd)==dd) break;
      }
      /* операции проводятся, если верхний левый угол 2х2 минора еще не диагональный */
      if(m!=l) {
        /* увеличить счетчик итераций и посмотреть, не слишком ли много. Функция
           nerror завершает программу с диагностикой ошибки. */
        if(++iter>=MAXITER) nerror("Too many iterations in tqli");
        /* сформировать сдвиг */
        g=(d[l+1]-d[l])/(2.*e[l]); r=hypot(1.,g);
        /* здесь d_m - k_s */
        if(g>=0.) g+=fabs(r);
        else g-=fabs(r);
        g=d[m]-d[l]+e[l]/g;
        /* инициализация s,c,p */
        s=c=1.; p=0.;
        /* плоская ротация оригинального QL алгоритма, сопровождаемая ротациями
           Гивенса для восстановления трехдиагональной формы */
        for(i=m-1;i>=l;i--) {
          f=s*e[i]; b=c*e[i];
          e[i+1]=r=hypot(f,g);
          /* что делать при малом или нулевом знаменателе */
          if(r==0.) {d[i+1]-=p; e[m]=0.; break;}
          /* основные действия на ротации */
          s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2.c*b; d[i+1]=g+(p=s*r); g=c*r-b;
          /* Содержимое следующего ниже цикла необходимо опустить, если
             не требуются значения собственных векторов */
          for(k=1;k<=n;k++) {
            f=z[k][i+1]; z[k][i+1]=s*z[k][i]+c*f; z[k][i]=c*z[k][i]-s*f;
          }
        }
        /* безусловный переход к новой итерации при нулевом знаменателе и недоведенной
           до конца последовательности ротаций */
        if(r==0. && i>=l) continue;
        /* новые значения на диагонали и "под ней" */
        d[l]-=p; e[l]=g; e[m]=0.;
      }
    } while(m!=l);
  }
}