C语言科学计算入门之矩阵乘法的相关计算_C 语言

1.矩阵相乘
矩阵相乘应满足的条件:
(1) 矩阵A的列数必须等于矩阵B的行数,矩阵A与矩阵B才能相乘;
(2) 矩阵C的行数等于矩阵A的行数,矩阵C的列数等于矩阵B的列数;
(3) 矩阵C中第i行第j列的元素等于矩阵A的第i行元素与矩阵B的第j列元素对应乘积之和,即

如:

则:

2. 常用矩阵相乘算法
    用A的第i行分别和B的第j列的各个元素相乘求和,求得C的第i行j列的元素,这种算法中,B的访问是按列进行访问的,代码如下:

void arymul(int a[4][5], int b[5][3], int c[4][3])
{
 int i, j, k;
 int temp;
 for(i = 0; i < 4; i++){
 for(j = 0; j < 3; j++){
  temp = 0;
  for(k = 0; k < 5; k++){
  temp += a[i][k] * b[k][j];
  }
  c[i][j] = temp;
  printf("%d/t", c[i][j]);
 }
 printf("%d/n");
 }
}

3. 改进的算法
    矩阵A、B、C都按行(数据的存储顺序)访问,以提高存储器访问效率,对于A的第i行中,第j列的元素分别和B的第j行的元素相乘,对于B中相同的列k在上述计算过程中求和,从而得到C第i行k列的数据,代码如下:

void arymul1(int a[4][5], int b[5][3], int c[4][3])
{
 int i, j, k;
 int temp[3] = {0};
 for(i = 0; i < 4; i++){
 for(k = 0; k < 3; k ++)
  temp[k] = 0;
 for(j = 0; j < 5; j++){//当前行的每个元素
  for(k = 0; k < 3; k++){
  temp[k] += a[i][j] * b[j][k];
  }
 }
 for(k = 0; k < 3; k++){
  c[i][k] = temp[k];
  printf("%d/t", c[i][k]);
 }
 printf("%d/n");
 }
}

这种算法很容易转到稀疏矩阵的相乘算法。

PS:斯特拉森算法的实现
斯特拉森方法,是由v.斯特拉森在1969年提出的一个方法。

我们先讨论二阶矩阵的计算方法。
对于二阶矩阵

a11 a12 b11 b12
A = a21 a22 B = b21 b22

先计算下面7个量(1)

x1 = (a11 + a22) * (b11 + b22);
x2 = (a21 + a22) * b11;
x3 = a11 * (b12 - b22);
x4 = a22 * (b21 - b11);
x5 = (a11 + a12) * b22;
x6 = (a21 - a11) * (b11 + b12);
x7 = (a12 - a22) * (b21 + b22);

再设C = AB。根据矩阵相乘的规则,C的各元素为(2)

c11 = a11 * b11 + a12 * b21
c12 = a11 * b12 + a12 * b22
c21 = a21 * b11 + a22 * b21
c22 = a21 * b12 + a22 * b22

比较(1)(2),C的各元素可以表示为(3)

c11 = x1 + x4 - x5 + x7
c12 = x3 + x5
c21 = x2 + x4
c22 = x1 + x3 - x2 + x6

根据以上的方法,我们就可以计算4阶矩阵了,先将4阶矩阵A和B划分成四块2阶矩阵,分别利用公式计算它们的乘积,再使用(1)(3)来计算出最后结果。

ma11 ma12 mb11 mb12
A4 = ma21 ma22 B4 = mb21 mb22

其中

a11 a12 a13 a14 b11 b12 b13 b14
ma11 = a21 a22 ma12 = a23 a24 mb11 = b21 b22 mb12 = b23 b24 

a31 a32 a33 a34 b31 b32 b33 b34
ma21 = a41 a42 ma22 = a43 a44 mb21 = b41 b42 mb22 = b43 b44 

实现

// 计算2X2矩阵
void Multiply2X2(float& fOut_11, float& fOut_12, float& fOut_21, float& fOut_22,
float f1_11, float f1_12, float f1_21, float f1_22,
float f2_11, float f2_12, float f2_21, float f2_22)
{
const float x1((f1_11 + f1_22) * (f2_11 + f2_22));
const float x2((f1_21 + f1_22) * f2_11);
const float x3(f1_11 * (f2_12 - f2_22));
const float x4(f1_22 * (f2_21 - f2_11));
const float x5((f1_11 + f1_12) * f2_22);
const float x6((f1_21 - f1_11) * (f2_11 + f2_12));
const float x7((f1_12 - f1_22) * (f2_21 + f2_22));
fOut_11 = x1 + x4 - x5 + x7;
fOut_12 = x3 + x5;
fOut_21 = x2 + x4;
fOut_22 = x1 - x2 + x3 + x6;
}
// 计算4X4矩阵
void Multiply(CLAYMATRIX& mOut, const CLAYMATRIX& m1, const CLAYMATRIX& m2)
{
float fTmp[7][4];
// (ma11 + ma22) * (mb11 + mb22)
Multiply2X2(fTmp[0][0], fTmp[0][1], fTmp[0][2], fTmp[0][3],
m1._11 + m1._33, m1._12 + m1._34, m1._21 + m1._43, m1._22 + m1._44,
m2._11 + m2._33, m2._12 + m2._34, m2._21 + m2._43, m2._22 + m2._44);
// (ma21 + ma22) * mb11
Multiply2X2(fTmp[1][0], fTmp[1][1], fTmp[1][2], fTmp[1][3],
m1._31 + m1._33, m1._32 + m1._34, m1._41 + m1._43, m1._42 + m1._44,
m2._11, m2._12, m2._21, m2._22);
// ma11 * (mb12 - mb22)
Multiply2X2(fTmp[2][0], fTmp[2][1], fTmp[2][2], fTmp[2][3],
m1._11, m1._12, m1._21, m1._22,
m2._13 - m2._33, m2._14 - m2._34, m2._23 - m2._43, m2._24 - m2._44);
// ma22 * (mb21 - mb11)
Multiply2X2(fTmp[3][0], fTmp[3][1], fTmp[3][2], fTmp[3][3],
m1._33, m1._34, m1._43, m1._44,
m2._31 - m2._11, m2._32 - m2._12, m2._41 - m2._21, m2._42 - m2._22);
// (ma11 + ma12) * mb22
Multiply2X2(fTmp[4][0], fTmp[4][1], fTmp[4][2], fTmp[4][3],
m1._11 + m1._13, m1._12 + m1._14, m1._21 + m1._23, m1._22 + m1._24,
m2._33, m2._34, m2._43, m2._44);
// (ma21 - ma11) * (mb11 + mb12)
Multiply2X2(fTmp[5][0], fTmp[5][1], fTmp[5][2], fTmp[5][3],
m1._31 - m1._11, m1._32 - m1._12, m1._41 - m1._21, m1._42 - m1._22,
m2._11 + m2._13, m2._12 + m2._14, m2._21 + m2._23, m2._22 + m2._24);
// (ma12 - ma22) * (mb21 + mb22)
Multiply2X2(fTmp[6][0], fTmp[6][1], fTmp[6][2], fTmp[6][3],
m1._13 - m1._33, m1._14 - m1._34, m1._23 - m1._43, m1._24 - m1._44,
m2._31 + m2._33, m2._32 + m2._34, m2._41 + m2._43, m2._42 + m2._44);
// 第一块
mOut._11 = fTmp[0][0] + fTmp[3][0] - fTmp[4][0] + fTmp[6][0];
mOut._12 = fTmp[0][1] + fTmp[3][1] - fTmp[4][1] + fTmp[6][1];
mOut._21 = fTmp[0][2] + fTmp[3][2] - fTmp[4][2] + fTmp[6][2];
mOut._22 = fTmp[0][3] + fTmp[3][3] - fTmp[4][3] + fTmp[6][3];
// 第二块
mOut._13 = fTmp[2][0] + fTmp[4][0];
mOut._14 = fTmp[2][1] + fTmp[4][1];
mOut._23 = fTmp[2][2] + fTmp[4][2];
mOut._24 = fTmp[2][3] + fTmp[4][3];
// 第三块
mOut._31 = fTmp[1][0] + fTmp[3][0];
mOut._32 = fTmp[1][1] + fTmp[3][1];
mOut._41 = fTmp[1][2] + fTmp[3][2];
mOut._42 = fTmp[1][3] + fTmp[3][3];
// 第四块
mOut._33 = fTmp[0][0] - fTmp[1][0] + fTmp[2][0] + fTmp[5][0];
mOut._34 = fTmp[0][1] - fTmp[1][1] + fTmp[2][1] + fTmp[5][1];
mOut._43 = fTmp[0][2] - fTmp[1][2] + fTmp[2][2] + fTmp[5][2];
mOut._44 = fTmp[0][3] - fTmp[1][3] + fTmp[2][3] + fTmp[5][3];
}

比较
在标准的定义算法中我们需要进行n * n * n次乘法运算,新算法中我们需要进行7log2n次乘法,对于最常用的4阶矩阵:   原算法 新算法
加法次数 48 72(48次加法,24次减法)
乘法次数 64 49
需要额外空间 16 * sizeof(float) 28 * sizeof(float)
新算法要比原算法多了24次减法运算,少了15次乘法。但因为浮点乘法的运算速度要远远慢于加/减法运算,所以新算法的整体速度有所提高。

以上是小编为您精心准备的的内容,在的博客、问答、公众号、人物、课程等栏目也有的相关内容,欢迎继续使用右上角搜索按钮进行搜索c语言
, 矩阵
, 入门
C语言入门
c语言矩阵乘法、c语言实现矩阵乘法、c语言矩阵乘法函数、矩阵的乘法c语言、c语言矩阵乘法优化,以便于您获取更多的相关知识。

时间: 2024-09-27 02:37:57

C语言科学计算入门之矩阵乘法的相关计算_C 语言的相关文章

矩阵乘法的运算量计算(华为OJ)

题目地址:https://www.nowcoder.com/practice/15e41630514445719a942e004edc0a5b?tpId=37&&tqId=21293&rp=1&ru=/activity/oj&qru=/ta/huawei/question-ranking 题目内容 矩阵乘法的运算量与矩阵乘法的顺序强相关. 例如: A是一个50×10的矩阵,B是10×20的矩阵,C是20×5的矩阵 计算ABC有两种顺序:((AB)C)或者(A(BC)

深入解析C语言中的内存分配相关问题_C 语言

C内存分配区域 程序代码区存放函数体的二进制代码 全局数据区全局变量和静态变量的存储是放在一起的.初始化的全局变量和静态变量在一块区域,未初始化的全局变量和未初始化的静态变量在相邻的另一块区域.常量数据存放在另一个区域里.这些数据在程序结束后由系统释放.我们所说的BSS段(bss segment)通常是指用来存放程序中未初始化的全局变量的一块内存区域.BSS是英文Block Started by Symbol的简称 栈区由编译器自动分配释放,存放函数的参数值,局部变量的值等.其操作方式类似于数据

C语言中操作进程信号的相关函数使用详解_C 语言

C语言signal()函数:设置信号处理方式头文件: #include <signal.h> 定义函数: void (*signal(int signum, void(* handler)(int)))(int); 函数说明:signal()会依参数signum 指定的信号编号来设置该信号的处理函数. 当指定的信号到达时就会跳转到参数handler 指定的函数执行. 如果参数handler 不是函数指针, 则必须是下列两个常数之一: 1.SIG_IGN 忽略参数signum 指定的信号. 2.

利用简洁的C语言代码解决跳台阶问题与约瑟夫环问题_C 语言

跳台阶问题 题目: 一个台阶总共有 n 级,如果一次可以跳 1 级,也可以跳 2 级. 求总共有多少总跳法,并分析算法的时间复杂度. 分析: 也是比较基础的题目,通过递归可以方便的求解 代码实现如下(GCC编译通过): #include "stdio.h" #include "stdlib.h" int function(int n); int main(void) { int tmp; tmp = function(5); printf("%3d\n&q

c语言中数组名a和&amp;amp;a详细介绍_C 语言

最近又把学习c语言提上日程上来了~~~先把我打算看的书都写下来吧,<C语言深度剖析>,<c和指针>系类,<c语言陷阱和缺陷> 先说说a和&a的区别(有三点,三个方向):1.是a和&a的本质,都是什么类型的.2.从2维数组的角度看.3.从指针运算的角度看. 声明:虽然数组名不是指针,但是用的很像指针,我们暂且把它叫做一个指针吧. 第一个问题:int a[10];  a ,&a和&a[0] 都是分别是什么?先说明a ,&a和&

C语言中的数组和指针汇编代码分析实例_C 语言

今天看<程序员面试宝典>时偶然看到讲数组和指针的存取效率,闲着无聊,就自己写了段小代码,简单分析一下C语言背后的汇编,可能很多人只注重C语言,但在实际应用当中,当出现问题时,有时候还是通过分析汇编代码能够解决问题.本文只是为初学者,大牛可以飘过~ C源代码如下: 复制代码 代码如下: #include "stdafx.h" int main(int argc, char* argv[]) {        char a=1;        char c[] = "

C语言安全编码之数组索引位的合法范围_C 语言

C语言中的数组索引必须保证位于合法的范围内! 示例代码如下: enum {TABLESIZE = 100}; int *table = NULL; int insert_in_table(int pos, int value) { if(!table) { table = (int *)malloc(sizeof(int) *TABLESIZE); } if(pos >= TABLESIZE) { return -1; } table[pos] = value; return 0; } 其中:p

举例讲解C语言程序中对二叉树数据结构的各种遍历方式_C 语言

二叉树遍历的基本思想 二叉树的遍历本质上其实就是入栈出栈的问题,递归算法简单且容易理解,但是效率始终是个问题.非递归算法可以清楚的知道每步实现的细节,但是乍一看不想递归算法那么好理解,各有各的好处吧.接下来根据下图讲讲树的遍历. 1.先序遍历:先序遍历是先输出根节点,再输出左子树,最后输出右子树.上图的先序遍历结果就是:ABCDEF  2.中序遍历:中序遍历是先输出左子树,再输出根节点,最后输出右子树.上图的中序遍历结果就是:CBDAEF 3.后序遍历:后序遍历是先输出左子树,再输出右子树,最后

C语言中结构体(struct)的几种初始化方法_C 语言

本文给大家总结的struct数据有3种初始化方法      1.顺序      2.C风格的乱序      3.C++风格的乱序 下面通过示例代码详细介绍这三种初始化方法. 1)顺序 这种方法很常见,在一般的介绍C的书中都有介绍.顺序初始化的特点是: 按照成员定义的顺序,从前到后逐个初始化:允许只初始化部分成员:在被初始化的成员之前,不能有未初始化的成员. 示例: struct User oneUser = {10, "Lucy", "/home/Lucy"}; 2