常用矩陣計算C語言程式碼
阿新 • • 發佈：2019-01-06
參考資料：
行列式：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;
}