加速递归行列式算法

Question

如何加速这个递归函数？当它达到10x10矩阵时，解决问题需要一分钟左右的时间。我还包括事件函数，以便您可以看到何时进行计算。

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

保持相同的算法可以做得更好一些，但至少是O(n!)（可能更糟糕），所以无论您对其进行多少优化，高阶矩阵都会变慢。注意我在MSVC 2010中做了基准时间，并且仅用于粗略比较。每次更改都是累积的，随着您进入列表并与原始算法进行比较。

1. 跳过山口检查 - 作为苏尔特建议，删除该得到我们的1％的速度增长。
2. 添加3x3的案例 - 添加另一个直接检查一个3x3矩阵得到我们最，55％
3. 变化POW（） - 在pow()电话更改为(r % 2 ? -1.0 : 1.0)会让我们多一点点，57 ％
4. 更改来切换 - 订单支票更改为开关会让我们多一点点，58％
5. 添加4x4的案例 - 添加另一个直接检查4x4矩阵得到更多，85％

的事情，不工作包括：

1. 的memcpy - 作为苏尔特建议这实际上失去一个很好的协议的速度，-100％
2. 线程 - 创建order线程根本无法正常工作，-160％

我希望使用线程可以让我们显着提高性能，但即使使用所有优化它比原来慢。我认为所有内存的复制使它不是很平行。

添加了3x3和4x4的情况下效果最好，并且是速度超过x6的主要原因。从理论上讲，你可以添加更多明确的案例（可能通过创建一个程序来输出所需的代码）来进一步降低速度。当然，在某种程度上，这种方式打破了使用递归算法的目的。

为了获得更多性能，您可能需要考虑另一种算法。从理论上讲，你可以通过管理自己的堆栈来将递归函数改变为迭代函数，但这是一项相当大的工作，无法保证性能提高。

Answer 2

这可能是一个问题branch mispredict（see also）。测试

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

据我所知，并不需要。

该测试在每个循环中失败1/order次，并且对于小的order占优势，这就是为什么较大的N不会受到如此影响的原因。时机仍然很大O（N！^ 3）（afaik）所以不要期待奇迹。

 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++; 
     } 

该算法在进入L2缓存时会进一步减速，最迟在N = 64时。

此外，矩阵副本可能无效，对于大型order而言，这可能更有效，代价是低效率低于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++; 
     } 

对原始代码进行测试，以获得我确定索引正确的行列式！