Обратная матрица

 public  double[][] MatrixDuplicate(double[][] matrix)
        {
            
            // Предполагается, что матрица не нулевая
            double[][] result = MatrixCreate(matrix.Length, matrix[0].Length);
            for (int i = 0; i < matrix.Length; ++i) // Копирование значений
                for (int j = 0; j < matrix[i].Length; ++j)
                    result[i][j] = matrix[i][j];
            return result;
        }
       public double[][] MatrixDecompose(int rows, int cols)
      {
          double[][] matrix;
           int[] perm;
           int toggle;
  // Разложение LUP Дулитла. Предполагается,
  // что матрица квадратная.
   // для удобства
  double[][] result = MatrixCreate(rows,cols);
  int n = result.Length;
  Print1(result,rows,cols);
  perm = new int[n];
  for (int i = 0; i < n; ++i) { perm[i] = i; }
  toggle = 1;
  for (int j = 0; j < n - 1; ++j) // каждый столбец
  {
    double colMax = Math.Abs(result[j][j]); // Наибольшее значение в столбце j
    int pRow = j;
    for (int i = j + 1; i < n; ++i)
    {
      if (result[i][j] > colMax)
      {
        colMax = result[i][j];
        pRow = i;
      }
    }
    if (pRow != j) // перестановка строк
    {//тут начинается перенумерация столбцов на обратный порядок
      double[] rowPtr = result[pRow];
      result[pRow] = result[j];
      result[j] = rowPtr;
      int tmp = perm[pRow]; // Меняем информацию о перестановке
      perm[pRow] = perm[j];
      perm[j] = tmp;
      toggle = -toggle; // переключатель перестановки строк
    }
    if (Math.Abs(result[j][j]) < 1.0E-20)
      return null;
    for (int i = j + 1; i < n; ++i)
    {
      result[i][j] /= result[j][j];
      for (int k = j + 1; k < n; ++k)
        result[i][k] -= result[i][j] * result[j][k];
    }
  } // основной цикл по столбцу j
  return result;
}

    static double[][] MatrixInverse(double[][] matrix)
    {
      int n = matrix.Length;
      double[][] result = MatrixDuplicate(matrix);
      int[] perm;
      int toggle;
      double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
      if (lum == null)
        throw new Exception("Unable to compute inverse");
      double[] b = new double[n];
      for (int i = 0; i < n; ++i)
      {
        for (int j = 0; j < n; ++j)
        {
          if (i == perm[j])
            b[j] = 1.0;
          else
            b[j] = 0.0;
        }
        double[] x = HelperSolve(lum, b);
        for (int j = 0; j < n; ++j)
          result[j][i] = x[j];
      }
      return result;
    }

    static double[] HelperSolve(double[][] luMatrix,
     double[] b)
    {
      // Решаем luMatrix * x = b
      int n = luMatrix.Length;
      double[] x = new double[n];
      b.CopyTo(x, 0);
      for (int i = 1; i < n; ++i)
      {
        double sum = x[i];
        for (int j = 0; j < i; ++j)
          sum -= luMatrix[i][j] * x[j];
        x[i] = sum;
      }
      x[n - 1] /= luMatrix[n - 1][n - 1];
      for (int i = n - 2; i >= 0; --i)
      {
        double sum = x[i];
        for (int j = i + 1; j < n; ++j)
          sum -= luMatrix[i][j] * x[j];
        x[i] = sum / luMatrix[i][i];
      }
      return x;
    }

        static double[,] One(double[,] hhh )
        { //'единичная матрица
            int h = 2;
            double[,] ggg = new double[h+1, h+1];// еденичная матрица
           
 
 
 
            //   Dim A_edin(k + 1, k + 1) As Double 'затем она перейдёт в обратную
            for (int g = 0; g < h+1 ; g++)
            {
                for (int j = 0; j < h+1 ; j++)
                {
                    if (g == j)
                    { ggg[g, j] = 1; }
                    else
                    { ggg[g, j] = 0; }
                }
            }
 
            //  'единичная матрица
 
            //'обратная матрица
 
 
            double arg=0 ;
            int i1=0;
            //  Dim arg As Double
            //   Dim i1 As Integer
            for (int j = 0; j < h ; )
            {
                for (int i = 0; i < h ; )
                {
                    if (i == j)
                    { j++; i++; }                   
                        arg = hhh[i, j] / hhh[j, j];
                 
                    //   If i = j Then GoTo 50
 
                    for (i1 = 1; i1 < h; )
                    {
                        //  For i1 = 1 To k + 1
                        hhh[i, i1] = hhh[i, i1] - hhh[j, i1] * arg;
                        //          A(i, i1) = A(i, i1) - A(j, i1) * arg
                        ggg[i, i1] = ggg[i, i1] - ggg[j, i1] * arg;
                        i1++;
                    }
                    //        A_edin(i, i1) = A_edin(i, i1) - A_edin(j, i1) * arg
                    //  Next i1
                    i++;
                    j++;
                }
            }
            //    50:         Next i                      Next j
            for (int j = 0; j < h ; )
            {
                for (int i = 0; i < h ; )
                { //  For j = 1 To k + 1
                    //     For i = 1 To k + 1
                    //       Dim arg_2 As Double
                    double arg_2 = 1;
                    if (i == j)
                    {
                        arg_2 = hhh[i, j];
                        for (i1 = 1; i1 < h + 1; )
                        {
                            hhh[i, i1] = hhh[i, i1] / arg_2;
                            ggg[i, i1] = ggg[i, i1] / arg_2;
                            i1++;
 
                        }
                    }
                    i++;
                    j++;
                }
 
 
                /*     If i = j Then
                         arg_2 = A(i, j)
                         For i1 = 1 To k + 1
                             A(i, i1) = A(i, i1) / arg_2 'стала единичной
                             A_edin(i, i1) = A_edin(i, i1) / arg_2 'стала обратной матрице А
                         Next i1
                     End If
 
                 Next i
             Next j
             'обратная матрица
 
     */
 
            }
            return ggg;
        }