using System;
using System.IO;
namespace prog
{
class program
{
static float[,] Input(out float n)
{
StreamReader file = new StreamReader("matrix.txt");
string s = file.ReadToEnd();
file.Close();
string[] строка = s.Split('\n');
string[] столбцы = строка[0].Split(' ');
float[,] a = new float[строка.Length, столбцы.Length];
float t = 0;
n = 0;
for (int i = 0; i < строка.Length; i++)
{
столбцы = строка[i].Split(' ');
for (int j = 0; j < столбцы.Length; j++)
{
t = float.Parse(столбцы[j]);
a[i, j] = t;
}
}
return a;
}
static void Print(float[,] 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,i;
float p;
float[,] myArray = Input(out p);
Console.WriteLine("Исходный массив:");
Print(myArray);
n = myArray.GetLength(0);
float[] d = new float[n+ 1];
float[] e = new float[n + 1];
for (i = 0; i < n + 1; i++)
d[i]=e[i]=0;
tred2(myArray, n, d, e);
Print(myArray);
tqli(d,e,n,myArray);
Console.ReadLine();
}
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[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=((float)(f>=0f?-Math.Sqrt(h):Math.Sqrt(h)));
e[i]=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][i]=a[i][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[i]=a[i,l];
d[i]=h;
}
/* если не нужны собственные вектора, опустите следующую инструкцию */
//d[1]=0.;
/* эту опускать не надо */
e[1]=0f;
/* Все содержание цикла, кроме одной инструкции, можно опустить, если не
требуются собственные вектора */
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]=1;
for(j=1;j<=1;j++)
a[j,i]=a[i,j]=0f;
}
}
static void tqli(float []d, float []e, int n, float [,]z)
{
int MAXITER =30;
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]=0f;
/* главный цикл идет по строкам матрицы */
for(l=1;l<=n;l++)
{
/* обнуляем счетчик итераций для этой строки */
iter=0;
/* цикл проводится, пока минор 2х2 в левом верхнем углу начиная со строки l
не станет диагональным */
do {
/* найти малый поддиагональный элемент, дабы расщепить матрицу */
for(m=l;m<=n-1;m++)
{
dd=Math.Abs(d[m])+Math.Abs(d[m+1]);
if((float)(Math.Abs(e[m]+dd))==dd) break;
}
/* операции проводятся, если верхний левый угол 2х2 минора еще не диагональный */
if(m!=l) {
/* увеличить счетчик итераций и посмотреть, не слишком ли много. */
if(++iter>=MAXITER)
Console.WriteLine("Too many iterations in program");
/* сформировать сдвиг */
g=(d[l+1]-d[l])/(2f*e[l]); r=(float)hypot(1f,g);
/* здесь d_m - k_s */
if(g>=0f) g+=Math.Abs(r);
else g-=Math.Abs(r);
g=d[m]-d[l]+e[l]/g;
/* инициализация s,c,p */
s=c=1f; p=0f;
/* плоская ротация оригинального QL алгоритма, сопровождаемая ротациями
Гивенса для восстановления трехдиагональной формы */
for(i=m-1;i>=l;i--)
{
f=s*e[i]; b=c*e[i];
e[i+1]=r=(float)hypot(f,g);
/* что делать при малом или нулевом знаменателе */
if(r==0f) {d[i+1]-=p; e[m]=0f; break;}
/* основные действия на ротации */
s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2f*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==0f && i>=l)
continue;
/* новые значения на диагонали и "под ней" */
d[l]-=p; e[l]=g; e[m]=0f;
}
} while(m!=l);
}
}
static double hypot(double x, double y)
{
double a = Math.Sqrt(x * x + y * y);
return a;
}
}
}
Посмотрел код, полагаю, что в цикле после комментария