C语言实现常见的矩阵运算函数

1.矩阵转置函数

[cpp] view plain copy

print?

1. void matrix_t(double **a_matrix, const double **b_matrix, int krow, int kline)  
2. ////////////////////////////////////////////////////////////////////////////  
3. //  a_matrix:转置后的矩阵  
4. //  b_matrix:转置前的矩阵  
5. //  krow    :行数  
6. //  kline   :列数  
7. ////////////////////////////////////////////////////////////////////////////  
8. {  
9.     int k, k2;     
10.   
11.     for (k = 0; k < krow; k++)  
12.     {  
13.         for(k2 = 0; k2 < kline; k2++)  
14.         {  
15.             a_matrix[k2][k] = b_matrix[k][k2];  
16.         }  
17.     }  
18. }  

2.矩阵加（减）法函数

[cpp] view plain copy

print?

1. void matrix_a(double **a_matrix, const double **b_matrix, const double **c_matrix,   
2.                     int krow, int kline, int ktrl)  
3. ////////////////////////////////////////////////////////////////////////////  
4. //  a_matrix=b_matrix+c_matrix  
5. //   krow   :行数  
6. //   kline  :列数  
7. //   ktrl   :大于0: 加法  不大于0:减法  
8. ////////////////////////////////////////////////////////////////////////////  
9. {  
10.     int k, k2;  
11.   
12.     for (k = 0; k < krow; k++)  
13.     {  
14.         for(k2 = 0; k2 < kline; k2++)  
15.         {  
16.             a_matrix[k][k2] = b_matrix[k][k2]  
17.                 + ((ktrl > 0) ? c_matrix[k][k2] : -c_matrix[k][k2]);   
18.         }  
19.     }  
20. }  

3.矩阵乘法函数

[cpp] view plain copy

print?

1. void matrix_m(double **a_matrix, const double **b_matrix, const double **c_matrix,  
2.                 int krow, int kline, int kmiddle, int ktrl)  
3. ////////////////////////////////////////////////////////////////////////////  
4. //  a_matrix=b_matrix*c_matrix  
5. //  krow  :行数  
6. //  kline :列数  
7. //  ktrl  : 大于0:两个正数矩阵相乘 不大于0:正数矩阵乘以负数矩阵  
8. ////////////////////////////////////////////////////////////////////////////  
9. {  
10.     int k, k2, k4;  
11.     double stmp;  
12.   
13.     for (k = 0; k < krow; k++)       
14.     {  
15.         for (k2 = 0; k2 < kline; k2++)     
16.         {  
17.             stmp = 0.0;  
18.             for (k4 = 0; k4 < kmiddle; k4++)    
19.             {  
20.                 stmp += b_matrix[k][k4] * c_matrix[k4][k2];  
21.             }  
22.             a_matrix[k][k2] = stmp;  
23.         }  
24.     }  
25.     if (ktrl <= 0)     
26.     {  
27.         for (k = 0; k < krow; k++)  
28.         {  
29.             for (k2 = 0; k2 < kline; k2++)  
30.             {  
31.                 a_matrix[k][k2] = -a_matrix[k][k2];  
32.             }  
33.         }  
34.     }  
35. }  

4.矩阵求逆函数

[cpp] view plain copy

print?

1. int  matrix_inv(double **a_matrix, int ndimen)  
2. ////////////////////////////////////////////////////////////////////////////  
3. //  a_matrix:矩阵  
4. //  ndimen :维数  
5. ////////////////////////////////////////////////////////////////////////////  
6. {  
7.     double tmp, tmp2, b_tmp[20], c_tmp[20];  
8.     int k, k1, k2, k3, j, i, j2, i2, kme[20], kmf[20];  
9.     i2 = j2 = 0;  
10.   
11.     for (k = 0; k < ndimen; k++)    
12.     {  
13.         tmp2 = 0.0;  
14.         for (i = k; i < ndimen; i++)    
15.         {  
16.             for (j = k; j < ndimen; j++)    
17.             {  
18.                 if (fabs(a_matrix[i][j] ) <= fabs(tmp2))   
19.                     continue;  
20.                 tmp2 = a_matrix[i][j];  
21.                 i2 = i;  
22.                 j2 = j;  
23.             }    
24.         }  
25.         if (i2 != k)   
26.         {  
27.             for (j = 0; j < ndimen; j++)     
28.             {  
29.                 tmp = a_matrix[i2][j];  
30.                 a_matrix[i2][j] = a_matrix[k][j];  
31.                 a_matrix[k][j] = tmp;  
32.             }  
33.         }  
34.         if (j2 != k)   
35.         {  
36.             for (i = 0; i < ndimen; i++)    
37.             {  
38.                 tmp = a_matrix[i][j2];  
39.                 a_matrix[i][j2] = a_matrix[i][k];  
40.                 a_matrix[i][k] = tmp;  
41.             }      
42.         }  
43.         kme[k] = i2;  
44.         kmf[k] = j2;  
45.         for (j = 0; j < ndimen; j++)    
46.         {  
47.             if (j == k)     
48.             {  
49.                 b_tmp[j] = 1.0 / tmp2;  
50.                 c_tmp[j] = 1.0;  
51.             }  
52.             else   
53.             {  
54.                 b_tmp[j] = -a_matrix[k][j] / tmp2;  
55.                 c_tmp[j] = a_matrix[j][k];  
56.             }  
57.             a_matrix[k][j] = 0.0;  
58.             a_matrix[j][k] = 0.0;  
59.         }  
60.         for (i = 0; i < ndimen; i++)    
61.         {  
62.             for (j = 0; j < ndimen; j++)    
63.             {  
64.                 a_matrix[i][j] = a_matrix[i][j] + c_tmp[i] * b_tmp[j];  
65.             }    
66.         }  
67.     }  
68.     for (k3 = 0; k3 < ndimen;  k3++)     
69.     {  
70.         k  = ndimen - k3 - 1;  
71.         k1 = kme[k];  
72.         k2 = kmf[k];  
73.         if (k1 != k)     
74.         {  
75.             for (i = 0; i < ndimen; i++)    
76.             {  
77.                 tmp = a_matrix[i][k1];  
78.                 a_matrix[i][k1] = a_matrix[i][k];  
79.                 a_matrix[i][k] = tmp;  
80.             }    
81.         }  
82.         if (k2 != k)     
83.         {  
84.             for(j = 0; j < ndimen; j++)    
85.             {  
86.                 tmp = a_matrix[k2][j];  
87.                 a_matrix[k2][j] = a_matrix[k][j];  
88.                 a_matrix[k][j] = tmp;  
89.             }  
90.         }  
91.     }  
92.     return (0);  
93. }  

5.矩阵乔里斯基分解函数

[cpp] view plain copy

print?

1. void chol(double **a_matrix, const double **b_matrix, int ndimen)  
2. ////////////////////////////////////////////////////////////////////////////  
3. //  输入参数:  
4. //      b_matrix:  对称正定方阵    ndimen: 矩阵维数  
5. //  返回值:  
6. //      a_matrix: 下三角矩阵  
7. ////////////////////////////////////////////////////////////////////////////  
8. {  
9.     int i, j, r;  
10.     double m = 0;     
11.     static double **c_matrix;  
12.     static int flag = 0;  
13.   
14.     if (flag == 0)  
15.     {  
16.         flag = 1;  
17.         c_matrix = (double **)malloc(ndimen * sizeof(double *));  
18.   
19.         for (i = 0; i < ndimen; i++)  
20.             c_matrix[i] = (double *)malloc(ndimen * sizeof(double));  
21.     }  
22.   
23.     for (i = 0; i < ndimen; i++)  
24.     {  
25.         for (j = 0; j < ndimen; j++)   
26.             c_matrix[i][j] = 0;  
27.     }  
28.   
29.     c_matrix[0][0] = sqrt(b_matrix[0][0]);  
30.   
31.     for (i = 1; i < ndimen; i++)  
32.     {  
33.         if (c_matrix[0][0] != 0)   
34.             c_matrix[i][0] = b_matrix[i][0] / c_matrix[0][0];  
35.     }  
36.   
37.     for (i = 1; i < ndimen; i++)  
38.     {  
39.         for (r = 0; r < i; r++)      m = m + c_matrix[i][r] * c_matrix[i][r];  
40.   
41.         c_matrix[i][i] = sqrt(b_matrix[i][i] - m);  
42.         m = 0.0;  
43.   
44.         for (j = i + 1; j < ndimen; j++)  
45.         {  
46.             for (r = 0; r < i; r++)      m = m + c_matrix[i][r] * c_matrix[j][r];  
47.             c_matrix[j][i] = (b_matrix[i][j] - m) / c_matrix[i][i];  
48.             m = 0;  
49.         }  
50.     }  
51.   
52.     for (i = 0; i < ndimen; i++)  
53.     {  
54.         for (j = 0; j < ndimen; j++)   
55.             a_matrix[i][j] = c_matrix[i][j];  
56.     }  
57. }