Speed up recursive determinant algorithm

Question

How do I speed up this recursive function? When it reaches a 10x10 matrix, it takes up a minute or so just to solve a problem. I included the event function as well so you can see when the calculation would take place.

void determinantsFrame::OnCalculateClick(wxCommandEvent &event)
{
    double elem[MAX][MAX]; double det; string test; bool doIt = true;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            test = (numbers[i][j]->GetValue()).mb_str();
            if (test == "")
            {
                doIt = false;
                break;
            }

            for (int k = 0; k < test.length(); k++)
                if (isalpha(test[k]) || test[k] == ' ')
                {
                    doIt = false;
                    break;
                }
                else if (ispunct(test[k]))
                {
                    if (test[k] == '.' && test.length() == 1)
                        doIt = false;
                    else if (test[k] == '.' && test.length() != 1)
                        doIt = true;
                    else if (test[k] != '.')
                        doIt = false;
                }

            if (doIt == false)
                break;
        }
        if (doIt == false)
            break;
    }

    if (doIt)
    {
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                elem[i][j] = static_cast<double>(wxAtof(numbers[i][j]->GetValue()));

        det = determinant(elem, n);
        wxMessageBox(wxString::Format(wxT("The determinant is: %.4lf"),det));
    }
    else
        wxMessageBox(wxT("You may have entered an invalid character. Please try again"));
}

double determinantsFrame::determinant(double matrix[MAX][MAX], int order) // Here's the recursive algorithm
{
    double det = 0; double temp[MAX][MAX]; int row, col;

    if (order == 1)
        return matrix[0][0];
    else if (order == 2)
        return ((matrix[0][0] * matrix[1][1]) - (matrix[0][1] * matrix[1][0]));
    else
    {
        for (int r = 0; r < order; r++)
        {
            col = 0; row = 0;
            for (int i = 1; i < order; i++)
            {
                for (int j = 0; j < order; j++)
                {
                    if (j == r)
                        continue;

                    temp[row][col] = matrix[i][j];
                    col++;

                    if (col == order - 1)
                        col = 0;
                }
                row++;
            }
            det = det + (matrix[0][r] * pow(-1, r) * determinant(temp, order - 1));
        }
        return det;
    }
}

Answer 1

It could be a branch mispredict problem (see also). The test

if (col == order - 1)
    col = 0;

Is not needed as far as I can see.

The test fails 1/order times per loop and dominates for small order, which is why larger N aren't so affected. The timing is still large O(N!^3) (afaik) so don't expect miracles.

        col = 0; row = 0;
        for (int i = 1; i < order; i++) {
            for (int j = 0; j < order; j++) {
                if (j == r)
                    continue;

                temp[row][col] = matrix[i][j];
                col++;

                //if (col == order - 1)
                //    col = 0;
            }
            col = 0; // no need to test
            row++;
        }

The algorithm will get a further slowdown when it hit L2 cache, at latest at N=64.

Also the matrix copy might be ineffective, this could be far more effective for large order at the cost of low effectiveness at low order.

    for (int r = 0; r < order; r++) {
        row = 0;
        for (int i = 1; i < order; i++) {
            memcpy(temp[row], matrix[i], r*sizeof(double)); // if r==0 will this work?
            memcpy(&temp[row][r], &matrix[i][r+1], (order-r-1)*sizeof(double));
            // amount of copied elements r+(order-r-1)=order-1.

            row++;
        }

Make a test with the original code to get the determinant that I got the indexes right!

Answer 2

You can do a bit better with keeping the same algorithm but it is at least O(n!) (probably worse) so higher order matrices will be slow no matter how much you optimize it. Note I did the benchmark times in MSVC 2010 and are there only for rough comparison purposes. Each change is cumulative as you go down the list and is compared to the original algorithm.

1. Skip Col Check -- As Surt suggested, removing this gets us a speed increase of 1%.
2. Add 3x3 Case -- Adding another explicit check for a 3x3 matrix gets us the most, 55%
3. Change pow() -- Changing the pow() call to (r % 2 ? -1.0 : 1.0) gets us a little bit more, 57%
4. Change to switch -- Changing the order check to a switch gets us a little bit more, 58%
5. Add 4x4 Case -- Adding another explicit check for a 4x4 matrix gets more, 85%

Things that don't work include:

1. memcpy -- As Surt suggested this actually looses a good deal of speed, -100%
2. Threads -- Creating order threads doesn't work well at all, -160%

I was hoping that using threads could get us a significant performance increase but even with all the optimization it is slower than the original. I think the copying of all the memory is making it not very parallel.

Added the 3x3 and 4x4 cases has the most effect and are the primary reason for the over x6 increase in speed. In theory you could add more explicit cases (probably by creating a program to output the required code) to reduce the speed even further. Of course, at some point this kind of defeats the purpose of using a recursive algorithm to begin with.

To get more performance you would probably have to consider a different algorithm. In theory you can change the recursive function into an iterative one by managing your own stack but it is considerable work and you aren't guaranteed a performance increase anyways.