最近频繁用到矩阵相乘运算, 所以想写一个函数。遇到的问题不少！
关于矩阵乘法规则，请查看： [050] 求两个矩阵的乘积矩阵
算法参考了徐士良的《常用算法程序集》里的一个实矩阵相乘函数，如下：
void brmul(a, b, m, n, k, c)
int m, n, k;
double a[], b[], 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];
}
}
参数说明：
━━━━━━━━━━━━━━━━━━━━━━━━━━━━
形参与函数类型 说明
————————————————————————————
double a[m][n] 存放矩阵A的元素
————————————————————————————
double b[n][k] 存放矩阵B的元素
————————————————————————————
int m 矩阵A与乘积矩阵C的行数
————————————————————————————
int n 矩阵A的列数，矩阵B的行数
————————————————————————————
int k 矩阵B与乘积矩阵C的列数
————————————————————————————
double c[m][k] 返回乘积矩阵C＝AB的元素
━━━━━━━━━━━━━━━━━━━━━━━━━━━━
一个验证程序：
#include <stdio.h>
void brmul(a, b, m, n, k, c)
int m, n, k;
double a[], b[], 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];
}
}
int main()
{
int i, j;
static double a[4][5] = { { 1.0, 3.0, -2.0, 0.0, 4.0},
{-2.0, -1.0, 5.0, -7.0, 2.0},
{ 0.0, 8.0, 4.0, 1.0, -5.0},
{ 3.0, -3.0, 2.0, -4.0, 1.0}};
static double c[4][3], b[5][3] = { { 4.0, 5.0, -1.0},
{ 2.0, -2.0, 6.0},
{ 7.0, 8.0, 1.0},
{ 0.0, 3.0, -5.0},
{ 9.0, 8.0, -6.0}};
brmul(a, b, 4, 5, 3, c);
for (i = 0; i <= 3; i++)
{
for (j = 0; j <= 2; j++)
printf("%14.7e", c[i][j]);
printf("\n");
}
printf("\n");
getch();
return 0;
}
运行结果(TC):
===================================================================
3.200000e+01 1.500000e+01 -9.000000e+00
4.300000e+01 2.700000e+01 2.400000e+01
-1.000000e+00 -2.100000e+01 7.700000e+01
2.900000e+01 3.300000e+01 -5.000000e+00
===================================================================
★ 这个程序可以在TC, Win-TC里运行,但在VC, Dev-C++里都会在函数声明那里报错,于是注意到了这个函数声明部分有些特别, 以前还真没看过这种形式, 参数的类型声明在函数定义和函数体之间,并且这里所用的形式还有不同,比如参数里用的a, 在类型声明时却用的是 double a[] 。想必是这里出了问题，开始怀疑这种方式不是C的标准而是Turbo C里自己定义的，后来查了一下谭老的书才明白，原来这是一种对形参类型声明的传统方式。将书中部分内容摘录如下：
对形参的声明的传统方式：
在老版本的C语言中，对形参类型的声明是放在函数定义的第2行，也就是不在第1行的括号内指定形参的类型，而在括号外单独指定。一般把这种方法称为传统的对形参的声明方式，而把现在常用的这种方法称为现代的方式。Turbo C和目前使用的多数C版本对这两种方法都允许使用，两种用法等价，ANSI新标准推荐后者，即现代方式。它与PASCAL语言中所用的方法是类似的。本书中的程序采用新标准推荐的现代方式。但由于有些过去写的书籍和程序使用传统方式，因此读者应对它有所了解，以便能方便地阅读它们。
又一次感到谭老的书内容还是比较全的。
★ 那么是不是将上面的声明改成现代方式就可以在VC里运行了呢？经试验是不行的，改成如下形式可以在TC里运行：
void brmul(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];
}
}
验证:
#include <stdio.h>
void brmul(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];
}
}
int main()
{
int i, j;
static double a[4][5] = { { 1.0, 3.0, -2.0, 0.0, 4.0},
{-2.0, -1.0, 5.0, -7.0, 2.0},
{ 0.0, 8.0, 4.0, 1.0, -5.0},
{ 3.0, -3.0, 2.0, -4.0, 1.0}};
static double c[4][3], b[5][3] = { { 4.0, 5.0, -1.0},
{ 2.0, -2.0, 6.0},
{ 7.0, 8.0, 1.0},
{ 0.0, 3.0, -5.0},
{ 9.0, 8.0, -6.0}};
brmul(a, b, 4, 5, 3, c);
for (i = 0; i <= 3; i++)
{
for (j = 0; j <= 2; j++)
printf("%14.7e", c[i][j]);
printf("\n");
}
printf("\n");
getch();
return 0;
}
TC运行结果:
===================================================================
3.200000e+01 1.500000e+01 -9.000000e+00
4.300000e+01 2.700000e+01 2.400000e+01
-1.000000e+00 -2.100000e+01 7.700000e+01
2.900000e+01 3.300000e+01 -5.000000e+00
===================================================================
★ 但在VC里调用这个函数时参数传递会出错，即上例中 brmul(a, b, 4, 5, 3, c); 这一句。又在Dev-C++里试了下，也是同样的错误。我想应该是参数传递时出现的问题，即这种参数传递方式也应该类似的是传统的一种方式，而非新标准所推荐的。 记得看过一篇文章讲的是数组名和指针的关系，作者阐述的是数组名并不等同于指针的观点。联系一下，我猜上例中的 brmul(a, b, 4, 5, 3, c); 一句就是将数组名用作指针的，而这种用法是老版本C所支持的，所以可以在TC里运行。而在VC，Dev-C++里执行的应该是新标准，不支持这种用法。那么用指针改写一下是不是就可以了呢? 作如下验证，调用时用 brmul(*a, *b, 4, 5, 3, *c); , 由结果看应该是正确的（由于是在VC里，所以将输出的位宽加了一下）：
#include <stdio.h>
void brmul(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];
}
}
int main()
{
int i, j;
static double a[4][5] = { { 1.0, 3.0, -2.0, 0.0, 4.0},
{-2.0, -1.0, 5.0, -7.0, 2.0},
{ 0.0, 8.0, 4.0, 1.0, -5.0},
{ 3.0, -3.0, 2.0, -4.0, 1.0}};
static double c[4][3], b[5][3] = { { 4.0, 5.0, -1.0},
{ 2.0, -2.0, 6.0},
{ 7.0, 8.0, 1.0},
{ 0.0, 3.0, -5.0},
{ 9.0, 8.0, -6.0}};
brmul(*a, *b, 4, 5, 3, *c);
for (i = 0; i <= 3; i++)
{
for (j = 0; j <= 2; j++)
printf("%16.7e", c[i][j]);
printf("\n");
}
printf("\n");
return 0;
}
运行结果（VC)：
===================================================================
3.2000000e+001 1.5000000e+001 -9.0000000e+000
4.3000000e+001 2.7000000e+001 2.4000000e+001
-1.0000000e+000 -2.1000000e+001 7.7000000e+001
2.9000000e+001 3.3000000e+001 -5.0000000e+000
===================================================================
★ 又在Dev-C++里运行了一下,通过! 最终得出的实矩阵相乘函数为(改了个函数名):
void MatMul(double a[], double b[], int m, int n, int k, double c[])
{
int i, j, t;
int u;
for (i = 0; i<= m - 1; i++)
for (j = 0; j <= k - 1; j++)
{
u = i * k + j;
c[u] = 0.0;
for (t = 0; t <= n - 1; t++)
c[u] = c[u] + a[i * n + t] * b[t * k + j];
}
}
参数说明：
━━━━━━━━━━━━━━━━━━━━━━━━━━━━
形参与函数类型 说明
————————————————————————————
double a[m][n] 存放矩阵A的元素
————————————————————————————
double b[n][k] 存放矩阵B的元素
————————————————————————————
int m 矩阵A与乘积矩阵C的行数
————————————————————————————
int n 矩阵A的列数，矩阵B的行数
————————————————————————————
int k 矩阵B与乘积矩阵C的列数
————————————————————————————
double c[m][k] 返回乘积矩阵C＝AB的元素
━━━━━━━━━━━━━━━━━━━━━━━━━━━━
调用方法: MatMul(*a, *b, m, n, k, *c);
######################################补充########################################
★ 又琢磨了一下, 函数参数中用的如 double a[] 的形式, 既然没有指定维数, 那实际上也就是传了个地址进去, 于是改写为 double *a 的形式:
void MatMul(double *a, double *b, int m, int n, int k, double *c)
{
int i, j, t;
int u;
for (i = 0; i<= m - 1; i++)
for (j = 0; j <= k - 1; j++)
{
u = i * k + j;
c[u] = 0.0;
for (t = 0; t <= n - 1; t++)
c[u] = c[u] + a[i * n + t] * b[t * k + j];
}
}
再次验证：
#include <stdio.h>
void MatMul(double *a, double *b, int m, int n, int k, double *c)
{
int i, j, t;
int u;
for (i = 0; i<= m - 1; i++)
for (j = 0; j <= k - 1; j++)
{
u = i * k + j;
c[u] = 0.0;
for (t = 0; t <= n - 1; t++)
c[u] = c[u] + a[i * n + t] * b[t * k + j];
}
}
int main()
{
int i, j;
static double a[4][5] = { { 1.0, 3.0, -2.0, 0.0, 4.0},
{-2.0, -1.0, 5.0, -7.0, 2.0},
{ 0.0, 8.0, 4.0, 1.0, -5.0},
{ 3.0, -3.0, 2.0, -4.0, 1.0}};
static double c[4][3], b[5][3] = { { 4.0, 5.0, -1.0},
{ 2.0, -2.0, 6.0},
{ 7.0, 8.0, 1.0},
{ 0.0, 3.0, -5.0},
{ 9.0, 8.0, -6.0}};
MatMul(*a, *b, 4, 5, 3, *c);
for (i = 0; i <= 3; i++)
{
for (j = 0; j <= 2; j++)
printf("%16.7e", c[i][j]);
printf("\n");
}
printf("\n");
return 0;
}
运行结果（VC）：
====================================================================
3.2000000e+001 1.5000000e+001 -9.0000000e+000
4.3000000e+001 2.7000000e+001 2.4000000e+001
-1.0000000e+000 -2.1000000e+001 7.7000000e+001
2.9000000e+001 3.3000000e+001 -5.0000000e+000
====================================================================
评论