主要在R下和Python下使用GPU加速計算，使用成熟的工具來提高自己的工作效率，現在對其中一些方法做一下簡單總結。

R的GPU加速計算包gputools

　　1）gputools，R下的GPU加速計算的函式包，包含常見的運算操作。
　　https://cran.r-project.org/web/packages/gputools/
　　2）iFes，Incremental Feature Selection Algorithm accelerated by GPU。

Python的GPU加速計算包cudamat及gnumpy

Python在GPU下實現的演算法

　　在pypi下，搜尋gpu，肯定有。已知的有cudatree 0.6（隨機森林），dpmix 0.3（高斯混合）。
　　百度/Google下，搜尋Python+GPU+algorithm。

補充程式碼

1. R的gputools測試程式碼

## Objection: Test GPU Compution Function in R ##
## time: 2015.08.27
## author: yjm
## [email protected]~# nvidia-smi  ## check the gpu information
library(gputools)
help(package = 'gputools')
## we can find the Computing Function using GPU ##
## chooseGpu; ## default=0
 

chooseGpu(deviceId=0)
############################################
## cpuMatMult (matrix multiplication); 
matA <- matrix(runif(2000*3000), 2000, 3000)
matB <- matrix(runif(3000*4000), 3000, 4000)
t1 = Sys.time()
y1 = cpuMatMult(matA, matB)
t2 = Sys.time()
t2-t1
#### without GPU also
y2 = matA %*% matB
t3 = Sys.time()
t3-t2
y3 = crossprod(t(matA), matB)
t4 = Sys.time()
t4-t3
#############################################
 

## getGpuId; 
?getGpuId
getGpuId()
## default is device 0
## gpuCor (compute correlation coefficient--"pearson & kendall"); Matrix coef
numAvars <- 5
numBvars <- 10
numSamples <- 30
A <- matrix(runif(numAvars*numSamples), numSamples, numAvars)
B <- matrix(runif(numBvars*numSamples), numSamples, numBvars)
gpuCor(A, B, method="pearson")
gpuCor(A, B, method="kendall")
A[3,2] <- NA
gpuCor(A, B, use="pairwise.complete.obs", method="pearson")
#### without GPU
cor.test(A[,1], B[,1], method='pearson')
#########################################################
## gpuCrossprod (cross-product)
matA <- matrix(runif(3000*2000), 3000, 2000)
matB <- matrix(runif(3000*4000), 3000, 4000)
t1 = Sys.time()
y1 = gpuCrossprod(matA, matB)
t2 = Sys.time()
t2-t1
#### without GPU
y2 = t(matA) %*% matB
t3 = Sys.time()
t3-t2
y3 = crossprod(matA, matB)
t4 = Sys.time()
t4-t3
#############################################################
## gpuDist(matrix, method) (comput distance between vec) each row is a vector
numVectors <- 500
dimension <- 1000
Vectors <- matrix(runif(numVectors*dimension), numVectors, dimension)
t1 = Sys.time()
y1 = gpuDist(Vectors, "euclidean")
t2 = Sys.time()
t2-t1
#### without GPU
y2 = dist(Vectors, "euclidean")
t3 = Sys.time()
t3-t2
#gpuDist(Vectors, "maximum")
#gpuDist(Vectors, "manhattan")
#gpuDist(Vectors, "minkowski", 4)
###################################################
## gpuDistClust
## gpuGlm (glm with gpu)
## gpuGranger (granger causality tests)
## gpuHclust (Hierarchical Clustering)
## gpuLm (lm with gpu)
## gpuLm.defaultTol
## gpuLm.fit 
## gpuLsfit (least sqares fit)
## gpuMatMult (matrix multiplication with GPU)
matA <- matrix(runif(2000*3000), 2000, 3000)
matB <- matrix(runif(3000*4000), 3000, 4000)
t1 = Sys.time()
y1 = gpuMatMult(matA, matB)
t2 = Sys.time()
t2-t1
#### without GPU
y2 = matA %*% matB
t3 = Sys.time()
t3-t2
y3 = crossprod(t(matA), matB)
t4 = Sys.time()
t4-t3
###############################################
## gpuMi () # mutual information computing # each col represent a random variable
x <- matrix(runif(60), 20, 3)
y <- matrix(runif(60), 20, 3)
# do something interesting
y[,2] <- 3.0 * (x[,1] + x[,3])
z <- gpuMi(x, y, bins = 10, splineOrder = 3)
print(z)
## gpuQr (QR decomposition)
## gpuSolve 
## gpuTcrossprod (matrix Transposed Cross-product with GPU)
matA <- matrix(runif(2000*3000), 2000, 3000)
matB <- matrix(runif(4000*3000), 4000, 3000)
t1 = Sys.time()
y1 = gpuTcrossprod(matA, matB)
t2 = Sys.time()
t2 - t1
#### without GPU
y2 = matA %*% t(matB)
t3 = Sys.time()
t3 - t2
y3 = tcrossprod(matA, matB)
t4 = Sys.time()
t4 - t3
## gpuTtest (T-test with gpu)
quit()

## conclusion: 
## when nrow and ncol > 1000, GPU compution is larger quick than cpu.(more 100 times)

2. Python的cudamat包測試程式碼

# coding: utf-8

# In[179]:
## here want to test cudmat package
## time: 2015.09.07
## author: yjm
## cudamat.py檔案裡面有作者的註釋，可以檢視 

#import pycuda as cuda
#import pycuda.autoinit
#from pycuda.compiler import SourceModule
import cudamat as cm
import nose
import numpy as np

# 1)
cm.cublas_init()
#print(cm.CUDAMatrix.ones.shape)
# 2)
#cm.cublas_shutdown() ## 這裡面對cm.CUDAMatrix.ones 置零

# In[180]:

## 重置形狀 ##
m = 256
n = 128
cm1 = np.array(np.random.rand(n, m)*10, dtype = np.float32, order = 'F')
cm2 = np.array(np.random.rand(m, 1)*10, dtype = np.float32, order = 'F')
gm1 = cm.CUDAMatrix(cm1)
gm2 = cm.CUDAMatrix(cm2)
# maybe a error here in CUDAMatrix 
# soloution ：Try to fix the value "gpu_memory" of your .pbtxt file to "2G" or "2.5G"
print('----顯示CPU下各個變數的大小----')
print(cm1.shape)
print(cm2.shape)
print('----顯示GPU下各個變數的大小----')
print(gm1.shape)
print(gm2.shape)
print('test reshape in gpu')
gm1.reshape((m, n))
print(gm1.shape)
print('test assign')
#gm2.assign(gm1)
gm1.reshape((n,m))
print(gm1.shape)

# In[81]:

## GPU點乘 ## 以及GPU變數的轉置 ##
m = 256
n = 128
cm1 = np.array(np.random.rand(n, m)*10, dtype = np.float32, order = 'F')
cm2 = np.array(np.random.rand(m, n)*10, dtype = np.float32, order = 'F')
gm1 = cm.CUDAMatrix(cm1)
gm2 = cm.CUDAMatrix(cm2)
print(gm1.shape)
print(gm2.shape)
gm = cm.dot(gm1, gm2) ## here is dot on GPU ##
print(gm.shape)
gm = cm.dot(gm2.T, gm1.T) ## here is transpose on GPU ##
print(gm.shape)

# In[63]:

## assign ## ??? what is assign ??? 設定值
cm2 = np.array(np.random.rand(n, m)*10, dtype=np.float32, order='F')
gm2 = cm.CUDAMatrix(cm2)
print('----cm2-----')
print(cm2[1:5, 1:5])
print('----cm1-----')
print(cm1[1:5, 1:5])
gm1.assign(gm2)
gm1.copy_to_host()
print('----after gm1.assign(gm2)----')
gm1.copy_to_host()
print(gm1.numpy_array[1:5, 1:5])

# In[64]:

## assign ##
a = np.array(np.random.rand(m, n)*10, dtype=np.float32, order = 'F')
m1 = cm.CUDAMatrix(a)
m1.assign(np.pi)
m1.copy_to_host()
print(m1.numpy_array[1:5, 1:5])

# In[74]:

## 獲得行的切片 ## 直接在GPU上進行的切片操作 ##
m = 256
n = 128
start = 11
end = 15

a = np.array(np.random.rand(m, n)*10, dtype=np.float32, order='F')
b = np.array(np.random.rand(end-start, n)*10, dtype=np.float32, order='F')
c = np.array(a[start:end,:], order='F')

m1 = cm.CUDAMatrix(a)
m2 = cm.CUDAMatrix(b)
print(m1.shape)
print(m2.shape)

m1.get_row_slice(start, end, target = m2)
m3 = m1.get_row_slice(start, end)
m1.copy_to_host()
m2.copy_to_host()
m3.copy_to_host()
print('--after m1.get_row_slice(start, end, target = m2)---')
print(m1.shape)
print(m2.shape)
print(m3.shape)
print(m1.numpy_array[start:end, 1:5])
print(m2.numpy_array[:, 1:5])
print(m3.numpy_array[:, 1:5])

# In[95]:

## 將列向量 對應位置值 加到矩陣的每一列 ## add_col_vec ## 
## sample as following:
## a=[1, 2, 3]; b =[1,1,1; 2,2,2;3,3,3] ## then b.add_col_vec(a) = [2,3,4; 3,4,5; 4,5,6]
## 將行向量 對應位置值 加到矩陣的每一行 ## add_row_vec ##
m = 256
n = 128
a = np.array(np.random.rand(m, n)*10, dtype=np.float32, order='F')
b = np.array(np.random.rand(m, 1)*10, dtype=np.float32, order='F')
c = a + b
m1 = cm.CUDAMatrix(a)
m2 = cm.CUDAMatrix(b)
print('---a-m1[1:5, 1:5]---')
print(a[1:5, 1:5])
print('---b-m2[1:5, :]---')
print(b[1:5, :])
print('---c-m1+m2---')
print(c[1:5, 1:5])

print('----after m1.add_col_ve(m2)----')
#m1.add_col_vec(m2, target = m3) 
## 搞不懂這地方為什麼要加target ？？？## 這一句也沒有什麼作用啊！！！
## soga ！ target 是為了儲存結果 ##
m1.add_col_vec(m2)  ## 將m2列向量直接疊加到矩陣m1的每一列 ##
m1.copy_to_host()
print(m1.numpy_array[1:5, 1:5])

# In[103]:

## GPU矩陣所有列資料加上 向量*一個數值 ## 看不懂這個函式有神馬意思 ？？？##
m = 256
n = 128
a = np.array(np.random.rand(m, n)*10, dtype=np.float32, order='F')
b = np.array(np.random.rand(m, 1)*10, dtype=np.float32, order='F')
m1 = cm.CUDAMatrix(a)
m2 = cm.CUDAMatrix(b)
print('---a--m1[1:5, 1:5]----')
print(a[1:5, 1:5])
print('---b--m2[1:5, :]----')
print(b[1:5, :])
m1.add_col_mult(m2, np.pi)
m1.copy_to_host()
print(m1.numpy_array[1:5, 1:5])

# In[104]:

## mult_by_row; mult_by_row; div_by_col; div_by_row; ## who knows what's this! ##

# In[175]:

# 求和 ## sum
m = 256
n = 128
a = np.array(np.random.rand(m, n)*10, dtype=np.float32, order='F')
rowSumRes = np.array(np.random.rand(m, 1)*10, dtype=np.float32, order='F') ## 儲存行sum結果##
colSumRes = np.array(np.random.rand(1, n)*10, dtype=np.float32, order='F') ##
m1 = cm.CUDAMatrix(a)
growSumRes = cm.CUDAMatrix(rowSumRes)
gcolSumRes = cm.CUDAMatrix(colSumRes)
print('---a--m1[1:5, 1:5]----')
print(a[1:5, 1:5])
mult = 1 ## 倍數控制 ##
m1.sum(axis = 1, target = growSumRes, mult = mult)
m1.sum(axis = 0, target = gcolSumRes, mult = mult)
growSumRes.copy_to_host()
print(growSumRes.numpy_array[1:5,:])
print(gcolSumRes.numpy_array[:,1:5])

# In[177]:

## 求均值 ## mean
m = 256
n = 128
a = np.array(np.random.rand(m, n)*10, dtype=np.float32, order='F')
rowMeaRes = np.array(np.random.rand(m, 1)*10, dtype=np.float32, order='F') ## 儲存行sum結果##
colMeaRes = np.array(np.random.rand(1, n)*10, dtype=np.float32, order='F') ##
m1 = cm.CUDAMatrix(a)
growMeaRes = cm.CUDAMatrix(rowSumRes)
gcolMeaRes = cm.CUDAMatrix(colSumRes)
print('---a--m1[1:5, 1:5]----')
print(a[1:5, 1:5])
mult = 1 ## 倍數控制 ##
m1.mean(axis = 1, target = growMeaRes)
m1.mean(axis = 0, target = gcolMeaRes)
growMeaRes.copy_to_host()
print(growMeaRes.numpy_array[1:5, :])
print(gcolMeaRes.numpy_array[:, 1:5])

# In[181]:

## 最大值最小值 ## max/min
## max(axis, target = None)
## min(axis, target = None)
## 跟求和sum 和 求均值mean 一樣用法。
## sign 函式是個什麼東東 ？？？##

# In[185]:

## sigmoid 函式 ## apply_sigmoid(target=None)
m = 256
n = 128
a = np.array(np.random.randn(m, n)*10, dtype=np.float32, order='F')
b = np.array(np.random.randn(m, n)*10, dtype=np.float32, order='F')

c = 1. / (1. + np.exp(-a))

m1 = cm.CUDAMatrix(a)
m2 = cm.CUDAMatrix(b)
m1.apply_sigmoid(target = m2)
m1.apply_sigmoid()

m1.copy_to_host()
m2.copy_to_host()
print(m1.numpy_array[1:5, 1:5])
print(m2.numpy_array[1:5, 1:5])
print(c[1:5, 1:5])

## 雙曲正切函式 ## tanh ##
## gm.apply_tanh(target = gm2)
## gm.apply_tanh()
## 軟閾值函式 ## soft_threshold ##

## log 函式 ## exp 函式 ## sqrt函式 ## 指數函式 ## where 函式 ##
## log(mat, target = None)
## exp(mat, target = None)
## pow(mat, p, target = None)
## where(condition_mat, if_mat, else_mat, target=None)