Программная реализация метода Крылова вычисления собственных значений

#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);

}