Собственные значения матрицы

#include <conio.h>
#include "matrix.h"
#include "equation.h"
cvector Krylov_eigenval(cmatrix &x)
{
if(x.getn()!=x.getm())
 throw xmsg("Для неквадратных матриц собственные 
значения не определены");
/*Согласно данному методу, нам необходимо составить 
последовательность линейно независимых векторов по 
специальному закону. Для удобства использования мы 
будем хранить их как вектор-столбцы, то есть как
последовательность одностолбцовых матриц:*/
cmatrix *B; /*указатель под последовательность 
векторов*/
//пытаемся выделить память
try
 {
 B=new cmatrix[x.getm()+1];
}
/*перехват этого исключения происходит при не-
хватке памяти*/
catch(xalloc)
{
/*не пытаясь исправить эту ситуацию, выбрасываем, в 
свою очередь, стандартное исключение с диагностиче-
ским текстом */
 throw xmsg("Не хватает памяти");
}
/*заполняем массив В векторами (одностолбцовыми 
матрицами) заданного размера */
for(long i=0;i<x.getm()+1;i++)
 B[i]=cmatrix(x.getm(),1);
/*создаём пока пустую матрицу, хранящую в своих 
столбцах строящуюся последовательность векторов*/
cmatrix test(x.getm(),x.getm());
do
{
 //формируем случайный вектор
 for(long i=0;i<B[0].getm();i++)
 B[0][i][0]=random(x.getm()+1);
 /*действуя на вектор линейным оператором, за-
данным матрицей, находим следующий вектор и т.д.*/
 for(long i=0;i<x.getm();i++)
 B[i+1]=x*B[i];
 /*переносим вектора из массива в соответствую-
щие столбцы матрицы*/
 for(long i=0;i<test.getm();i++)
 for(long j=0;j<test.getm();j++)
 test[j][i]=B[i][j][0];
 randomize();/*переустанавливаем генератор слу-
чайных чисел*/
/*процесс построения повторяем до тех пор, пока не 
выполнится условие линейной независимости векторов 
- неравенство нулю детерминанта матрицы, составлен-
ной из построенных векторов*/
}while(det(test)==(complex)0);
/*разлагаем последний вектор по предыдущим (в си-
лу их линейной независимости) как по базису и нахо-
дим коэффициенты разложения*/
cmatrix coeffs=SLAE_Orto(test,B[x.getm()]);
delete[] B;//освобождаем память из-под массива
 //формируем полином из вычисленных коэффициентов
 cpolynom p=x.getm()+1;
 p[x.getm()]=1;
 for(long i=0;i<coeffs.getm();i++)
 
p[i]=-coeffs[i][0];
/*находим корни полинома (собственные значения) мо-
дифицированным методом Ньютона*/
return newton(p);
}

float Matr(float A[][N], float norm) {
    float w0[100], w[100], summ = 0, w0norm[100], e, d, d0;
    int i, j, k;
    w0[0] = 1;
    for (i = 1; i < N; i++)
        w0[i] = 0;
    do {
        for (i = 0; i < N; i++)
            summ = summ + w0[i] * w0[i];
        d0 = sqrt(summ);
        for (i = 0; i < N; i++)
            w0norm[i] = w0[i] / d0;
        for (i = 0; i < N; i++) {
            w[i] = 0;
            for (j = 0; j < N; j++)
                w[i] = w[i] + A[i][j] * w0norm[j];
        }
        summ = 0;
        for (i = 0; i < N; i++)
            summ = summ + w[i] * w[i];
        d = sqrt(summ);
        e = fabs(d - d0);
        for (i = 0; i < N; i++)
            w0[i] = w[i];
        summ = 0;
    } while (e > 0.001);
 
    std::cout << "\nL = " << d << std::endl;
    norm = sqrt(d);
    return norm;
}

float L, E;
int i, j;
 
void matrix_m_vector(float a[N][N], float v[N], float res[N]) {
    for (i = 0; i < N; i++) {
        res[i] = 0;
        for (j = 0; j < N; j++) {
            res[i] += a[i][j] * v[j];
        }
    }
}
 
void vector_m_vector(float v1[N], float v2[N], float res[N]) {
    for (i = 0; i < N; i++) {
        *res += v1[i] * v2[i];
    }
}
 
float NORM(float A[N][N], float norm) {
 
    float *Xn = new float[N];
    float *Xn1 = new float[N];
    
    cout << "\ninsert priblijenie\n";
    for (i = 0; i < N; i++)
        cin >> Xn[i];
    cout << "\nE = ";
    cin >> E;
 
    while (1) {
        matrix_m_vector(A, Xn, Xn1);
        float num = 0, den = 0;
        vector_m_vector(Xn1, Xn, &num);
        vector_m_vector(Xn, Xn, &den);
        float Lprev = L;
        L = num / den;
        if ((fabs(L - Lprev)) < E) break;
        memcpy(Xn, Xn1, sizeof(Xn));
    }
    cout << "max = " << L << endl;
    return L;
}