Программа tred2 - C++ - Ответ 1923885

/* Редукция Хаусхолдера действительной симметричной матрицы 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.;
  }
}