Программа tred2

static int[,] Input(out int n)
        {
            StreamReader file = new StreamReader("matrix.txt");
            string s = file.ReadToEnd();
            file.Close();
            string[] строка = s.Split('\n');
            string[] столбцы = строка[0].Split(' ');
            int[,] a = new int[строка.Length, столбцы.Length];
            int t = 0;
            n = 0;
            for (int i = 0; i < строка.Length; i++)
            {
                столбцы = строка.Split(' ');
                for (int j = 0; j < столбцы.Length; j++)
                {
                    t = Convert.ToInt32(столбцы[j]);
                    a[i, j] = t;
                }
            }
            return a;
        }
 
        static void Print(int[,] a)
        {
            for (int i = 0; i < a.GetLength(0); ++i, Console.WriteLine())
                for (int j = 0; j < a.GetLength(1); ++j)
                    Console.Write("{0} ", a[i, j]);
        }
        static void Main()
        {
            int n;
            int[,] myArray = Input(out n);
            Console.WriteLine("Исходный массив:");
            Print(myArray);
            tred2(myArray,myArray.GetLength,,);
            //tqli(,,myArray.GetLength,myArray);
 
            Console.ReadLine();
        }
        /* Редукция Хаусхолдера действительной симметричной матрицы a[1...n][1...n].
           На выходе a заменяется ортогональной матрицей трансформации q.
           d[1...n] возвращает диагональ трехдиагональной матрицы.
           e[1...n] возвращает внедиагональные элементы, причем e[1]=0.
           Некоторые инструкции программы могут быть опущены (это указано в комментариях),
           если требуется отыскать только собственные значения, а не вектора. В
           этом случае на выходе массив a будет содержать мусор.
        */
        // #include <math.h>
        static 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 = 0f;
                /* сложный процесс везде, кроме последней стадии */
                if (l > 1)
                {
                    /* вычислить шкалу */
                    for (k = 1; k <= 1; k++)
                        scale += Math.Abs(a[i, k]);
                    /* малая величина шкалы -> пропустить преобразование */
                    if (scale == 0)
                        e = 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 = ((float)(f >= 0f ? -Math.Sqrt(h) : Math.Sqrt(h)));
                        e = scale * g;
                        h -= f * g;
                        /* записать u на место i-го ряда a */
                        a[i, l] = f - g;
                        /* вычисление u/h, Au, p, K */
                        f = 0f;
                        for (j = 1; j <= 1; j++)
                        {
                            /* следующая инструкция не нужна, если не требуются вектора,
             она содержит загрузку u/h в столбец a */
                            //a[j]=a[j]/h;
                            /* сформировать элемент Au (в g) */
                            g = 0f;
                            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 = a[i, l];
                d = h;
            }
            /* если не нужны собственные вектора, опустите следующую инструкцию */
            //d[1]=0.;
            /* эту опускать не надо */
            e[1] = 0f;
            /* Все содержание цикла, кроме одной инструкции, можно опустить, если не
    требуются собственные вектора */
            for (i = 1; i <= n; i++)
            {
                //l=i-1;
                ///* этот блок будет пропущен при i=1 */
                //if(d!=0.)
                //{
                //    for(j=1;j<=l;j++)
                //    {
                //        g=0.;
                //        /* формируем PQ, используя u и u/H */
                //        for(k=1;k<=l;k++)
                //            g+=a[k]*a[k][j];
                //        for(k=1;k<=l;k++)
                //            a[k][j]-=g*a[k];
                //    }
                //}
                /* эта инструкция остается */
                d = a[i, i];
                /* ряд и колонка матрицы a преобразуются к единичной, для след. итерации */
                a[i, i] = 1;
                for (j = 1; j <= 1; j++)
                    a[j, i] = a[i, j] = 0f;
            }
 
 
        }