Метод бисопряженных градиентов для решения СЛАУ

            double eps = 0.0001;
            int Rows = grid1.RowDefinitions.Count;
            int m1columns_m1rows = comboBox.SelectedIndex + 2;
            double[] Xk = new double[m1columns_m1rows];//Корни
            double[] Rk = new double[m1columns_m1rows];//вектор невязки. Отличие полученных свободных членов от требуемых
            double[] Sz = new double[m1columns_m1rows];
            double[] Zk = new double[m1columns_m1rows];//Вектор спуска
            double[] Pk = new double[m1columns_m1rows];
            double[] St = new double[m1columns_m1rows];
            double[,] t = Transpose(matrix_left); 
            double[] Sk = new double[m1columns_m1rows];
            double alpha, beta, mf;
            double Spr, Spr1, Spz;
            int i, j;
            double tmp,a,b,c;
            double[] prev_x = new double[m1columns_m1rows];
 
            /* Вычисляем сумму квадратов элементов вектора right_matrix(правой части)*/
            for (mf = 0, i = 0; i < Rows; i++)
            {
                mf += matrix_right[0, i]* matrix_right[0, i];
            }
 
            /* Задаем начальное приближение корней. В Хk хранятся значения корней
            * к-й итерации. */
            for (i = 0; i < Rows; i++)
            {
                Xk[i] = 0;
            }
 
            /* Задаем начальное значение r0 и z0. */
            for (i = 0; i < Rows; i++)
            {
                for (Sz[i] = 0, j = 0; j < Rows; j++)
                {
                    Sz[i] += matrix_left[i, j] * Xk[j];
                }
                Rk[i] = matrix_right[0, i] - Sz[i];
                Zk[i] = Rk[i];
                Sk[i] = Rk[i];
                Pk[i] = Rk[i];
            }
            do
            {
                /* Вычисляем числитель и знаменатель для коэффициента
                * alpha = (pk-1,rk-1)/(Azk-1,sk-1) */
                Spz = 0;
                Spr = 0;
 
                for (i = 0; i < Rows; i++)
                {
                    for (Sz[i] = 0,St[i] = 0, j = 0; j < Rows; j++)
                    {
                        Sz[i] += matrix_left[i, j] * Zk[j];
                        St[i] += t[i, j] * Sk[j];
                    }
                    Spz += Sz[i] * Sk[i];
                    Spr += Pk[i] * Rk[i];
                }
                alpha = Spr / Spz;             /*  alpha    */ //величина смещения
 
 
                /* Вычисляем вектор решения: xk = xk-1+ alpha * zk-1,
                вектор невязки: rk = rk-1 - alpha * A * zk-1 и числитель для beta равный (rk,rk) */
                Spr1 = 0;
                a = 0;
                b = 0;
                for (i = 0; i < Rows; i++)
                {
                    Xk[i] += alpha * Zk[i];
                    Rk[i] -= alpha * Sz[i];
                    Pk[i] -= alpha * St[i];
                    Spr1 += Pk[i] * Rk[i];
                }
 
                /* Вычисляем  beta  beta = (pk * rk)/(pk-1 * rk-1) (*/
                beta = Spr1 / Spr;
 
                /* Вычисляем вектор спуска: zk = rk+ beta * zk-1 */
                for (i = 0; i < Rows; i++)
                {
                    Zk[i] = Rk[i] + beta * Zk[i];
                    Sk[i] = Pk[i] + beta * Sk[i];
                }
 
            }
            /* Проверяем условие выхода из итерационного цикла  */
            while (Spr1/mf > eps);
            int R;
            for (i = 0; i < Rows; i++)
            {
                R = i + 1;
                textBox3.Text += "x[" + (i + 1) + "]= " + Xk[i].ToString() + "\n";
            }

   public static double[,] Transpose(double[,] m)
        {
            double[,] t = new double[m.GetLength(1), m.GetLength(0)];
            for (int i = 0; i < m.GetLength(0); i++)
                for (int j = 0; j < m.GetLength(1); j++)
                    t[j, i] = m[i, j];
 
            return t;
        }