Метод Ланцоша для симметричных СЛАУ. C#

Пытаюсь реализовать данный метод http://window.edu.ru/resource/183/41183/files/nstu144.pdf, стр 52. Вроде бы сделал все как в алгоритме, но результат выдает неверный. Не могу найти ошибку, нужна помощь!

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;

namespace SymmetricMatrix
{
    public partial class Form1 : Form
    {
        int n;
        List<List<double>> matrix;
        List<double> b;
        public Form1()
        {
            InitializeComponent();
            matrix = new List<List<double>>();
            b = new List<double>();
            n = 3;
            Init();
            LanczosMethod();

        }

        private void Init()
        {
            Random rnd = new Random();
            for (int i = 0; i < n; i++)
            {
                matrix.Add(new List<double>());
                for (int j = 0; j < n; j++)
                    matrix[i].Add(0);
            }
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    if (i == j)
                        matrix[i][j] = rnd.Next(-10, 10)/* + rnd.NextDouble()*/;
                    else if (i < j)
                    {
                        double randomValue = rnd.Next(-10, 10)/* + rnd.NextDouble()*/;
                        matrix[i][j] = randomValue;
                        matrix[j][i] = randomValue;
                    }
                }
            }

            for (int i = 0; i < n; i++)
                b.Add(rnd.Next(-10, 10));
        }

        private List<double> MatrToVec(List<List<double>> matr, List<double> vec)
        {
            List<double> result = new List<double>();
            for (int i = 0; i < n; i++)
                result.Add(0);
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    result[i] += matr[i][j] * vec[j];
            }
            return result;
        }

        private double Norma(List<double> vec)
        {
            double res = 0;
            for (int i = 0; i < n; i++)
                res += vec[i] * vec[i];
            return Math.Sqrt(res);
        }

        private List<double> Division(List<double> vec, double x)
        {
            List<double> result = new List<double>();
            for (int i = 0; i < n; i++)
                result.Add(0);
            for (int i = 0; i < n; i++)
                result[i] = vec[i] / x;
            return result;
        }

        private double ScalarMult(List<double> vec1, List<double> vec2)
        {
            double res = 0;
            for (int i = 0; i < n; i++)
                res += vec1[i] * vec2[i];
            return res;
        }

        private List<double> Sub(List<double> vec1, List<double> vec2)
        {
            List<double> result = new List<double>();
            for (int i = 0; i < n; i++)
                result.Add(0);
            for (int i = 0; i < n; i++)
                result[i] = vec1[i] - vec2[i];
            return result;
        }

        private List<double> Sum(List<double> vec1, List<double> vec2)
        {
            List<double> result = new List<double>();
            for (int i = 0; i < n; i++)
                result.Add(0);
            for (int i = 0; i < n; i++)
                result[i] = vec1[i] + vec2[i];
            return result;
        }

        private List<double> Mult(List<double> vec, double x)
        {
            List<double> result = new List<double>();
            for (int i = 0; i < n; i++)
                result.Add(0);
            for (int i = 0; i < n; i++)
               result[i] = vec[i] * x;
            return result;
        }

        private List<double> fill(List<double> vec, double x)
        {
            for (int i = 0; i < n; i++)
                vec.Add(x);
            return vec;
        }

        private void LanczosMethod()
        {
            int m = n;
            List<double> x0 = new List<double>();
            List<double> r0 = new List<double>();
            List<double> alpha = new List<double>();
            List<double> beta = new List<double>();
            List<List<double>> v = new List<List<double>>();
            List<List<double>> w = new List<List<double>>();
            for (int i = 0; i < m + 2; i++)
            {
                v.Add(new List<double>());
                w.Add(new List<double>());
            }
            for (int i = 0; i < m; i++)
                x0.Add(1);
            r0 = Sub(b, MatrToVec(matrix, x0));
            fill(v[0], 0);
            beta.Add(Norma(r0));
            beta.Add(0);
            alpha.Add(0);
            v[1] = Division(r0, beta[0]);
            for (int j = 1; j <=m; j++)
            {
                w[j] = Sub(MatrToVec(matrix, v[j]), Mult(v[j - 1], beta[j]));
                alpha.Add(ScalarMult(w[j], v[j]));
                beta.Add(Norma(w[j]));
                if (beta[j + 1] == 0)
                {
                    m = j;
                    break;
                }
                v[j + 1] = Division(w[j], beta[j + 1]);
            }

            List<double> y = new List<double>();
            for (int i = 0; i < m; i++)
                y.Add(0);
            List<double> f = new List<double>();
            f.Add(beta[0]);
            for (int i = 1; i <m; i++)
                f.Add(0);

            List<double> a1 = new List<double>();
            List<double> b1 = new List<double>();
            List<double> c1 = new List<double>();

            for (int i = 1; i <= m; i++)  //создание трехдиагональной матрицы
            {
                b1.Add(beta[i + 1]);
                a1.Add(beta[i]);
                c1.Add(alpha[i]);
            }
            b1[m - 1] = 0;

            double tmp;                // метод прогонки
            for (int i = 1; i < m; i++)
            {
                tmp = a1[i] / c1[i - 1];
                c1[i] = c1[i] - tmp * b1[i - 1];
                f[i] = f[i] - tmp * f[i - 1];
            }

            y[n - 1] = f[n - 1] / c1[n - 1];

            for (int i = m - 2; i >= 0; i--)
                y[i] = (f[i] - b1[i] * y[i + 1]) / c1[i];


            List<double> x = new List<double>();
            for (int i = 0; i < m; i++)
                x.Add(0);

            for (int i = 1; i <= m; i++)        //результат
                x = Sum(x, Mult(v[i],y[i-1]));
            x = Sum(x, x0);

        }
    }
}

Матрица:

5 - 1  -7 |  9
-1  5  9  | -6
-7  9  -8 |  2

Ожидаемый результат:

 0,8533
 0,2527
-0,7122

Результат:

 0.9670
-0.5740
-2.1193