矩阵乘法

最新推荐文章于 2024-06-09 22:45:40 发布
最新推荐文章于 2024-06-09 22:45:40 发布
阅读量1.6k 收藏 13
点赞数 5
分类专栏： 算法导论 文章标签： 矩阵乘法

   编程实现矩阵乘法，并考虑当矩阵规模较大时的优化方法:

(1)一般的解法:

 

//矩阵乘法，3个for循环搞定   
void Mul(int** matrixA, int** matrixB, int** matrixC)   
{   
    for(int i = 0; i < n; ++i)    
    {   
        for(int j = 0; j < n; ++j)    
        {   
            matrixC[i][j] = 0;   
            for(int k = 0; k < n; ++k)    
            {   
                matrixC[i][j] += matrixA[i][k] * matrixB[k][j];   
            }   
        }   
    }   
}
(2)strassen算法(分块+分治思想):


算法思想:



strassen算法算法流程:



通用矩阵相乘算法时间复杂度是技术分享，而Strassen算法复杂度只是技术分享。但随着n的变大，比如当n >> 100时，Strassen算法是比通用矩阵相乘算法变得更有效率。

具体的代码实现如下:

#include <iostream>
#include <ctime>
#include <Windows.h>
using namespace std;


template<typename T>
class Strassen_class{
public:
	void ADD(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize);  //矩阵相加
	void SUB(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize);  //矩阵相减
	void MUL(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize);//朴素算法实现
	void FillMatrix(T** MatrixA, T** MatrixB, int length);//A,B矩阵赋值
	void PrintMatrix(T **MatrixA, int MatrixSize);//打印矩阵
	void Strassen(int N, T **MatrixA, T **MatrixB, T **MatrixC);//Strassen算法实现
};


//矩阵的相加
template<typename T>
void Strassen_class<T>::ADD(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize)
{
	for (int i = 0; i < MatrixSize; i++)
	{
		for (int j = 0; j < MatrixSize; j++)
		{
			MatrixResult[i][j] = MatrixA[i][j] + MatrixB[i][j];
		}
	}
}

// 矩阵的相减
template<typename T>
void Strassen_class<T>::SUB(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize)
{
	for (int i = 0; i < MatrixSize; i++)
	{
		for (int j = 0; j < MatrixSize; j++)
		{
			MatrixResult[i][j] = MatrixA[i][j] - MatrixB[i][j];
		}
	}
}
// 矩阵的相乘
template<typename T>
void Strassen_class<T>::MUL(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize)
{
	for (int i = 0; i < MatrixSize; i++)
	{
		for (int j = 0; j < MatrixSize; j++)
		{
			MatrixResult[i][j] = 0;
			for (int k = 0; k < MatrixSize; k++)
			{
				MatrixResult[i][j] = MatrixResult[i][j] + MatrixA[i][k] * MatrixB[k][j];
			}
		}
	}
}

/*
c++使用二维数组，申请动态内存方法
申请
int **A;
A = new int *[desired_array_row];
for ( int i = 0; i < desired_array_row; i++)
A[i] = new int [desired_column_size];

释放
for ( int i = 0; i < your_array_row; i++)
delete [] A[i];
delete[] A;

*/
template<typename T>
void Strassen_class<T>::Strassen(int N, T **MatrixA, T **MatrixB, T **MatrixC)
{

	int HalfSize = N / 2;
	int newSize = N / 2;

	if (N <= 64)    //分治门槛，小于这个值时不再进行递归计算，而是采用常规矩阵计算方法
	{
		MUL(MatrixA, MatrixB, MatrixC, N);
	}
	else
	{
		T** A11;
		T** A12;
		T** A21;
		T** A22;

		T** B11;
		T** B12;
		T** B21;
		T** B22;

		T** C11;
		T** C12;
		T** C21;
		T** C22;

		T** M1;
		T** M2;
		T** M3;
		T** M4;
		T** M5;
		T** M6;
		T** M7;
		T** AResult;
		T** BResult;

		//making a 1 diminsional pointer based array.
		A11 = new T *[newSize];
		A12 = new T *[newSize];
		A21 = new T *[newSize];
		A22 = new T *[newSize];

		B11 = new T *[newSize];
		B12 = new T *[newSize];
		B21 = new T *[newSize];
		B22 = new T *[newSize];

		C11 = new T *[newSize];
		C12 = new T *[newSize];
		C21 = new T *[newSize];
		C22 = new T *[newSize];

		M1 = new T *[newSize];
		M2 = new T *[newSize];
		M3 = new T *[newSize];
		M4 = new T *[newSize];
		M5 = new T *[newSize];
		M6 = new T *[newSize];
		M7 = new T *[newSize];

		AResult = new T *[newSize];
		BResult = new T *[newSize];

		int newLength = newSize;

		//making that 1 diminsional pointer based array , a 2D pointer based array
		for (int i = 0; i < newSize; i++)
		{
			A11[i] = new T[newLength];
			A12[i] = new T[newLength];
			A21[i] = new T[newLength];
			A22[i] = new T[newLength];

			B11[i] = new T[newLength];
			B12[i] = new T[newLength];
			B21[i] = new T[newLength];
			B22[i] = new T[newLength];

			C11[i] = new T[newLength];
			C12[i] = new T[newLength];
			C21[i] = new T[newLength];
			C22[i] = new T[newLength];

			M1[i] = new T[newLength];
			M2[i] = new T[newLength];
			M3[i] = new T[newLength];
			M4[i] = new T[newLength];
			M5[i] = new T[newLength];
			M6[i] = new T[newLength];
			M7[i] = new T[newLength];

			AResult[i] = new T[newLength];
			BResult[i] = new T[newLength];


		}
		//splitting input Matrixes, into 4 submatrices each.
		for (int i = 0; i < N / 2; i++)
		{
			for (int j = 0; j < N / 2; j++)
			{
				A11[i][j] = MatrixA[i][j];
				A12[i][j] = MatrixA[i][j + N / 2];
				A21[i][j] = MatrixA[i + N / 2][j];
				A22[i][j] = MatrixA[i + N / 2][j + N / 2];

				B11[i][j] = MatrixB[i][j];
				B12[i][j] = MatrixB[i][j + N / 2];
				B21[i][j] = MatrixB[i + N / 2][j];
				B22[i][j] = MatrixB[i + N / 2][j + N / 2];

			}
		}

		//here we calculate M1..M7 matrices .
		//M1[][]
		ADD(A11, A22, AResult, HalfSize);
		ADD(B11, B22, BResult, HalfSize);                //p5=(a+d)*(e+h)
		Strassen(HalfSize, AResult, BResult, M1); //now that we need to multiply this , we use the strassen itself .


		//M2[][]
		ADD(A21, A22, AResult, HalfSize);              //M2=(A21+A22)B11   p3=(c+d)*e
		Strassen(HalfSize, AResult, B11, M2);       //Mul(AResult,B11,M2);

		//M3[][]
		SUB(B12, B22, BResult, HalfSize);              //M3=A11(B12-B22)   p1=a*(f-h)
		Strassen(HalfSize, A11, BResult, M3);       //Mul(A11,BResult,M3);

		//M4[][]
		SUB(B21, B11, BResult, HalfSize);           //M4=A22(B21-B11)    p4=d*(g-e)
		Strassen(HalfSize, A22, BResult, M4);       //Mul(A22,BResult,M4);

		//M5[][]
		ADD(A11, A12, AResult, HalfSize);           //M5=(A11+A12)B22   p2=(a+b)*h
		Strassen(HalfSize, AResult, B22, M5);       //Mul(AResult,B22,M5);


		//M6[][]
		SUB(A21, A11, AResult, HalfSize);
		ADD(B11, B12, BResult, HalfSize);             //M6=(A21-A11)(B11+B12)   p7=(c-a)(e+f)
		Strassen(HalfSize, AResult, BResult, M6);    //Mul(AResult,BResult,M6);

		//M7[][]
		SUB(A12, A22, AResult, HalfSize);
		ADD(B21, B22, BResult, HalfSize);             //M7=(A12-A22)(B21+B22)    p6=(b-d)*(g+h)
		Strassen(HalfSize, AResult, BResult, M7);     //Mul(AResult,BResult,M7);

		//C11 = M1 + M4 - M5 + M7;
		ADD(M1, M4, AResult, HalfSize);
		SUB(M7, M5, BResult, HalfSize);
		ADD(AResult, BResult, C11, HalfSize);

		//C12 = M3 + M5;
		ADD(M3, M5, C12, HalfSize);

		//C21 = M2 + M4;
		ADD(M2, M4, C21, HalfSize);

		//C22 = M1 + M3 - M2 + M6;
		ADD(M1, M3, AResult, HalfSize);
		SUB(M6, M2, BResult, HalfSize);
		ADD(AResult, BResult, C22, HalfSize);

		//at this point , we have calculated the c11..c22 matrices, and now we are going to
		//put them together and make a unit matrix which would describe our resulting Matrix.
		//组合小矩阵到一个大矩阵
		for (int i = 0; i < N / 2; i++)
		{
			for (int j = 0; j < N / 2; j++)
			{
				MatrixC[i][j] = C11[i][j];
				MatrixC[i][j + N / 2] = C12[i][j];
				MatrixC[i + N / 2][j] = C21[i][j];
				MatrixC[i + N / 2][j + N / 2] = C22[i][j];
			}
		}

		// 释放矩阵内存空间
		for (int i = 0; i < newLength; i++)
		{
			delete[] A11[i]; delete[] A12[i]; delete[] A21[i];
			delete[] A22[i];

			delete[] B11[i]; delete[] B12[i]; delete[] B21[i];
			delete[] B22[i];
			delete[] C11[i]; delete[] C12[i]; delete[] C21[i];
			delete[] C22[i];
			delete[] M1[i]; delete[] M2[i]; delete[] M3[i]; delete[] M4[i];
			delete[] M5[i]; delete[] M6[i]; delete[] M7[i];
			delete[] AResult[i]; delete[] BResult[i];
		}
		delete[] A11; delete[] A12; delete[] A21; delete[] A22;
		delete[] B11; delete[] B12; delete[] B21; delete[] B22;
		delete[] C11; delete[] C12; delete[] C21; delete[] C22;
		delete[] M1; delete[] M2; delete[] M3; delete[] M4; delete[] M5;
		delete[] M6; delete[] M7;
		delete[] AResult;
		delete[] BResult;

	}//end of else

}

//填充矩阵
template<typename T>
void Strassen_class<T>::FillMatrix(T** MatrixA, T** MatrixB, int length)
{
	for (int row = 0; row < length; row++)
	{
		for (int column = 0; column < length; column++)
		{

			MatrixB[row][column] = (MatrixA[row][column] = rand() % 5);

		}

	}
}

//打印矩阵
template<typename T>
void Strassen_class<T>::PrintMatrix(T **MatrixA, int MatrixSize)
{
	cout << endl;
	for (int row = 0; row < MatrixSize; row++)
	{
		for (int column = 0; column < MatrixSize; column++)
		{


			cout << MatrixA[row][column] << "\t";
			if ((column + 1) % ((MatrixSize)) == 0)
				cout << endl;
		}

	}
	cout << endl;
}

int main()
{
	Strassen_class<int> stra;//定义Strassen_class类对象
	int MatrixSize = 0;

	int** MatrixA;    //存放矩阵A
	int** MatrixB;    //存放矩阵B
	int** MatrixC;    //存放结果矩阵

	clock_t startTime_For_Normal_Multipilication;
	clock_t endTime_For_Normal_Multipilication;

	clock_t startTime_For_Strassen;
	clock_t endTime_For_Strassen;
	srand(time(0));

	cout << "\n请输入矩阵大小(必须是2的幂指数值(例如:32,64,512,..): ";
	cin >> MatrixSize;

	int N = MatrixSize;//for readiblity.

	//申请内存
	MatrixA = new int *[MatrixSize];
	MatrixB = new int *[MatrixSize];
	MatrixC = new int *[MatrixSize];

	for (int i = 0; i < MatrixSize; i++)
	{
		MatrixA[i] = new int[MatrixSize];
		MatrixB[i] = new int[MatrixSize];
		MatrixC[i] = new int[MatrixSize];
	}

	stra.FillMatrix(MatrixA, MatrixB, MatrixSize);  //矩阵赋值

	//*******************conventional multiplication test
	cout << "朴素矩阵算法开始时钟:  " << (startTime_For_Normal_Multipilication = clock());

	stra.MUL(MatrixA, MatrixB, MatrixC, MatrixSize);//朴素矩阵相乘算法 T(n) = O(n^3)

	cout << "\n朴素矩阵算法结束时钟: " << (endTime_For_Normal_Multipilication = clock());

	cout << "\n矩阵运算结果... \n";
	stra.PrintMatrix(MatrixC, MatrixSize);

	//*******************Strassen multiplication test
	cout << "\nStrassen算法开始时钟: " << (startTime_For_Strassen = clock());

	stra.Strassen(N, MatrixA, MatrixB, MatrixC); //strassen矩阵相乘算法

	cout << "\nStrassen算法结束时钟: " << (endTime_For_Strassen = clock());


	cout << "\n矩阵运算结果... \n";
	stra.PrintMatrix(MatrixC, MatrixSize);

	cout << "矩阵大小 " << MatrixSize;
	cout << "\n朴素矩阵算法: " << (endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication) << " Clocks.." << (endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication) / CLOCKS_PER_SEC << " Sec";
	cout << "\nStrassen算法:" << (endTime_For_Strassen - startTime_For_Strassen) << " Clocks.." << (endTime_For_Strassen - startTime_For_Strassen) / CLOCKS_PER_SEC << " Sec\n";
	system("Pause");

	return 0;
}

运行结果如下:


三名狂客
关注
  • 5
    点赞
  • 踩
  • 13
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
矩阵的相乘
weixin_44274975的博客
03-15 152
def matrix_multi(m1,m2): result = np.zeros((m1.shape[0],m2.shape[1])) if m1.shape[1] != m2.shape[0]: return #遍历第一个矩阵的行 for i in range(m1.shape[0]): #遍历第二个矩阵的列 for j in range(m2.shape[1]): sum = 0 for ...
python矩阵和向量乘法总结
QWxixi的博客
08-30 5986
python中常见矩阵和向量乘法总结
【架构分析】GPU执行GEMM矩阵运算实例演示
最新发布
HaoBBNuanMM的博客
06-09 465
Cutlass是 NVIDIA 提供的一套用于高效实现矩阵乘法和卷积操作的 C++ 库。它以 CUDA 为基础，提供了高度优化的数学运算，尤其适用于GPU上的高性能并行计算。本文以GEMM矩阵运算作为实例，展示Cutlass在GPU上执行GEMM运算的过程。
GEMM算法及优化流程详解
独行侠
02-16 2万+
前言 神经网络前向耗时主要由卷积的耗时决定，参考賈杨青毕业论文，那么如何对卷积加速便成了重要的一个点，主流的加速方法有 以下几种： im2col+GEMM：目前几乎所有的主流计算框架包括 Caffe, MXNet 等都实现了该方法. 该方法把整个卷积过程转化成了GEMM过程，而GEMM在各种 BLAS 库中都是被极致优化的，一般来说，速度较快。 Winograd： Winograd 是存在...
【矩阵计算】矩阵乘法其一：基础符号和算法
一直学，你始终是学生，但是只要创造一招，你便是宗师。
09-22 5320
矩阵符号矩阵操作向量符号向量操作Saxpy算法Gaxpy算法外积矩阵分割和冒号符号矩阵-矩阵乘法复数矩阵 矩阵符号 如果用表示所有实数的集合，那么我们用表示所有的实数矩阵组成的向量空间，即： 其中，大写字母（如）表示矩阵，带下标的小写字母（如）表示矩阵中的元素。除了用表示矩阵中第行第列的元素之外，也可以用和表示。 矩阵操作 矩阵转置（transposition）: 矩阵加法（addition...
【线性代数】矩阵的三种相乘方式
小白兔的窝
07-24 3万+
普通乘积 矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的列数（column）和第二个矩阵的行数（row）相同时才有意义[1]  。一般单指矩阵乘积时，指的便是一般矩阵乘积。一个m×n的矩阵就是m×n个数排成m行n列的一个数阵。由于它把许多数据紧凑的集中到了一起，所以有时候可以简便地表示一些复杂的模型。 设A为    的矩阵，B为    的矩阵，那么称  
Java实现的矩阵乘法示例
08-26
Java实现的矩阵乘法示例 矩阵乘法是线性代数中的一种基本操作，它将两个矩阵相乘，得到一个新的矩阵。在Java中，矩阵乘法可以通过循环遍历和数组操作来实现。下面是Java实现的矩阵乘法示例的知识点总结： 一、矩阵...
并行矩阵乘法
09-24
并行矩阵乘法是计算机科学中的一个重要概念，特别是在高性能计算和大数据处理领域。它通过将一个大任务分解成多个小任务，然后在多个处理器或核心上同时执行这些小任务，来显著提高计算效率。在本项目中，我们将讨论...
OpenBLAS项目与矩阵乘法优化_OpenBLAS项目与矩阵乘法优化_openblas_
09-30
OpenBLAS，全称为Open Basic Linear Algebra Subprograms，是一个开源的高性能线性代数库，专为基本线性代数操作而设计，如矩阵乘法。这个项目旨在提供一个免费且可高度优化的替代方案，以提升计算密集型应用的性能...
矩阵乘法MPI并行程序报告.pdf
06-06
【矩阵乘法MPI并行程序】是分布式计算中的一种常见应用，主要利用Message Passing Interface (MPI)库来实现。MPI是一种标准，允许程序员在多处理器系统或者计算机集群上编写并行程序，它提供了丰富的通信原语，如MPI...
矩阵乘法_矩阵乘法_
09-30
在计算机科学和数学中，矩阵乘法是一种基本的操作，它涉及到二维数组，即矩阵的运算。这个主题在编程中尤其重要，特别是在数值计算、图形学、机器学习等领域。本篇文章将详细探讨矩阵乘法的概念、算法以及如何用...
BLAS之GEMM和GEMV
永远飞翔的鸟
10-17 7772
BLAS简介 BLAS全称是Basic Linear Algebra Subprograms是规定了一套低级的执行常见线性代数操作的规范。其实现经常针对特殊的机器进行优化，比较著名的·BLAS库有ACML, ATLAS, MKL, OpenBLAS。许多常见的数值软件均采用兼容BLAS规范的实现库来进行线性代数计算，比如Matlab, Numpy, Mathematica`。 其中，Level 1 BLAS主要提供向量操作 Level 2 BLAS提供矩阵向量操作（gemv) gemv
GEMM矩阵相乘与深度学习
weixin_43728590的博客
07-08 4211
GEMM矩阵相乘与深度学习1. GEMM矩阵相乘1.1 GEMM算法基础1.2 GEMM算法优化1.2.1 循环重排充分利用缓存1.2.2 平铺（Tiling）充分利用缓存1.2.3 展开（Unrolling）1.2.4 内存对齐1.2.5 向量化1.2.6 矩阵分块1.2.7 双缓冲2. GEMM与卷积计算2.1 计算卷积的方法2.2 Img2col2.2.1 CNN中张量的存储2.2.2 卷积运算转化为GEMM 1. GEMM矩阵相乘 1.1 GEMM算法基础   GEMM（General Matrix
矩阵与向量的乘积
weixin_30528371的博客
11-07 772
以下内容来源于：https://www.zhihu.com/people/August_666/posts 先上运算，再解读： 一个矩阵乘以一个列向量相当于矩阵的列向量的线性组合。 一个行向量乘以矩阵，相当于矩阵的行向量的线性组合。 方程组： 在二维平面中，相当于找两条直线的交点。 写成如下形式： 把方程组看成是Ax=b，相当于是寻找矩阵A的...
Python 矩阵与矩阵以及矩阵与向量的乘法
weixin_34392843的博客
09-25 6762
  import numpy as np numpy模块的array相乘时，有两种方式：一是矩阵形式，二是挨个相乘。 需要用矩阵形式相乘时，则要用np.dot()函数。 #矩阵与矩阵相乘a = np.array([[1,2,3],[4,5,6],[7,8,9]])c = a.copy()print(a * c)print(np.dot(a, c))#a*c 得出的结果是a和c中每个元素依次相乘，...
MKL 函数cblas_?gemm参数分析
u013806541的博客
04-03 3604
参考链接：https://software.intel.com/en-us/node/520775 函数原型： void cblas_sgemm (const CBLAS_LAYOUT Layout, const CBLAS_TRANSPOSE transa, constCBLAS_TRANSPOSE transb, const MKL_INT m, const MKL_INT n, cons
什么是GEMM?该怎么去学习GEMM?
YDM6211的博客
11-22 821
而且猿代码科技也与各个地区的国家超算中心有紧密的合作，例如“无锡国家超算中心”、“天津国家超算中心”等，现向全社会人员以及各个高校在校生提供100个名额，如果你如果你想挣到程序员的高薪，又想有稳定的发展，那请点击下方链接，详细咨询，前50名可得先进计算基础课程，提供面试机会，助你进入目标单位。在同一个CPU循环中，SIMD可在多个值上同时执行相同的运算/指令。我们需要了解一下硬件的基础的体系结构，了解硬件如何工作，才能在软件上实现算法的时候，尽可能去考虑利用硬件的一些特性，来做到高效的、极致的优化。
带你深入理解矩阵乘法
热门推荐
zh_94的博客
08-11 7万+
为了不浪费大家宝贵的时间，开头我先简要说明一下这篇博文对哪些读者可能会有帮助 1、你是正在学习矩阵的乘法运算，觉得矩阵的乘法掌握起来很困难 2、你已经学会了矩阵乘法，但如果你在计算矩阵乘法时还在使用“一行乘一列得一数”的方法，那我强烈建议你看看后面的内容。 因为，我将带你更加深刻地理解矩阵，与之而来是对矩阵乘法的全新计算方式。这不仅让你在计算矩阵乘法时更快，而且更省心。 “矩阵就是数...
Python矩阵乘法
05-08
在Python中，我们可以使用NumPy库来实现矩阵乘法。具体来说，我们可以使用numpy.dot()函数或@运算符实现矩阵乘法。下面是一个简单的示例代码： ``` python import numpy as np # 定义两个矩阵 A = np.array([[1, 2...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值