Matrix.cs

Go to the documentation of this file.

/*
    Matrix class in C#
    Written by Ivan Kuckir (ivan.kuckir@gmail.com, http://blog.ivank.net)
    Faculty of Mathematics and Physics
    Charles University in Prague
    (C) 2010
    - updated on 14.6.2012 - parsing improved. Thanks to Andy!
    - updated on 3.10.2012 - there was a terrible bug in LU, SoLE and Inversion. Thanks to Danilo Neves Cruz for reporting that!
    
    Updated by Sam McKennoch Early 2014
    - update to correct Resharper Issues
    - updated RandomMatrix to accept double dispersion and probability parameter
    - added GetRow
    - added ConstantMatrix and double + Matrix
    - added ==, !=, Equals, GetHashCode
    - added xml comments throughout
    - Minor refactoring
     
    This code is distributed under MIT licence.
    
    
        Permission is hereby granted, free of charge, to any person
        obtaining a copy of this software and associated documentation
        files (the "Software"), to deal in the Software without
        restriction, including without limitation the rights to use,
        copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the
        Software is furnished to do so, subject to the following
        conditions:

        The above copyright notice and this permission notice shall be
        included in all copies or substantial portions of the Software.

        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
        EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
        OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
        NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
        HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
        WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
        FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
        OTHER DEALINGS IN THE SOFTWARE.
*/

namespace SpikingNeuronNetwork.Lib
{
    using System;
    using System.Collections.Generic;
    using System.Text.RegularExpressions;
    
    public class Matrix
    {
        private const double Epsilon = 1e-5;
        private Matrix L;
        private Matrix U;
        private double _detOfP = 1;
        private int[] _pi;
        private double[,] Mat { get; set; }

        public Matrix(int iRows, int iCols)
        {
            Mat = new double[iRows, iCols];
        }

        public int Rows
        {
            get { return Mat.GetLength(0); }
        }

        public int Cols
        {
            get { return Mat.GetLength(1); }
        }

        public bool IsSquare
        {
            get { return (Rows == Cols); }
        }

        public double this[int iRow, int iCol] // Access this matrix as a 2D array
        {
            get { return Mat[iRow, iCol]; }
            set { Mat[iRow, iCol] = value; }
        }
        
        public override bool Equals(object obj)
        {
            if (ReferenceEquals(null, obj)) return false;
            if (ReferenceEquals(this, obj)) return true;
            return obj.GetType() == GetType() && Equals((Matrix)obj);
        }

        protected bool Equals(Matrix m1)
        {
            return this == m1;
        }

        public override int GetHashCode()
        {
            unchecked
            {
                var hashCode = (Mat != null ? Mat.GetHashCode() : 0);
                hashCode = (hashCode * 397) ^ Rows;
                hashCode = (hashCode * 397) ^ Cols;
                return hashCode;
            }
        }

        public Matrix GetRow(int k)
        {
            var m = new Matrix(1, Cols);
            for (var i = 0; i < Cols; i++) m[0, i] = Mat[k, i];
            return m;
        }

        public Matrix GetCol(int k)
        {
            var m = new Matrix(Rows, 1);
            for (var i = 0; i < Rows; i++) m[i, 0] = Mat[i, k];
            return m;
        }

        public List<int> FindAllIndexOf(double val)
        {
            var indices = new List<int>();
            for (var row = 0; row < Rows; row ++)
            {
                for (var col = 0; col < Cols; col++)
                {
                    if (Math.Abs(Mat[row, col] - val) < Epsilon)
                    {
                        indices.Add(col + row*Cols);
                    }
                }                
            }
            return indices;
        }

        public List<int> FindAllIndexOfNot(double val)
        {
            var indices = new List<int>();
            for (var row = 0; row < Rows; row++)
            {
                for (var col = 0; col < Cols; col++)
                {
                    if (Math.Abs(Mat[row, col] - val) >= Epsilon)
                    {
                        indices.Add(col + row * Cols);
                    }
                }
            }
            return indices;
        }

        public void SetCol(Matrix v, int k)
        {
            for (var i = 0; i < Rows; i++) Mat[i, k] = v[i, 0];
        }

        public void MakeLU()
        {
            if (!IsSquare) throw new MatrixException("The matrix is not square!");
            L = IdentityMatrix(Rows, Cols);
            U = Clone();

            _pi = new int[Rows];
            for (var i = 0; i < Rows; i++) _pi[i] = i;

            var k0 = 0;

            for (var k = 0; k < Cols - 1; k++)
            {
                double p = 0;
                for (var i = k; i < Rows; i++) // find the row with the biggest pivot
                {
                    if (Math.Abs(U[i, k]) > p)
                    {
                        p = Math.Abs(U[i, k]);
                        k0 = i;
                    }
                }
                if (Math.Abs(p - 0) < Epsilon) // samé nuly ve sloupci
                    throw new MatrixException("The matrix is singular!");

                var pom1 = _pi[k];
                _pi[k] = _pi[k0];
                _pi[k0] = pom1; // switch two rows in permutation matrix

                double pom2;
                for (var i = 0; i < k; i++)
                {
                    pom2 = L[k, i];
                    L[k, i] = L[k0, i];
                    L[k0, i] = pom2;
                }

                if (k != k0) _detOfP *= -1;

                for (var i = 0; i < Cols; i++) // Switch rows in U
                {
                    pom2 = U[k, i];
                    U[k, i] = U[k0, i];
                    U[k0, i] = pom2;
                }

                for (var i = k + 1; i < Rows; i++)
                {
                    L[i, k] = U[i, k]/U[k, k];
                    for (var j = k; j < Cols; j++)
                        U[i, j] = U[i, j] - L[i, k]*U[k, j];
                }
            }
        }

        public Matrix SolveWith(Matrix v)
        {
            if (Rows != Cols) throw new MatrixException("The matrix is not square!");
            if (Rows != v.Rows) throw new MatrixException("Wrong number of results in solution vector!");
            if (L == null) MakeLU();

            var b = new Matrix(Rows, 1);
            for (var i = 0; i < Rows; i++) b[i, 0] = v[_pi[i], 0]; // switch two items in "v" due to permutation matrix

            var z = SubsForth(L, b);
            var x = SubsBack(U, z);

            return x;
        }

        public Matrix Invert()
        {
            if (L == null) MakeLU();

            var inv = new Matrix(Rows, Cols);

            for (var i = 0; i < Rows; i++)
            {
                var ei = ZeroMatrix(Rows, 1);
                ei[i, 0] = 1;
                var col = SolveWith(ei);
                inv.SetCol(col, i);
            }
            return inv;
        }

        public double Det()
        {
            if (L == null) MakeLU();
            var det = _detOfP;
            for (var i = 0; i < Rows; i++) det *= U[i, i];
            return det;
        }

        public Matrix GetPermutationMatrix()
        {
            if (L == null) MakeLU();

            var matrix = ZeroMatrix(Rows, Cols);
            for (var i = 0; i < Rows; i++) matrix[_pi[i], i] = 1;
            return matrix;
        }

        public Matrix Clone()
        {
            var matrix = new Matrix(Rows, Cols);
            for (var i = 0; i < Rows; i++)
                for (var j = 0; j < Cols; j++)
                    matrix[i, j] = Mat[i, j];
            return matrix;
        }

        public static Matrix SubsForth(Matrix mA, Matrix b)
        {
            if (mA.L == null) mA.MakeLU();
            var n = mA.Rows;
            var x = new Matrix(n, 1);

            for (var i = 0; i < n; i++)
            {
                x[i, 0] = b[i, 0];
                for (var j = 0; j < i; j++) x[i, 0] -= mA[i, j]*x[j, 0];
                x[i, 0] = x[i, 0]/mA[i, i];
            }
            return x;
        }

        public static Matrix SubsBack(Matrix mA, Matrix b)
        {
            if (mA.L == null) mA.MakeLU();
            var n = mA.Rows;
            var x = new Matrix(n, 1);

            for (var i = n - 1; i > -1; i--)
            {
                x[i, 0] = b[i, 0];
                for (var j = n - 1; j > i; j--) x[i, 0] -= mA[i, j]*x[j, 0];
                x[i, 0] = x[i, 0]/mA[i, i];
            }
            return x;
        }

        public static Matrix ZeroMatrix(int iRows, int iCols)
        {
            var matrix = new Matrix(iRows, iCols);
            for (var i = 0; i < iRows; i++)
                for (var j = 0; j < iCols; j++)
                    matrix[i, j] = 0;
            return matrix;
        }

        public static Matrix ConstantMatrix(int iRows, int iCols, double val)
        {
            var matrix = new Matrix(iRows, iCols);
            for (var i = 0; i < iRows; i++)
                for (var j = 0; j < iCols; j++)
                    matrix[i, j] = val;
            return matrix;
        }

        public static Matrix IdentityMatrix(int iRows, int iCols)
        {
            var matrix = ZeroMatrix(iRows, iCols);
            for (var i = 0; i < Math.Min(iRows, iCols); i++)
                matrix[i, i] = 1;
            return matrix;
        }

        public static Matrix RandomMatrix(int iRows, int iCols, double dispersion, double probability = 1.0)
        {
            var random = new Random();
            var matrix = new Matrix(iRows, iCols);
            for (var i = 0; i < iRows; i++)
                for (var j = 0; j < iCols; j++)
                    matrix[i, j] = (random.NextDouble() <= probability) ? dispersion*2*(random.NextDouble() - 0.5) : 0;
            return matrix;
        }

        public static Matrix Parse(string ps)
        {
            var s = NormalizeMatrixString(ps);
            var rows = Regex.Split(s, "\r\n");
            var nums = rows[0].Split(' ');
            var matrix = new Matrix(rows.Length, nums.Length);
            try
            {
                for (var i = 0; i < rows.Length; i++)
                {
                    nums = rows[i].Split(' ');
                    for (var j = 0; j < nums.Length; j++) matrix[i, j] = double.Parse(nums[j]);
                }
            }
            catch (FormatException)
            {
                throw new MatrixException("Wrong input format!");
            }
            return matrix;
        }

        public override string ToString()
        {
            string s = "";
            for (int i = 0; i < Rows; i++)
            {
                for (int j = 0; j < Cols; j++) s += String.Format("{0,5:0.00}", Mat[i, j]) + " ";
                s += "\r\n";
            }
            return s;
        }

        public static Matrix Transpose(Matrix m)
        {
            var t = new Matrix(m.Cols, m.Rows);
            for (int i = 0; i < m.Rows; i++)
                for (int j = 0; j < m.Cols; j++)
                    t[j, i] = m[i, j];
            return t;
        }

        public static Matrix Power(Matrix m, int pow)
        {
            if (pow == 0) return IdentityMatrix(m.Rows, m.Cols);
            if (pow == 1) return m.Clone();
            if (pow == -1) return m.Invert();

            Matrix x;
            if (pow < 0)
            {
                x = m.Invert();
                pow *= -1;
            }
            else x = m.Clone();

            // ToDo: This would be much faster with a decomposition
            Matrix ret = IdentityMatrix(m.Rows, m.Cols);
            while (pow != 0)
            {
                if ((pow & 1) == 1) ret *= x;
                x *= x;
                pow >>= 1;
            }
            return ret;
        }

        private static void SafeAplusBintoC(Matrix mA, int xa, int ya, Matrix mB, int xb, int yb, Matrix mC, int size)
        {
            for (int i = 0; i < size; i++) // rows
                for (int j = 0; j < size; j++) // cols
                {
                    mC[i, j] = 0;
                    if (xa + j < mA.Cols && ya + i < mA.Rows) mC[i, j] += mA[ya + i, xa + j];
                    if (xb + j < mB.Cols && yb + i < mB.Rows) mC[i, j] += mB[yb + i, xb + j];
                }
        }

        private static void SafeAminusBintoC(Matrix mA, int xa, int ya, Matrix mB, int xb, int yb, Matrix mC, int size)
        {
            for (int i = 0; i < size; i++) // rows
                for (int j = 0; j < size; j++) // cols
                {
                    mC[i, j] = 0;
                    if (xa + j < mA.Cols && ya + i < mA.Rows) mC[i, j] += mA[ya + i, xa + j];
                    if (xb + j < mB.Cols && yb + i < mB.Rows) mC[i, j] -= mB[yb + i, xb + j];
                }
        }

        private static void SafeACopytoC(Matrix mA, int xa, int ya, Matrix mC, int size)
        {
            for (int i = 0; i < size; i++) // rows
                for (int j = 0; j < size; j++) // cols
                {
                    mC[i, j] = 0;
                    if (xa + j < mA.Cols && ya + i < mA.Rows) mC[i, j] += mA[ya + i, xa + j];
                }
        }

        private static void AplusBintoC(Matrix mA, int xa, int ya, Matrix mB, int xb, int yb, Matrix mC, int size)
        {
            for (int i = 0; i < size; i++) // rows
                for (int j = 0; j < size; j++) mC[i, j] = mA[ya + i, xa + j] + mB[yb + i, xb + j];
        }

        private static void AminusBintoC(Matrix mA, int xa, int ya, Matrix mB, int xb, int yb, Matrix mC, int size)
        {
            for (int i = 0; i < size; i++) // rows
                for (int j = 0; j < size; j++) mC[i, j] = mA[ya + i, xa + j] - mB[yb + i, xb + j];
        }

        private static void ACopytoC(Matrix mA, int xa, int ya, Matrix mC, int size)
        {
            for (int i = 0; i < size; i++) // rows
                for (int j = 0; j < size; j++) mC[i, j] = mA[ya + i, xa + j];
        }

        private static Matrix StrassenMultiply(Matrix mA, Matrix mB) // Smart matrix multiplication
        {
            if (mA.Cols != mB.Rows) throw new MatrixException("Wrong dimension of matrix!");

            Matrix mR;

            int msize = Math.Max(Math.Max(mA.Rows, mA.Cols), Math.Max(mB.Rows, mB.Cols));

            if (msize < 32)
            {
                mR = ZeroMatrix(mA.Rows, mB.Cols);
                for (int i = 0; i < mR.Rows; i++)
                    for (int j = 0; j < mR.Cols; j++)
                        for (int k = 0; k < mA.Cols; k++)
                            mR[i, j] += mA[i, k]*mB[k, j];
                return mR;
            }

            int size = 1;
            int n = 0;
            while (msize > size)
            {
                size *= 2;
                n++;
            }

            var h = size/2;


            var mField = new Matrix[n,9];

            /*
             *  8x8, 8x8, 8x8, ...
             *  4x4, 4x4, 4x4, ...
             *  2x2, 2x2, 2x2, ...
             *  . . .
             */

            for (var i = 0; i < n - 4; i++) // rows
            {
                var z = (int) Math.Pow(2, n - i - 1);
                for (int j = 0; j < 9; j++) mField[i, j] = new Matrix(z, z);
            }

            SafeAplusBintoC(mA, 0, 0, mA, h, h, mField[0, 0], h);
            SafeAplusBintoC(mB, 0, 0, mB, h, h, mField[0, 1], h);
            StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 1], 1, mField); // (A11 + A22) * (B11 + B22);

            SafeAplusBintoC(mA, 0, h, mA, h, h, mField[0, 0], h);
            SafeACopytoC(mB, 0, 0, mField[0, 1], h);
            StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 2], 1, mField); // (A21 + A22) * B11;

            SafeACopytoC(mA, 0, 0, mField[0, 0], h);
            SafeAminusBintoC(mB, h, 0, mB, h, h, mField[0, 1], h);
            StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 3], 1, mField); //A11 * (B12 - B22);

            SafeACopytoC(mA, h, h, mField[0, 0], h);
            SafeAminusBintoC(mB, 0, h, mB, 0, 0, mField[0, 1], h);
            StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 4], 1, mField); //A22 * (B21 - B11);

            SafeAplusBintoC(mA, 0, 0, mA, h, 0, mField[0, 0], h);
            SafeACopytoC(mB, h, h, mField[0, 1], h);
            StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 5], 1, mField); //(A11 + A12) * B22;

            SafeAminusBintoC(mA, 0, h, mA, 0, 0, mField[0, 0], h);
            SafeAplusBintoC(mB, 0, 0, mB, h, 0, mField[0, 1], h);
            StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 6], 1, mField); //(A21 - A11) * (B11 + B12);

            SafeAminusBintoC(mA, h, 0, mA, h, h, mField[0, 0], h);
            SafeAplusBintoC(mB, 0, h, mB, h, h, mField[0, 1], h);
            StrassenMultiplyRun(mField[0, 0], mField[0, 1], mField[0, 1 + 7], 1, mField); // (A12 - A22) * (B21 + B22);

            mR = new Matrix(mA.Rows, mB.Cols); // result

            // C11
            for (int i = 0; i < Math.Min(h, mR.Rows); i++) // rows
                for (int j = 0; j < Math.Min(h, mR.Cols); j++) // cols
                    mR[i, j] = mField[0, 1 + 1][i, j] + mField[0, 1 + 4][i, j] - mField[0, 1 + 5][i, j] +
                              mField[0, 1 + 7][i, j];

            // C12
            for (int i = 0; i < Math.Min(h, mR.Rows); i++) // rows
                for (int j = h; j < Math.Min(2*h, mR.Cols); j++) // cols
                    mR[i, j] = mField[0, 1 + 3][i, j - h] + mField[0, 1 + 5][i, j - h];

            // C21
            for (int i = h; i < Math.Min(2*h, mR.Rows); i++) // rows
                for (int j = 0; j < Math.Min(h, mR.Cols); j++) // cols
                    mR[i, j] = mField[0, 1 + 2][i - h, j] + mField[0, 1 + 4][i - h, j];

            // C22
            for (int i = h; i < Math.Min(2*h, mR.Rows); i++) // rows
                for (int j = h; j < Math.Min(2*h, mR.Cols); j++) // cols
                    mR[i, j] = mField[0, 1 + 1][i - h, j - h] - mField[0, 1 + 2][i - h, j - h] +
                              mField[0, 1 + 3][i - h, j - h] + mField[0, 1 + 6][i - h, j - h];

            return mR;
        }

        // function for square matrix 2^N x 2^N

        private static void StrassenMultiplyRun(Matrix mA, Matrix mB, Matrix mC, int l, Matrix[,] f)
            // A * B into C, level of recursion, matrix field
        {
            int size = mA.Rows;
            int h = size/2;

            if (size < 32)
            {
                for (int i = 0; i < mC.Rows; i++)
                    for (int j = 0; j < mC.Cols; j++)
                    {
                        mC[i, j] = 0;
                        for (int k = 0; k < mA.Cols; k++) mC[i, j] += mA[i, k]*mB[k, j];
                    }
                return;
            }

            AplusBintoC(mA, 0, 0, mA, h, h, f[l, 0], h);
            AplusBintoC(mB, 0, 0, mB, h, h, f[l, 1], h);
            StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 1], l + 1, f); // (A11 + A22) * (B11 + B22);

            AplusBintoC(mA, 0, h, mA, h, h, f[l, 0], h);
            ACopytoC(mB, 0, 0, f[l, 1], h);
            StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 2], l + 1, f); // (A21 + A22) * B11;

            ACopytoC(mA, 0, 0, f[l, 0], h);
            AminusBintoC(mB, h, 0, mB, h, h, f[l, 1], h);
            StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 3], l + 1, f); //A11 * (B12 - B22);

            ACopytoC(mA, h, h, f[l, 0], h);
            AminusBintoC(mB, 0, h, mB, 0, 0, f[l, 1], h);
            StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 4], l + 1, f); //A22 * (B21 - B11);

            AplusBintoC(mA, 0, 0, mA, h, 0, f[l, 0], h);
            ACopytoC(mB, h, h, f[l, 1], h);
            StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 5], l + 1, f); //(A11 + A12) * B22;

            AminusBintoC(mA, 0, h, mA, 0, 0, f[l, 0], h);
            AplusBintoC(mB, 0, 0, mB, h, 0, f[l, 1], h);
            StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 6], l + 1, f); //(A21 - A11) * (B11 + B12);

            AminusBintoC(mA, h, 0, mA, h, h, f[l, 0], h);
            AplusBintoC(mB, 0, h, mB, h, h, f[l, 1], h);
            StrassenMultiplyRun(f[l, 0], f[l, 1], f[l, 1 + 7], l + 1, f); // (A12 - A22) * (B21 + B22);

            // C11
            for (int i = 0; i < h; i++) // rows
                for (int j = 0; j < h; j++) // cols
                    mC[i, j] = f[l, 1 + 1][i, j] + f[l, 1 + 4][i, j] - f[l, 1 + 5][i, j] + f[l, 1 + 7][i, j];

            // C12
            for (int i = 0; i < h; i++) // rows
                for (int j = h; j < size; j++) // cols
                    mC[i, j] = f[l, 1 + 3][i, j - h] + f[l, 1 + 5][i, j - h];

            // C21
            for (int i = h; i < size; i++) // rows
                for (int j = 0; j < h; j++) // cols
                    mC[i, j] = f[l, 1 + 2][i - h, j] + f[l, 1 + 4][i - h, j];

            // C22
            for (int i = h; i < size; i++) // rows
                for (int j = h; j < size; j++) // cols
                    mC[i, j] = f[l, 1 + 1][i - h, j - h] - f[l, 1 + 2][i - h, j - h] + f[l, 1 + 3][i - h, j - h] +
                              f[l, 1 + 6][i - h, j - h];
        }

        private static Matrix Multiply(double n, Matrix m) // Multiplication by constant n
        {
            var r = new Matrix(m.Rows, m.Cols);
            for (int i = 0; i < m.Rows; i++)
                for (int j = 0; j < m.Cols; j++)
                    r[i, j] = m[i, j]*n;
            return r;
        }

        private static Matrix Add(Matrix m1, Matrix m2) // Sčítání matic
        {
            if (m1.Rows != m2.Rows || m1.Cols != m2.Cols)
                throw new MatrixException("Matrices must have the same dimensions!");
            var r = new Matrix(m1.Rows, m1.Cols);
            for (int i = 0; i < r.Rows; i++)
                for (int j = 0; j < r.Cols; j++)
                    r[i, j] = m1[i, j] + m2[i, j];
            return r;
        }

        private static string NormalizeMatrixString(string matStr) // From Andy - thank you! :)
        {
            // Remove any multiple spaces
            while (matStr.IndexOf("  ", StringComparison.Ordinal) != -1)
                matStr = matStr.Replace("  ", " ");

            // Remove any spaces before or after newlines
            matStr = matStr.Replace(" \r\n", "\r\n");
            matStr = matStr.Replace("\r\n ", "\r\n");

            // If the data ends in a newline, remove the trailing newline.
            // Make it easier by first replacing \r\n’s with |’s then
            // restore the |’s with \r\n’s
            matStr = matStr.Replace("\r\n", "|");
            while (matStr.LastIndexOf("|", StringComparison.Ordinal) == (matStr.Length - 1))
                matStr = matStr.Substring(0, matStr.Length - 1);

            matStr = matStr.Replace("|", "\r\n");
            return matStr;
        }

        #region operators

        public static Matrix operator -(Matrix m)
        {
            return Multiply(-1, m);
        }

        public static Matrix operator +(double m1, Matrix m2)
        {
            var m3 = new Matrix(m2.Rows, m2.Cols);
            for (var i = 0; i < m2.Rows; i++)
                for (var j = 0; j < m2.Cols; j++)
                    m3[i, j] = m2[i,j] + m1;
            return m2;
        }

        public static Matrix operator +(Matrix m1, double m2)
        {
            return m2 + m1;
        }

        public static Matrix operator +(Matrix m1, Matrix m2)
        {
            return Add(m1, m2);
        }

        public static Matrix operator -(Matrix m1, Matrix m2)
        {
            return Add(m1, -m2);
        }

        public static Matrix operator *(Matrix m1, Matrix m2)
        {
            return StrassenMultiply(m1, m2);
        }

        public static Matrix operator *(double n, Matrix m)
        {
            return Multiply(n, m);
        }

        public static bool operator ==(Matrix m1, Matrix m2)
        {
            if (m1 == null || m2 == null || m1.Rows != m2.Rows || m1.Cols != m2.Cols)
                throw new MatrixException("Matrices must have the same dimensions!");
            for (var i = 0; i < m1.Rows; i++)
                for (var j = 0; j < m1.Cols; j++)
                    if (Math.Abs(m1[i, j] - m2[i, j]) > Epsilon) return false;
            return true;
        }

        public static bool operator !=(Matrix m1, Matrix m2)
        {
            return !(m1 == m2);
        }

        #endregion
    }
}