Возведение большой матрицы, читаемой с файла, в степень N, и запись обратно в файл на C#

using System;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.CompilerServices;
using System.Threading;
 
namespace ConsoleApplication18
{
    public class Matrix : IEnumerable<double>
    {
        const string DimentionExceptionMessage = "Invalid matrix dimensions";
        public static Matrix GetE(int n)
        {
            var result = new Matrix(n, n);
            for (int i = 0; i < n; i++)
            {
                result[i, i] = 1.0;
            }
            return result;
        }
 
        private readonly double[,] _value;
 
        public Matrix(int rows, int cols)
            : this(new double[rows, cols])
        {
        }
 
        public Matrix(double[,] value)
        {
            _value = (double[,])value.Clone();
        }
 
        public static Matrix operator *(Matrix one, Matrix other)
        {
            return new Matrix(Multiply(one._value, other._value));
        }
 
        public static Matrix operator *(double scalar, Matrix matrix)
        {
            var result = (double[,])matrix._value.Clone();
            for (int i = 0; i < result.GetLength(0); i++)
            {
                for (int j = 0; j < result.GetLength(1); j++)
                {
                    result[i, j] *= scalar;
                }
            }
            return result;
        }
 
        public static Matrix operator *(Matrix matrix, double scalar)
        {
            return scalar * matrix;
        }
 
        public static Matrix operator +(Matrix one, Matrix another)
        {
            return MatrixSum(one, another, 1);
        }
 
        public static Matrix operator -(Matrix one, Matrix another)
        {
            return MatrixSum(one, another, -1);
        }
 
        private static Matrix MatrixSum(Matrix one, Matrix another, int sign)
        {
            var m = GetDim(one, another, 0);
            int n = GetDim(one, another, 1);
            var res = new double[m, n];
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    res[i, j] = one[i, j] + another[i, j] * sign;
                }
            }
            return new Matrix(res);
        }
 
        private static int GetDim(Matrix one, Matrix another, int dim)
        {
            int dimOne = one._value.GetLength(dim);
            int dimAnother = another._value.GetLength(dim);
 
            if (dimOne != dimAnother)
                throw new ArgumentException(DimentionExceptionMessage);
            return dimOne;
        }
 
        private int GetDimIfSquare()
        {
            int x = _value.GetLength(0);
            int y = _value.GetLength(1);
            if (x != y)
                throw new ArgumentException(DimentionExceptionMessage);
            return x;
        }
 
        public static Matrix operator /(Matrix matrix, double scalar)
        {
            return matrix * (1.0 / scalar);
        }
 
        public double this[int i, int j]
        {
            get { return _value[i, j]; }
            set { _value[i, j] = value; }
        }
 
        public IEnumerator<double> GetEnumerator()
        {
            // ReSharper disable once LoopCanBeConvertedToQuery
            foreach (double d in _value)
            {
                yield return d;
            }
        }
 
        IEnumerator IEnumerable.GetEnumerator()
        {
            return GetEnumerator();
        }
 
        public static implicit operator Matrix(double[,] doubles)
        {
            return new Matrix(doubles);
        }
 
        public static explicit operator double[,](Matrix matrix)
        {
            return matrix._value;
        }
        public int GetLength(int dim)
        {
            return _value.GetLength(dim);
        }
 
        private static double[,] Multiply(double[,] a, double[,] b)
        {
            if (a.GetLength(1) != b.GetLength(0))
                throw new ArgumentException(DimentionExceptionMessage);
            return SafeMultiply(a, b, a.GetLength(0), b.GetLength(1), a.GetLength(1));
        }
 
 
        public Matrix Pow(int power)
        {
            int dim = GetDimIfSquare();
            var result = GetE(dim)._value;
            var tmp = _value;
            while (power != 0)
            {
                if (power % 2 != 0)
                {
                    result = SafeMultiply(tmp, result, dim, dim, dim);
                    power -= 1;
                }
                tmp = SafeMultiply(tmp, tmp, dim, dim, dim);
                power /= 2;
            }
            return result;
        }
 
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static double[,] SafeMultiply(double[,] a, double[,] b, int imax, int jmax, int kmax)
        {
            double[,] c = new double[imax, jmax];
            int blockCount = Environment.ProcessorCount;
            int blockSize = (int) Math.Ceiling((double) jmax/blockCount);
            int completedBlocks = 0;
            var reset = new ManualResetEvent(false);
            for (int block = 0; block < blockCount; block++)
            {
                int from = block*blockSize;
                int to = Math.Min(from + blockSize, jmax);
                ThreadPool.QueueUserWorkItem(state =>
                {
                    var bcolj = new double[kmax];
                    for (int j = from; j < to; j++)
                    {
                        for (var k = 0; k < kmax; k++)
                            bcolj[k] = b[k, j];
 
                        for (var i = 0; i < imax; i++)
                        {
                            double s = 0;
                            for (var k = 0; k < kmax; k++)
                                s += a[i, k] * bcolj[k];
                            c[i, j] = s;
                        }
                    }
                    if (Interlocked.Increment(ref completedBlocks) == blockCount)
                    {
                        reset.Set();
                    }
                });
            }
            reset.WaitOne();
            return c;
        }
    }
}