求逆矩阵以及两矩阵相乘的算法实现

下面的程序借鉴了多位高贤的代码，玄机逸士加以整理和修改而成。

// matrixcomputation.h

//计算行列式，参数一为存储行列式的数组，参数二为阶数
double calculateDeterminant(double *p,int n);

// 使用Gauss-Jordan消去法求n阶实矩阵的逆矩阵
// 返回结果存放在a中，n是矩阵的阶数
int inverseMatrix(double *a, int n); 

// 矩阵相乘
// 矩阵a[m*n]乘以矩阵b[n*k]，得到矩阵c[m*k]
// 矩阵c[m*k]为返回结果
void multiplyMatrix(double a[], double b[], int m, int n, int k, double c[]);

// matrixcomputation.cpp

#include "matrixcomputation.h"
#include <math.h>

//计算行列式，参数一为存储行列式的数组，参数二为阶数
double calculateDeterminant(double *p,int n) 
{
 if(n<1)
 {
 return -1;
 }
   //若行列式阶数为1，则返回其内的唯一元素值
 if(n==1) 
 {
 return *p;
 }
 //定义临时指针变量指向当前行列式余子式数组
 double *curr = new double[(n-1)*(n-1)]; 
 double result = 0;
 int i, r, c, cc;

   //求当前第一行第i列余子式(按第一行展开)
 for(i = 1; i <= n; ++i) 
 {
 for(r = 2; r <= n; ++r)
 {
 cc = -1;
 for(c = 1; c <= n; ++c)
 {
 if(c != i)
 {
 *(curr + (r-2) * (n-1) + (++cc)) = *(p + (r - 1) * n + (c-1));
 }
 }
 }
      //计算求代数余子式之和，使用递归调用计算其余子式的值
 result +=*(p + i - 1) * pow(-1, 1 + i) * calculateDeterminant(curr, n-1); 
 }
   //释放存储临时变量所占空间
 delete [] curr; 
 return result;
}

// 使用Gauss-Jordan消去法求n阶实矩阵的逆矩阵
// 返回结果存放在a中，n是矩阵的阶数
int inverseMatrix(double *a, int n) 
{ 
int *is, *js, i, j, k, l, u, v; 
double d, p; 
is = new int[n];
js = new int[n];

// 判断矩阵的的行列式是否为0，
// 如果为0，则说明没有逆矩阵
// 否则有逆矩阵
double *q = new double[n * n];
for(i = 0; i < n; i++)
{
for(j = 0; j < n; j++)
{
q[i * n + j] = a[i * n + j];
}
}

if(calculateDeterminant(q, n) == 0)
{
// printf("This matrix does not have inverse matrix...");
return 0;
}
delete [] q;

// 开始计算逆矩阵
for(k = 0; k <= n-1; k++) 
{ 
d = 0.0; 
for(i = k; i <= n-1; i++) 
{
for(j = k; j <= n-1; j++) 
{ 
l = i * n + j; 
p = fabs(a[l]); 
if(p > d) 
{ 
d = p; 
is[k] = i; 
js[k] = j;
} 
}
}

if(is[k] != k)
{
for(j = 0; j <= n-1; j++) 
{ 
u = k * n + j; 
v = is[k] * n + j; 
p = a[u]; 
a[u] = a[v]; 
a[v] = p; 
} 
}
if(js[k] != k) 
{
for(i = 0; i <= n-1; i++) 
{ 
u = i * n + k; 
v = i * n + js[k]; 
p = a[u]; 
a[u] = a[v]; 
a[v] = p; 
}
}
l = k * n + k; 
a[l] = 1.0/a[l]; 
for(j = 0; j <= n-1; j++)
{
if(j != k) 
{ 
u = k * n + j; 
a[u] = a[u] * a[l];
}
}
for(i = 0; i <= n-1; i++)
{
if(i != k)
{
for(j = 0; j <= n-1; j++) 
{
if (j != k) 
{ 
u = i * n + j; 
a[u] = a[u] - a[i * n + k] * a[k * n + j]; 
}
}
}
}
for(i = 0; i <= n-1; i++)
{
if(i != k) 
{ 
u = i * n + k; 
a[u] = -a[u] * a[l];
} 
}
} 
for(k = n-1; k >= 0; k--) 
{ 
if(js[k] != k) 
{
for(j = 0; j <= n-1; j++) 
{ 
u = k * n + j; 
v = js[k] * n + j; 
p = a[u]; 
a[u] = a[v]; 
a[v] = p; 
}
}
if(is[k] != k) 
for(i = 0; i <= n-1; i++) 
{ 
u = i * n + k; 
v = i * n + is[k]; 
p = a[u]; 
a[u] = a[v]; 
a[v] = p; 
} 
} 
delete [] is;
delete [] js;
return 1; 
} 

// 矩阵相乘
// 矩阵a[m*n]乘以矩阵b[n*k]，得到矩阵c[m*k]
// 矩阵c[m*k]为返回结果
void multiplyMatrix(double a[], double b[], int m, int n, int k, double c[]) 
{ 
int i,j,l,u; 
for(i = 0; i <= m-1; i++) 
{
for(j = 0; j <= k-1; j++) 
{
u = i * k + j; 
c[u] = 0.0; 
for(l = 0; l <= n-1; l++)
{
c[u] = c[u] + a[i * n + l] * b[l * k + j]; 
}
}
}
return; 
}

// 测试代码：calculate.cpp

#include <iostream>
#include "matrixcomputation.h"

using namespace std;

int main() 
{ 
const int NUM = 8;
int i, j; 

static double a[NUM][NUM]=
{ 
{16,11,10,16,24,40,51,61},
{12,12,14,19,26,58,60,55},
{14,13,16,24,40,57,69,56},
{14,17,22,29,51,87,80,62},
{18,22,37,56,68,109,103,77},
{24,35,55,64,81,104,113,92},
{49,64,78,87,103,121,120,101},
{72,92,95,98,112,100,103,99}
}; 

double b[NUM][NUM], c[NUM][NUM]; 

// b是a的一份拷贝
memcpy(b, a, NUM * NUM * sizeof(double));

// 计算逆矩阵
i=inverseMatrix(&a[0][0],NUM);

if (i!=0) 
{ 
cout << "Original matrix is:" << endl; 
for(i = 0; i < NUM; i++) 
{ 
for(j = 0; j < NUM; j++)
{
printf("%7.2f ",b[i][j]); 
}
cout << endl;
} 
cout << endl;

printf("Inverse matrix of the original matrix is:/n"); 
for(i = 0; i < NUM; i++) 
{ 
for(j = 0; j < NUM; j++)
{
printf("%7.2f ",a[i][j]); 
}
cout << endl;
} 
cout << endl;

// 输出原矩阵及其逆矩阵的乘积
printf("Multiplication of the original matrix and its inverse matrix is:/n"); 
multiplyMatrix(&b[0][0],&a[0][0],NUM,NUM,NUM,&c[0][0]); 
for(i = 0; i < NUM; i++) 
{ 
for(j = 0; j < NUM; j++)
{
printf("%7.2f ",c[i][j]); 
}
printf("/n");
} 
} 

return 0; 
} 

运算结果：
/*
Original matrix is:
 16.00 11.00 10.00 16.00 24.00 40.00 51.00 61.00 
 12.00 12.00 14.00 19.00 26.00 58.00 60.00 55.00 
 14.00 13.00 16.00 24.00 40.00 57.00 69.00 56.00 
 14.00 17.00 22.00 29.00 51.00 87.00 80.00 62.00 
 18.00 22.00 37.00 56.00 68.00 109.00 103.00 77.00 
 24.00 35.00 55.00 64.00 81.00 104.00 113.00 92.00 
 49.00 64.00 78.00 87.00 103.00 121.00 120.00 101.00 
 72.00 92.00 95.00 98.00 112.00 100.00 103.00 99.00 

Inverse matrix of the original matrix is:
 0.10 -0.12 0.14 -0.07 -0.16 -0.27 0.64 -0.31 
 -0.16 0.20 -0.09 0.07 0.18 0.23 -0.72 0.37 
 0.06 -0.06 0.02 -0.03 -0.21 -0.11 0.53 -0.28 
 -0.02 0.04 -0.03 -0.05 0.15 0.05 -0.25 0.13 
 0.04 -0.11 -0.00 0.07 -0.00 0.02 -0.04 0.02 
 0.01 0.00 -0.05 0.03 -0.00 -0.03 0.05 -0.03 
 -0.07 0.07 0.11 -0.06 -0.02 -0.02 0.04 -0.02 
 0.05 -0.01 -0.08 0.03 0.03 0.06 -0.11 0.05 

Multiplication of the original matrix and its inverse matrix is:
 1.00 -0.00 0.00 0.00 0.00 -0.00 0.00 -0.00 
 -0.00 1.00 0.00 0.00 0.00 0.00 -0.00 -0.00 
 0.00 -0.00 1.00 0.00 0.00 -0.00 0.00 -0.00 
 0.00 -0.00 0.00 1.00 -0.00 -0.00 0.00 -0.00 
 -0.00 -0.00 0.00 -0.00 1.00 -0.00 0.00 -0.00 
 0.00 -0.00 0.00 0.00 -0.00 1.00 0.00 -0.00 
 -0.00 -0.00 0.00 0.00 0.00 -0.00 1.00 -0.00 
 -0.00 -0.00 -0.00 -0.00 0.00 0.00 0.00 1.00 
*/