一、上三角變換、對角線求值

 1 void
 2 triangle_trans(double **square, int size, int pos)
 3 {
 4   double k1, k2;
 5   k1 = square[pos][pos];
 6   for(int i = pos+1; i < size; i++)
 7     {
 8       k2 = square[i][pos];
 9       for(int j=pos; j < size; j++)
10         {
11           square[i][j] += k2 * square[pos][j] * (-1
 
) / k1;
12         }
13     }
14 }

1 for(int i=0; i < size-1; i++)
2   {
3     triangle_trans(square, size, i);
4   }
5 for(i=0 ; i < size; i++)
6   {
7     ret *= square[i][i]; /* diagonal line */
8   }

　　!TODO: 算法簡要分析

二、行列式按行展開

　　遞歸分解至二階行列式

 1 static std::deque<int
 
> term_stack;
 2 
 3 double
 4 det_A(double **square, int bound_row, int size)
 5 {
 6   int i, j, k;
 7   double result = 0.0;
 8   
 9   if( size - bound_row == 2 )
10     {
11       /*
12        *      |                              |
13        *  D = | a11 [+0][+0]    a12 [+0][+1] |
14        *      | a21 [+1][+0]    a22 [+1][+1] |
 
15        *      |                              |
16        */
17       static const int offset[][2] = { {0, 0}, {0, 1},
18                                        {1, 0}, {1, 1} };
19       int row = bound_row, col_1=-1, col_2=-1;
20       
21       for(j = 0; j < size; j++)
22         {
23           char scope = 1;
24           for(std::deque<int>::iterator it = term_stack.begin(); it!=term_stack.end(); ++it)
25             if( j == (*it) )
26               {
27                 scope=0;
28                 break;
29               }
30           if(scope)
31             {
32               if(col_1 == -1)
33                 col_1 = j;
34               else if(col_2 == -1)
35                 col_2 = j;
36               else assert(0);
37             }
38         }
39 
40       assert(col_1 != -1 && col_2 != -1);
41 
42       return square[row][col_1] * square[row+1][col_2] - square[row][col_2] * square[row+1][col_1];
43     }
44   
45   int col=0;
46 
47   for(j=0; j < size; j++)
48     {
49       char scope = 1;
50       for(std::deque<int>::iterator it = term_stack.begin(); it!=term_stack.end(); ++it)
51         if( j == (*it) )
52           {
53             scope = 0;
54             break;
55           }
56       if(scope)
57         {
58           term_stack.push_front(j);
59           
60           /*
61            * Recursive call to implement row/column expand formula...
62            * D = (-1)^(i+j) * M(i,j)
63            */
64           double acc = square[bound_row][j] * det_A(square, bound_row+1, size);
65           if( (col++) % 2 )
66             acc = -acc;
67 
68           result += acc;
69           term_stack.pop_front();
70         }
71     }
72 
73   return result;
74 }

　　!TODO: 算法簡要分析

三、統計方法

　　通過簡單地調用clock()函數實現對算法運行時間的低精度統計，重復運行N次後求算術平均值。

 1 #define N 100
 2   for(int u=0;u<N;u++)
 3     {
 4       s = clock();
 5 #if USING_EXPAND
 6       term_stack.push_front(-1);
 7       ret = det_A(square, 0, size);
 8 #elif USING_TRIANGLE
 9       ret = 1;
10 
11       for(int i=0; i < size-1; i++)
12         triangle_trans(square, size, i);
13       for(i=0 ; i < size; i++)
14         ret *= square[i][i]; /* diagonal line */
15 #endif
16       d = clock() - s;
17       sum += (double)d * 1000 / CLOCKS_PER_SEC;
18 
19 #if PRINT_TRANS_SPEED
20       std::cout << (double)d / CLOCKS_PER_SEC << std::endl;
21 #endif
22     }
23   std::cout << sum/N << " ms" << std::endl;

四、性能對比

　　如下所示為兩種算法的平均絕對時間對比圖。橫坐標為目標行列式規模（n^2），縱坐標為運行100次所消耗的平均時間，單位為ms。（在Core i7-8550 @3.9GHz測試）

　　INF表示消耗時間超出實驗所允許的設定值（2min）。

　　如下所示為兩種算法的平均絕對時間對比圖。橫坐標為目標行列式規模（n^2），縱坐標為以基於三角變換的算法在147456規模的絕對耗時為基準，各規模絕對耗時與基準的比值。

　　由此可見，基於行（列）展開定理的算法在144（12階）規模下，運行時間就達到了統計設定的上限；而此時基於上三角變換的算法運行時間甚至還低於定時器的精度，開銷遠低於前者，性能差距懸殊。

　　即便在實驗所做的最大規模下，基於上三角變換的算法運行時間不到105ms。

基於上三角變換或行（列）展開定理的n階行列式求值算法及性能評估