MATLAB、cuBLAS、numpy矩阵乘法性能对比分析

import numpy as np
import time

s=10000

# float类型
a=np.ones((s,s),dtype=np.float32)
b=np.ones((s,s),dtype=np.float32)

t=time.time()
c=np.dot(a,b)
print(time.time()-t)

# double类型
a=np.ones((s,s),dtype=np.float64)
b=np.ones((s,s),dtype=np.float64)

t=time.time()
c=np.dot(a,b)
print(time.time()-t)

#include <cublas_v2.h>

template <typename T = float>
void testMatrixTime2(int M = 10000, int N = 10000, int K = 10000) {
	T* A, * B, * C;
	T* dev_A, * dev_B, * dev_C;
	T alpha = 1, beta = 0;
	// 初始化cuda句柄
	cublasHandle_t handle;
	cublasCreate(&handle);

	A = (T*)malloc(M * K * sizeof(T));
	B = (T*)malloc(K * N * sizeof(T));
	C = (T*)malloc(M * N * sizeof(T));

	cudaMalloc(&dev_A, M * K * sizeof(T));
	cudaMalloc(&dev_B, K * N * sizeof(T));
	cudaMalloc(&dev_C, M * N * sizeof(T));

	for (int i = 0; i < M * K; ++i) {
		A[i] = 1;
	}

	for (int i = 0; i < K * N; ++i) {
		B[i] = 1;
	}
	clock_t t = clock();

	cudaMemcpy(dev_A, A, M * K * sizeof(T), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_B, B, K * N * sizeof(T), cudaMemcpyHostToDevice);

	// double计算时改为cublasDgemm
	cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
		M, N, K,
		&alpha,
		dev_A, M,
		dev_B, K,
		&beta,
		dev_C, M
	);

	 cudaMemcpy(C, dev_C, M * N * sizeof(T), cudaMemcpyDeviceToHost);
	 cudaThreadSynchronize();
	 cout << "time" << (double)(clock() - t) / CLOCKS_PER_SEC << endl;
}

s=10000;

% 双精度
a=ones(s,s);
b=ones(s,s);
tic;
c=a*b;
toc;

% 单精度
a=single(a);
b=single(b);
tic;
c=a*b;
toc;