Openblas加速三维矩阵卷积操作

接上个博客（http://blog.csdn.net/samylee/article/details/73252715）所讲，这篇博文介绍如何用openblas加速三维矩阵。

注意：博主未执行累加运算，只是卷积的过程！

代码如下（以下程序经博主测试准确无误）：

头文件function.h

 

//function.h
//cblas加速三维矩阵卷积操作
//注意：stride=1
//作者：samylee

#ifndef FUNCTION_H
#define FUNCTION_H

#include <cblas.h>  
#include <iostream>

using namespace std;

//A加pad的计算
void comput_Apad(const int pad_w, const int Map, const int channel, float *A_pad, const float *A);

//pad_A的转换，以适用于Openblas
void convertA(float *A_convert, const int outM, const int convAw, const int pad_w, float *A_pad, const int channel);

//kernel转换以适用于openblas
void convertB(const int convAw, const int channel, float *B, float *B_convert);

//Openblas矩阵乘积计算
void Matrixmul_blas(const int convAh, const int convAw, float *A_convert, float *B_convert, float *C, const int channel);

//转换C为常用矩阵排列
void convertC(const int channel, const int convAh, float *C_convert, float *C);

//验证结果是否正确
void evaluate_blas(const int channel, const int Map, const float *A, const int Kernel, const float *B, const int outM, float *C_convert);

#endif // !FUNCTION_H

 

 

 

 

 

函数文件function.cpp

 

//function.cpp
//cblas加速三维矩阵卷积操作
//注意：stride=1
//作者：samylee

#include "function.h"

//A加pad的计算
void comput_Apad(const int pad_w, const int Map, const int channel, float *A_pad, const float *A)
{
	int pad_one_channel = pad_w*pad_w;
	int org_one_channel = Map*Map;
	for (int c = 0; c < channel; c++)
	{
		for (int i = 0; i < pad_w; i++)
		{
			for (int j = 0; j < pad_w; j++)
			{
				int col = c*pad_one_channel + i*pad_w + j;

				if (i == 0 || i == pad_w - 1)
				{
					A_pad[col] = 0;
				}
				else
				{
					if (j == 0 || j == pad_w - 1)
					{
						A_pad[col] = 0;
					}
					else
					{
						A_pad[col] = A[c*org_one_channel + (i - 1)*Map + j - 1];
					}
				}
			}
		}
	}
}

//pad_A的转换，以适用于Openblas
void convertA(float *A_convert, const int outM, const int convAw, const int pad_w, float *A_pad, const int channel)
{
	int pad_one_channel = pad_w*pad_w;
	int seg = channel * convAw;
	for (int c = 0; c < channel; c++)
	{
		for (int i = 0; i < outM; i++)
		{
			for (int j = 0; j < outM; j++)
			{
				int wh = c*convAw + i * outM * seg + j * seg;

				int col1 = c*pad_one_channel + i * pad_w + j;
				A_convert[wh] = A_pad[col1];
				A_convert[wh + 1] = A_pad[col1 + 1];
				A_convert[wh + 2] = A_pad[col1 + 2];

				int col2 = c*pad_one_channel + (i + 1) * pad_w + j;
				A_convert[wh + 3] = A_pad[col2];
				A_convert[wh + 4] = A_pad[col2 + 1];
				A_convert[wh + 5] = A_pad[col2 + 2];

				int col3 = c*pad_one_channel + (i + 2) * pad_w + j;
				A_convert[wh + 6] = A_pad[col3];
				A_convert[wh + 7] = A_pad[col3 + 1];
				A_convert[wh + 8] = A_pad[col3 + 2];
			}
		}
	}
}

//kernel转换以适用于openblas
void convertB(const int convAw, const int channel, float *B, float *B_convert)
{
	int block_A_convert = convAw*channel;
	for (int c = 0; c < channel; c++)
	{
		int block = c*block_A_convert;
		for (int i = 0; i < convAw; i++)
		{
			for (int j = 0; j < channel; j++)
			{
				if (c == j)
				{
					B_convert[block + i*channel + j] = B[c*convAw + i];
				}
				else
				{
					B_convert[block + i*channel + j] = 0;
				}

			}
		}
	}
}

//Openblas矩阵乘积计算
void Matrixmul_blas(const int convAh, const int convAw, float *A_convert, float *B_convert, float *C, const int channel)
{
	const enum CBLAS_ORDER Order = CblasRowMajor;
	const enum CBLAS_TRANSPOSE TransA = CblasNoTrans;
	const enum CBLAS_TRANSPOSE TransB = CblasNoTrans;
	const int M = convAh;//A的行数，C的行数
	const int N = channel;//B的列数，C的列数
	const int K = convAw * channel;//A的列数，B的行数
	const float alpha = 1;
	const float beta = 0;
	const int lda = K;//A的列
	const int ldb = N;//B的列
	const int ldc = N;//C的列

	cblas_sgemm(Order, TransA, TransB, M, N, K, alpha, A_convert, lda, B_convert, ldb, beta, C, ldc);
}

//转换C为常用矩阵排列
void convertC(const int channel, const int convAh, float *C_convert, float *C)
{
	for (int c = 0; c < channel; c++)
	{
		for (int i = 0; i < convAh; i++)
		{
			C_convert[c*convAh + i] = C[i*channel + c];
		}
	}
}

//验证结果是否正确
void evaluate_blas(const int channel, const int Map, const float *A, const int Kernel, const float *B, const int outM, float *C_convert)
{
	cout << "A is:" << endl;
	for (int c = 0; c < channel; c++)
	{
		for (int i = 0; i < Map; i++)
		{
			for (int j = 0; j < Map; j++)
			{
				cout << A[c*Map*Map + i*Map + j] << " ";
			}
			cout << endl;
		}
		cout << endl;
	}

	cout << "B is:" << endl;
	for (int c = 0; c < channel; c++)
	{
		for (int i = 0; i < Kernel; i++)
		{
			for (int j = 0; j < Kernel; j++)
			{
				cout << B[c * Kernel * Kernel + i*Kernel + j] << " ";
			}
			cout << endl;
		}
		cout << endl;
	}

	cout << "C is:" << endl;
	for (int c = 0; c < channel; c++)
	{
		for (int i = 0; i < outM; i++)
		{
			for (int j = 0; j < outM; j++)
			{
				cout << C_convert[c * outM * outM + i*outM + j] << " ";
			}
			cout << endl;
		}
		cout << endl;
	}
}

 

 

 

 

 

主函数文件main.cpp

 

//main.cpp
//cblas加速三维矩阵卷积操作
//注意：stride=1
//作者：samylee

#include "function.h"

int main()
{
	//卷积参数初始化
	const int pad = 1;
	const int stride = 1;

	//定义被卷积三维矩阵
	const int Map = 4;
	const int channel = 3;
	const float A[Map * Map * channel] = {
		1,2,3,4,
		1,2,3,4,
		1,2,3,4,
		1,2,3,4,

		1,2,3,4,
		1,2,3,4,
		1,2,3,4,
		1,2,3,4,

		1,2,3,4,
		1,2,3,4,
		1,2,3,4,
		1,2,3,4 };

	//定义三维卷积核
	const int Kernel = 3;
	float B[Kernel * Kernel * channel] = {
		1,1,1,
		1,1,1,
		1,1,1,

		2,2,2,
		2,2,2,
		2,2,2,

		3,3,3,
		3,3,3,
		3,3,3 };

	//计算卷积输出矩阵宽高
	const int outM = (Map - Kernel + 2 * pad) / stride + 1;

	//计算三维pad_A
	const int pad_w = Map + 2 * pad;
	float A_pad[pad_w*pad_w*channel];
	comput_Apad(pad_w, Map, channel, A_pad, A);

	//定义被卷积矩阵宽高
	const int convAw = Kernel*Kernel;
	const int convAh = outM*outM;

	//转换被卷积矩阵
	float A_convert[convAh*convAw*channel];
	convertA(A_convert, outM, convAw, pad_w, A_pad, channel);

	//转换卷积核以适用于cblas
	float B_convert[channel * Kernel * Kernel * channel];
	convertB(convAw, channel, B, B_convert);
	
	//定义卷积输出矩阵
	float C[convAh * channel];
	//cblas计算输出矩阵
	Matrixmul_blas(convAh, convAw, A_convert, B_convert, C, channel);

	//将输出转换为常用矩阵形式
	float C_convert[outM * outM * channel];
	convertC(channel, convAh, C_convert, C);
	
	//输出验证
	evaluate_blas(channel, Map, A, Kernel, B, outM, C_convert);

	system("pause");
	return EXIT_SUCCESS;
}

 

 

 

 

 

效果如下：

 

 

任何问题请加唯一QQ2258205918（名称samylee）！

或唯一VX：samylee_csdn