常用矩陣計算C語言代碼

　　參考資料：

　　行列式：http://zh.wikipedia.org/wiki/行列式#.E4.BB.A3.E6.95.B0.E4.BD.99.E5.AD.90.E5.BC.8F

　　伴随矩陣：http://zh.wikipedia.org/wiki/伴随矩陣

　　餘因子矩陣：http://zh.wikipedia.org/wiki/餘因子矩陣

　　逆矩陣：http://zh.wikipedia.org/wiki/逆矩陣

　　

　　 關于求一個矩陣的行列式，網上好多代碼其實都是有問題的，我看到好多求行列式的代碼都隻是簡單地把所有正對角線的元素乘起來再求和，然後在減去所有負對角矩陣的相應元素的乘積。這種方法在矩陣的階大于等于4的時候是有問題的，漏掉了好多因子。正确的做法是參照行列式的定義，可以檢視文章頂部的參考資料。

　　

　　 求矩陣行列式的正确代碼如下：

　　//--------------------------------------------------------
　　//功能：求矩陣 n X n 的行列式
　　//入口參數：矩陣首位址 p；矩陣行數 n
　　//傳回值：矩陣的行列式值
　　//--------------------------------------------------------
　　double determinant(double *p, int n) 
　　{
　　 int *list = new int[n];
　　 for (int i = 0; i < n; i++)
　　 list[i] = i;
　　 double ret = det(p, n, 0, list, 0.0);
　　 delete[] list;
　　 return ret;
　　}
　　
　　
　　double det(double *p, int n, int k, int list[], double sum)
　　{
　　 if(k >= n) 
　　 { 
　　 int order = inver_order(list, n);
　　 double item = (double)sgn(order);
　　 for (int i = 0; i < n; i++)
　　 {
　　 //item *= p[i][list[i]];
　　 item *= *(p + i * n + list[i]);
　　 }
　　 return sum + item;
　　 } 
　　 else 
　　 { 
　　 for(int i = k; i < n; i++) 
　　 { 
　　 swap(&list[k], &list[i]); 
　　 sum = det(p, n, k+1, list, sum); 
　　 swap(&list[k], &list[i]); 
　　 } 
　　 }
　　 return sum;
　　}
　　
　　void swap(int *a, int *b) 
　　{ 
　　 int m; 
　　 m = *a; 
　　 *a = *b; 
　　 *b = m; 
　　}
　　
　　//求逆序對的個數
　　int inver_order(int list[], int n)
　　{
　　 int ret = 0;
　　 for(int i = 1; i < n; i++) 
　　 for (int j = 0; j < i; j++)
　　 if (list[j] > list[i])
　　 ret++;
　　 return ret;
　　}
　　
　　int sgn(int order)
　　{
　　 return order % 2 ? -1 : 1;
　　}
           
   

　　當然還可以用LU分解法來求，在矩陣的階比較大時，用高斯消元法或者LU分解法求解具有一定的優勢。

　　

　　由于行列式是求矩陣的代數餘子式的基礎，代數餘子式又是求矩陣的伴随矩陣的基礎，求出伴随矩陣之後才可以求矩陣的逆矩陣。A矩陣的逆矩陣等于A矩陣的伴随矩陣除以A矩陣的行列式。

　　

　　求矩陣某個元素的代數餘子式的代碼如下：

　　//----------------------------------------------------
　　//功能：求k×k矩陣中元素A(mn)的代數餘子式
　　//入口參數：k×k矩陣首位址；元素A的下标m,n; 矩陣行數 k
　　//傳回值： k×k矩陣中元素A(mn)的代數餘子式
　　//----------------------------------------------------
　　double algebraic_cofactor(double *p, int m, int n, int k)
　　{
　　 int len = (k - 1) * (k - 1);
　　 double *cofactor = new double[len];
　　
　　 int count = 0;
　　 int raw_len = k * k;
　　 for (int i = 0; i < raw_len; i++)
　　 if (i / k != m && i % k != n)
　　 *(cofactor + count++) = *(p + i);
　　
　　 double ret = determinant(cofactor, k - 1);
　　 if ((m + n) % 2)
　　 ret = -ret;
　　 delete[] cofactor;
　　 return ret;
　　}
           

　　

　　求伴随矩陣的代碼如下：

　　//----------------------------------------------------
　　//功能：求k×k矩陣的伴随矩陣
　　//入口參數：m是k×k矩陣首位址；矩陣行數 k；輸出參數adj是伴随矩陣的入口位址
　　//傳回值： 無
　　//----------------------------------------------------
　　void adjoint_m(double *m, double *adj, int k)
　　{
　　 int len = k * k;
　　 int count = 0;
　　 for (int i = 0; i < len; i++)
　　 {
　　 *(adj + count++) = algebraic_cofactor(m, i % k, i / k, k);
　　 }
　　}
           

　　

　　求逆矩陣的代碼如下：

　　//----------------------------------------------------
　　//功能：求k×k矩陣的逆矩陣
　　//入口參數：m是k×k矩陣首位址；矩陣行數 k；輸出參數inv是逆矩陣的入口位址
　　//傳回值： 無
　　//----------------------------------------------------
　　void inverse_matrix(double *raw, double *inv, int k)
　　{
　　 double det = determinant(raw, k); //求行列式
　　 if (det == 0)
　　 {
　　 cout << "矩陣不可逆" << endl;
　　 return;
　　 }
　　 adjoint_m(raw, inv, k); //求伴随矩陣
　　 int len = k * k;
　　 for (int i = 0; i < len; i++)
　　 *(inv + i) /= det;
　　}
           

　　

　　兩個矩陣相乘的代碼如下：

　　//--------------------------------------------------------
　　//功能：求矩陣a和b的相乘結果
　　//入口參數：矩陣首位址 a和b；矩陣a行數ra和列數rc；矩陣b的行數rb和列數cb
　　//傳回值：矩陣a和b的相乘結果
　　//--------------------------------------------------------
　　double* m_multiply(double *a, double *b, double *c, int ra, int ca, int rb, int cb)
　　{
　　 if (ca != rb)
　　 {
　　 cout << "矩陣不可乘" << endl;
　　 return NULL;
　　 }
　　
　　 double *ret = c;
　　 if (NULL == ret)
　　 {
　　 ret = new double[ra * cb];
　　 }
　　 for (int i = 0; i < ra; i++)
　　 for (int j = 0; j < cb; j++)
　　 {
　　 //double sum = a[i][0] * b[0][j];
　　 double sum = *(a + i * ca) * (*(b + j));
　　 for (int k = 1; k < ca; k++)
　　 //sum += a[i][k] * b[k][j];
　　 sum += *(a + i*ca + k) * (*(b + k*cb + j));
　　 //c[i][j] = sum;
　　 *(ret + i*cb + j) = sum;
　　 }
　　
　　 return ret;
　　}
           

　　

　　測試程式代碼如下：

　　#include 
　　#include "matrix.h"
　　
　　using namespace std;
　　void main()
　　{
　　 //double a[] = { 2, 6, 3, 1, 0, 2, 5, 8, 4 };
　　 double b[] = { 1, 4, 7, 3, 0, 5, -1, 9, 11 };
　　 double a[] = { -3, 2, -5, -1, 0, -2, 3, -4, 1 };
　　
　　 double det = determinant(a, 3);
　　 cout << det << endl;
　　
　　 det = algebraic_cofactor(a, 1, 2, 3);
　　 cout << det << endl;
　　
　　 double *adj = new double[9];
　　 adjoint_m(a, adj, 3);
　　 for (int i = 0; i < 9; i++)
　　 cout << adj[i] << " ";
　　 cout << endl;
　　
　　 double *c = m_multiply(a, b, NULL, 3, 3, 3, 3);
　　 for (int i = 0; i < 9; i++)
　　 cout << c[i] << " ";
　　 cout << endl;
　　}