C# Strassen's Matrix Multiplication Algorithm. Works slowly

Question

C++ / C#

Please help to optimize the algorithm.

I implemented this algorithm, but for some reason it works for me much slower than the traditional n ^ 3 (tested on 1000x1000 matrices). The algorithm is implemented for square matrices, since no more is required.
For the algorithm, a Matrix class was created that contains an int array and three pointers characterizing which part of the array is being worked on (I , J - upper left corner, n - how many elements from this corner to take for work). This is done in order to be able to take submatrices without any special problems and push the array by reference.
If arrays with a size less than or equal to 32 are transferred to the recursion, then they are multiplied according to the Winograd algorithm.

Description of the algorithm itself on Wiki.

The algorithm multiplies correctly, tested on a variety of tests.

Code under the cut.

public class Matrix
{
    public int I { get; set; }
    public int J { get; set; }
    public int[,] mas { get; set; }
    public int n{get;set;}

    public Matrix(int Size)
    {
       mas = new int[Size, Size];
       I = 0; J = 0; n = Size;
    }
    
    public Matrix(int[,] mas,int I,int J,int n)
    {
        this.I = I;
        this.J = J;
        this.n = n;
        this.mas = mas;
    }

    public Matrix(int[,] mas)
    {
        this.I = 0;
        this.J = 0;
        this.mas = mas;
        this.n = mas.GetLength(0);
      
    }

    public Matrix GetLeftTopSubMatrix()
    {
        return new Matrix(mas,I,J,n/2);
    }

    public Matrix GetLeftDownSubMatrix()
    {
        return new Matrix(mas, I+n/2, J, n / 2);
    }

    public Matrix GetRightTopSubMatrix()
    {
        return new Matrix(mas, I, J+n/2, n / 2);
    }

    public Matrix GetRightDownSubMatrix()
    {
        return new Matrix(mas, I + n / 2, J + n / 2, n / 2);
    }

 
}

        public int[,] Shtrassen(Matrix mat1, Matrix mat2)
        {

            if (mat1.n > 2)
            {
                int n2 = mat1.n / 2;


                Matrix P1 = new Matrix(Shtrassen(PlusMatrix(mat1.GetLeftTopSubMatrix(), 
mat1.GetRightDownSubMatrix()), PlusMatrix(mat2.GetLeftTopSubMatrix(), mat2.GetRightDownSubMatrix())));
                Matrix P2 = new Matrix(Shtrassen(PlusMatrix(mat1.GetLeftDownSubMatrix(), 
mat1.GetRightDownSubMatrix()), mat2.GetLeftTopSubMatrix()));
                Matrix P3 = new Matrix(Shtrassen(mat1.GetLeftTopSubMatrix(), 
MinusMatrix(mat2.GetRightTopSubMatrix(), mat2.GetRightDownSubMatrix())));
                Matrix P4 = new Matrix(Shtrassen(mat1.GetRightDownSubMatrix(), 
MinusMatrix(mat2.GetLeftDownSubMatrix(), mat2.GetLeftTopSubMatrix())));
                Matrix P5 = new Matrix(Shtrassen(PlusMatrix(mat1.GetLeftTopSubMatrix(), 
mat1.GetRightTopSubMatrix()), mat2.GetRightDownSubMatrix()));
                Matrix P6 = new Matrix(Shtrassen(MinusMatrix(mat1.GetLeftDownSubMatrix(), 
mat1.GetLeftTopSubMatrix()), PlusMatrix(mat2.GetLeftTopSubMatrix(), mat2.GetRightTopSubMatrix())));
                Matrix P7 = new Matrix(Shtrassen(MinusMatrix(mat1.GetRightTopSubMatrix(), 
mat1.GetRightDownSubMatrix()), PlusMatrix(mat2.GetLeftDownSubMatrix(), mat2.GetRightDownSubMatrix())));



                Matrix R = PlusMatrix(PlusMatrix(P1, P4), MinusMatrix(P7, P5));
                Matrix S = PlusMatrix(P3, P5);
                Matrix T = PlusMatrix(P2, P4);
                Matrix U = PlusMatrix(PlusMatrix(P3, P6), MinusMatrix(P1, P2));


                Matrix Result = new Matrix(mat1.n);


                for (int i = 0; i < n2; i++) for (int j = 0; j < n2; j++)
                    {
                        Result.mas[i, j] = R.mas[i, j];
                        Result.mas[i, j + n2] = S.mas[i, j];
                        Result.mas[i + n2, j] = T.mas[i, j];
                        Result.mas[i + n2, j + n2] = U.mas[i, j];
                    }


                return Result.mas;
            }
            else
            {
                int[,] Result = Vinograd(mat1, mat2);
                return Result;
            }

            int[,] Result1 = new int[1, 1];

            return Result1;
        } 

        public Matrix PlusMatrix(Matrix mat1, Matrix mat2)
        {
            Matrix mat3 = new Matrix(mat1.n);

            for (int i = 0; i < mat1.n; i++) for (int j = 0; j < mat1.n; j++) mat3.mas[i, j] = 
mat1.mas[i + mat1.I, j + mat1.J] + mat2.mas[i + mat2.I, j + mat2.J];
            return mat3;
        }

Answer 1

didn’t look at the algorithm, but already purely by code:
1. int[,] mas - if you need speed, then use [][] or even better just a flat array
2. public int n {get; set;} - auto property - evil

Answer 2

Run a profiler and see where it works slowly.
Well, or stick times and output them to a file, and there already analyze what slows down.

Answer 3

In general, it can be seen that at each step of the recursion, many objects of the Matrix class are created. And memory allocation also takes time. Perhaps this is the problem.
How much slower is this implementation compared to what you would expect from it?

Answer 4

It works slower somewhere in 3.5 - 4 times. The constant of the algorithm is quite large, as are the memory requirements, but on a matrix size of 1000x1000, it should already exactly overtake the traditional multiplication algorithm.
There are suspicions that in the created objects of the Matrix class, the array is not pushed by reference. If this is the case, then in this case you will have to remove this class and work directly with the array, although in this case the code will look unreadable. Now I’ll try to do it ... I
drove it into the profiler, if I can correctly interpret its report, then 70% of the processor time is spent on multiplication by the Winograd algorithm in the case if the matrix size is less than 32

Answer 5

I can suggest the implementation of matrix multiplication on APL (I will do a dll, + I will have to throw a couple of dlls into the bin folder), if the matrices are large, then the speed increase is many times.
For example: 2000 by 2000 on APL - about 8 seconds, on C # there is almost a minute or even more. (I don’t honestly remember the exact results, but my friend and I set up such experiments) =)

Answer 6

Issue resolved. I did not immediately think to check the operation of the algorithm on matrices with a size greater than 1024x1024.
As a result, it turned out to be faster at 4096x4096. I optimized the code according to the advice, removed the Matrix wrapper class, as a result, the algorithm became faster than the traditional one already at 1024x1024, although it is written in the literature that it is preferable to use the Strassen algorithm even at sizes larger than 128x128. But faster, without using caching, I’m afraid I won’t succeed in C #.
Thank you for the suggestion to upload the dll) This is a training project, but unfortunately none of the teachers could help me with it.
Thank you very much for the help)