梳理caffe代码math_functions(一)

最新推荐文章于 2022-06-05 01:39:09 发布
ZhikangFu
最新推荐文章于 2022-06-05 01:39:09 发布
阅读量685
点赞数
分类专栏： caffe
caffe 专栏收录该内容
30 篇文章 2 订阅
订阅专栏
转载 http://blog.csdn.net/langb2014/article/details/50986678
先从caffe中使用的函数入手看看：
[cpp]view plaincopy 
    
 #include <boost/math/special_functions/next.hpp>  
 #include <boost/random.hpp>  
   
 #include <limits>  
   
 #include "caffe/common.hpp"  
 #include "caffe/util/math_functions.hpp"  
 #include "caffe/util/rng.hpp"  
   
 namespace caffe {  
 //math_function 定义了caffe 中用到的一些矩阵操作和数值计算的一些函数  
 /* 
  *功能： C=alpha*A*B+beta*C 
  *A,B,C 是输入矩阵（一维数组格式） 
  *CblasRowMajor :数据是行主序的（二维数据也是用一维数组储存的） 
  *TransA, TransB：是否要对A和B做转置操作（CblasTrans CblasNoTrans） 
  *M： A、C 的行数 
  *N： B、C 的列数 
  *K： A 的列数， B 的行数 
  *lda ： A的列数（不做转置）行数（做转置） 
  *ldb： B的列数（不做转置）行数（做转置） 
 */  
 template<>  
 void caffe_cpu_gemm<float>(const CBLAS_TRANSPOSE TransA,  
     const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,  
     const float alpha, const float* A, const float* B, const float beta,  
     float* C) {  
   int lda = (TransA == CblasNoTrans) ? K : M;  
   int ldb = (TransB == CblasNoTrans) ? N : K;  
   cblas_sgemm(CblasRowMajor, TransA, TransB, M, N, K, alpha, A, lda, B,  
       ldb, beta, C, N);  
 }  
 /*计算矩阵乘法的函数之一是 cblas_sgemm，使用单精度实数，另外还有对应双精度实数， 
 单精度复数和双精度复数的函数。在此以 cblas_sgemm为例。 
 函数定义为： 
 void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, 
 const enum CBLAS_TRANSPOSE TransB, const int M, const int N, 
 const int K, const float alpha, const float  *A, 
 const int lda, const float  *B, const int ldb, 
 const float beta, float  *C, const int ldc) 
 得到的结果是: 
  C = alpha*op( A )*op( B ) + beta*C 
 const enum CBLAS_ORDER Order，这是指的数据的存储形式，在CBLAS的函数中无论一维还是二维数据 
 都是用一维数组存储，这就要涉及是行主序还是列主序，在C语言中数组是用 行主序，fortran中是列 
 主序。我还是习惯于是用行主序，所以这个参数是用CblasRowMajor，如果是列主序的话就是 CblasColMajor。 
 const int M，矩阵A的行，矩阵C的行 
 const int N，矩阵B的列，矩阵C的列 
 const int K，矩阵A的列，矩阵B的行 
 const float alpha， const float beta，计算公式中的两个参数值，如果只是计算C=A*B，则alpha=1,beta=0 
 const float  *A， const float  *B， const float  *C，矩阵ABC的数据 
 const int lda， const int ldb， const int ldc，在BLAS的文档里，这三个参数分别为ABC的行数， 
 但是实际使用发现，在CBLAS里应该是列数。 
 */  
 template<>  
 void caffe_cpu_gemm<double>(const CBLAS_TRANSPOSE TransA,  
     const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,  
     const double alpha, const double* A, const double* B, const double beta,  
     double* C) {  
   int lda = (TransA == CblasNoTrans) ? K : M;  
   int ldb = (TransB == CblasNoTrans) ? N : K;  
   cblas_dgemm(CblasRowMajor, TransA, TransB, M, N, K, alpha, A, lda, B,  
       ldb, beta, C, N);  
 }  
 /* 
 功能： y=alpha*A*x+beta*y 
 其中X和Y是向量，A 是矩阵 
 M：A 的行数 
 N：A 的列数 
 cblas_sgemv 中的 参数1 表示对X和Y的每个元素都进行操作 
 */  
 template <>  
 void caffe_cpu_gemv<float>(const CBLAS_TRANSPOSE TransA, const int M,  
     const int N, const float alpha, const float* A, const float* x,  
     const float beta, float* y) {  
   cblas_sgemv(CblasRowMajor, TransA, M, N, alpha, A, N, x, 1, beta, y, 1);  
 }  
   
 template <>  
 void caffe_cpu_gemv<double>(const CBLAS_TRANSPOSE TransA, const int M,  
     const int N, const double alpha, const double* A, const double* x,  
     const double beta, double* y) {  
   cblas_dgemv(CblasRowMajor, TransA, M, N, alpha, A, N, x, 1, beta, y, 1);  
 }  
 /* 
 功能： Y=alpha*X+Y 
 N：为X和Y中element的个数 
 */  
 template <>  
 void caffe_axpy<float>(const int N, const float alpha, const float* X,  
     float* Y) { cblas_saxpy(N, alpha, X, 1, Y, 1); }  
   
 template <>  
 void caffe_axpy<double>(const int N, const double alpha, const double* X,  
     double* Y) { cblas_daxpy(N, alpha, X, 1, Y, 1); }  
 /* 
 功能：用常数 alpha 对 Y 进行初始化 
 函数 void *memset(void *buffer, char c, unsigned count) 一般为新申请的内存做初始化， 
 功能是将buffer所指向内存中的每个字节的内容全部设置为c指定的ASCII值, count为块的大小, 
 使用memset函数来初始化数组或者结构体比其他初始化方法更快一点 
 */  
 template <typename Dtype>  
 void caffe_set(const int N, const Dtype alpha, Dtype* Y) {  
   if (alpha == 0) {  
     memset(Y, 0, sizeof(Dtype) * N);  // NOLINT(caffe/alt_fn)  
     return;  
   }  
   for (int i = 0; i < N; ++i) {  
     Y[i] = alpha;  
   }  
 }  
   
 template void caffe_set<int>(const int N, const int alpha, int* Y);  
 template void caffe_set<float>(const int N, const float alpha, float* Y);  
 template void caffe_set<double>(const int N, const double alpha, double* Y);  
 //功能： 给 Y 的每个 element 加上常数 alpha  
 template <>  
 void caffe_add_scalar(const int N, const float alpha, float* Y) {  
   for (int i = 0; i < N; ++i) {  
     Y[i] += alpha;  
   }  
 }  
   
 template <>  
 void caffe_add_scalar(const int N, const double alpha, double* Y) {  
   for (int i = 0; i < N; ++i) {  
     Y[i] += alpha;  
   }  
 }  
 /* 
 函数 void *memcpy(void *dest, void *src, unsigned int count) 把src所指向 
 的内存区域 copy到dest所指向的内存区域, count为块的大小 
 表头文件: #include <string.h> 
 定义函数: void *memcpy(void *dest, const void *src, size_t n) 
 函数说明: memcpy()用来拷贝src所指的内存内容前n个字节到dest所指的内存地址上。与strcpy()不同的是,memcpy()会完整的复制n个字节,不会因为遇到字符串结束'\0'而结束 
 返回值:   返回指向dest的指针 
 */  
 template <typename Dtype>  
 void caffe_copy(const int N, const Dtype* X, Dtype* Y) {  
   if (X != Y) {  
     if (Caffe::mode() == Caffe::GPU) {  
 #ifndef CPU_ONLY  
       // NOLINT_NEXT_LINE(caffe/alt_fn)  
       CUDA_CHECK(cudaMemcpy(Y, X, sizeof(Dtype) * N, cudaMemcpyDefault));  
 #else  
       NO_GPU;  
 #endif  
     } else {  
       memcpy(Y, X, sizeof(Dtype) * N);  // NOLINT(caffe/alt_fn)  
     }  
   }  
 }  
   
 template void caffe_copy<int>(const int N, const int* X, int* Y);  
 template void caffe_copy<unsigned int>(const int N, const unsigned int* X,  
     unsigned int* Y);  
 template void caffe_copy<float>(const int N, const float* X, float* Y);  
 template void caffe_copy<double>(const int N, const double* X, double* Y);  
 /* 
 功能：X = alpha*X 
 N： X中element的个数 
 */  
 template <>  
 void caffe_scal<float>(const int N, const float alpha, float *X) {  
   cblas_sscal(N, alpha, X, 1);  
 }  
   
 template <>  
 void caffe_scal<double>(const int N, const double alpha, double *X) {  
   cblas_dscal(N, alpha, X, 1);  
 }  
 /* 
 功能：Y= alpha*X+beta*Y 
 */  
 template <>  
 void caffe_cpu_axpby<float>(const int N, const float alpha, const float* X,  
                             const float beta, float* Y) {  
   cblas_saxpby(N, alpha, X, 1, beta, Y, 1);  
 }  
   
 template <>  
 void caffe_cpu_axpby<double>(const int N, const double alpha, const double* X,  
                              const double beta, double* Y) {  
   cblas_daxpby(N, alpha, X, 1, beta, Y, 1);  
 }  
 /* 
 功能：这四个函数分别实现element-wise的加减乘除（y[i] = a[i] + - * \ b[i]） 
 */  
 template <>  
 void caffe_add<float>(const int n, const float* a, const float* b,  
     float* y) {  
   vsAdd(n, a, b, y);  
 }  
   
 template <>  
 void caffe_add<double>(const int n, const double* a, const double* b,  
     double* y) {  
   vdAdd(n, a, b, y);  
 }  
   
 template <>  
 void caffe_sub<float>(const int n, const float* a, const float* b,  
     float* y) {  
   vsSub(n, a, b, y);  
 }  
   
 template <>  
 void caffe_sub<double>(const int n, const double* a, const double* b,  
     double* y) {  
   vdSub(n, a, b, y);  
 }  
   
 template <>  
 void caffe_mul<float>(const int n, const float* a, const float* b,  
     float* y) {  
   vsMul(n, a, b, y);  
 }  
   
 template <>  
 void caffe_mul<double>(const int n, const double* a, const double* b,  
     double* y) {  
   vdMul(n, a, b, y);  
 }  
   
 template <>  
 void caffe_div<float>(const int n, const float* a, const float* b,  
     float* y) {  
   vsDiv(n, a, b, y);  
 }  
   
 template <>  
 void caffe_div<double>(const int n, const double* a, const double* b,  
     double* y) {  
   vdDiv(n, a, b, y);  
 }  
 /* 
 功能 : 同样是element-wise操作，分别是y[i] = a[i] ^ b， y[i] = a[i]^2，y[i] = exp(a[i] )，y[i] = |a[i] | 
 */  
 template <>  
 void caffe_powx<float>(const int n, const float* a, const float b,  
     float* y) {  
   vsPowx(n, a, b, y);  
 }  
   
 template <>  
 void caffe_powx<double>(const int n, const double* a, const double b,  
     double* y) {  
   vdPowx(n, a, b, y);  
 }  
   
 template <>  
 void caffe_sqr<float>(const int n, const float* a, float* y) {  
   vsSqr(n, a, y);  
 }  
   
 template <>  
 void caffe_sqr<double>(const int n, const double* a, double* y) {  
   vdSqr(n, a, y);  
 }  
   
 template <>  
 void caffe_exp<float>(const int n, const float* a, float* y) {  
   vsExp(n, a, y);  
 }  
   
 template <>  
 void caffe_exp<double>(const int n, const double* a, double* y) {  
   vdExp(n, a, y);  
 }  
   
 template <>  
 void caffe_log<float>(const int n, const float* a, float* y) {  
   vsLn(n, a, y);  
 }  
   
 template <>  
 void caffe_log<double>(const int n, const double* a, double* y) {  
   vdLn(n, a, y);  
 }  
   
 template <>  
 void caffe_abs<float>(const int n, const float* a, float* y) {  
     vsAbs(n, a, y);  
 }  
   
 template <>  
 void caffe_abs<double>(const int n, const double* a, double* y) {  
     vdAbs(n, a, y);  
 }  
 /* 
 功能：返回一个随机数 
 */  
 unsigned int caffe_rng_rand() {  
   return (*caffe_rng())();  
 }  
 /* 
 功能 ： 返回 b 最大方向上可以表示的最接近的数值。 
 */  
 template <typename Dtype>  
 Dtype caffe_nextafter(const Dtype b) {  
   return boost::math::nextafter<Dtype>(  
       b, std::numeric_limits<Dtype>::max());  
 }  
   
 template  
 float caffe_nextafter(const float b);  
   
 template  
 double caffe_nextafter(const double b);  
   
 template <typename Dtype>  
 void caffe_rng_uniform(const int n, const Dtype a, const Dtype b, Dtype* r) {  
   CHECK_GE(n, 0);  
   CHECK(r);  
   CHECK_LE(a, b);  
   boost::uniform_real<Dtype> random_distribution(a, caffe_nextafter<Dtype>(b));  
   boost::variate_generator<caffe::rng_t*, boost::uniform_real<Dtype> >  
       variate_generator(caffe_rng(), random_distribution);  
   for (int i = 0; i < n; ++i) {  
     r[i] = variate_generator();  
   }  
 }  
   
 template  
 void caffe_rng_uniform<float>(const int n, const float a, const float b,  
                               float* r);  
   
 template  
 void caffe_rng_uniform<double>(const int n, const double a, const double b,  
                                double* r);  
   
 template <typename Dtype>  
 void caffe_rng_gaussian(const int n, const Dtype a,  
                         const Dtype sigma, Dtype* r) {  
   CHECK_GE(n, 0);  
   CHECK(r);  
   CHECK_GT(sigma, 0);  
   boost::normal_distribution<Dtype> random_distribution(a, sigma);  
   boost::variate_generator<caffe::rng_t*, boost::normal_distribution<Dtype> >  
       variate_generator(caffe_rng(), random_distribution);  
   for (int i = 0; i < n; ++i) {  
     r[i] = variate_generator();  
   }  
 }  
   
 template  
 void caffe_rng_gaussian<float>(const int n, const float mu,  
                                const float sigma, float* r);  
   
 template  
 void caffe_rng_gaussian<double>(const int n, const double mu,  
                                 const double sigma, double* r);  
   
 template <typename Dtype>  
 void caffe_rng_bernoulli(const int n, const Dtype p, int* r) {  
   CHECK_GE(n, 0);  
   CHECK(r);  
   CHECK_GE(p, 0);  
   CHECK_LE(p, 1);  
   boost::bernoulli_distribution<Dtype> random_distribution(p);  
   boost::variate_generator<caffe::rng_t*, boost::bernoulli_distribution<Dtype> >  
       variate_generator(caffe_rng(), random_distribution);  
   for (int i = 0; i < n; ++i) {  
     r[i] = variate_generator();  
   }  
 }  
   
 template  
 void caffe_rng_bernoulli<double>(const int n, const double p, int* r);  
   
 template  
 void caffe_rng_bernoulli<float>(const int n, const float p, int* r);  
   
 template <typename Dtype>  
 void caffe_rng_bernoulli(const int n, const Dtype p, unsigned int* r) {  
   CHECK_GE(n, 0);  
   CHECK(r);  
   CHECK_GE(p, 0);  
   CHECK_LE(p, 1);  
   boost::bernoulli_distribution<Dtype> random_distribution(p);  
   boost::variate_generator<caffe::rng_t*, boost::bernoulli_distribution<Dtype> >  
       variate_generator(caffe_rng(), random_distribution);  
   for (int i = 0; i < n; ++i) {  
     r[i] = static_cast<unsigned int>(variate_generator());  
   }  
 }  
   
 template  
 void caffe_rng_bernoulli<double>(const int n, const double p, unsigned int* r);  
   
 template  
 void caffe_rng_bernoulli<float>(const int n, const float p, unsigned int* r);  
   
 template <>  
 float caffe_cpu_strided_dot<float>(const int n, const float* x, const int incx,  
     const float* y, const int incy) {  
   return cblas_sdot(n, x, incx, y, incy);  
 }  
 /* 
 功能： 返回 vector X 和 vector Y 的内积。 
 incx， incy ： 步长，即每隔incx 或 incy 个element 进行操作。 
 */  
 template <>  
 double caffe_cpu_strided_dot<double>(const int n, const double* x,  
     const int incx, const double* y, const int incy) {  
   return cblas_ddot(n, x, incx, y, incy);  
 }  
   
 template <typename Dtype>  
 Dtype caffe_cpu_dot(const int n, const Dtype* x, const Dtype* y) {  
   return caffe_cpu_strided_dot(n, x, 1, y, 1);  
 }  
   
 template  
 float caffe_cpu_dot<float>(const int n, const float* x, const float* y);  
   
 template  
 double caffe_cpu_dot<double>(const int n, const double* x, const double* y);  
 /* 
 功能：返回 x 和 y 之间的海明距离。（两个等长字符串之间的海明距离 
 是两个字符串对应位置的不同字符的个数。） 
 */  
 template <>  
 int caffe_cpu_hamming_distance<float>(const int n, const float* x,  
                                   const float* y) {  
   int dist = 0;  
   for (int i = 0; i < n; ++i) {  
     dist += __builtin_popcount(static_cast<uint32_t>(x[i]) ^  
                                static_cast<uint32_t>(y[i]));  
   }  
   return dist;  
 }  
   
 template <>  
 int caffe_cpu_hamming_distance<double>(const int n, const double* x,  
                                    const double* y) {  
   int dist = 0;  
   for (int i = 0; i < n; ++i) {  
     dist += __builtin_popcountl(static_cast<uint64_t>(x[i]) ^  
                                 static_cast<uint64_t>(y[i]));  
   }  
   return dist;  
 }  
 /* 
 功能：计算 vector x 的所有element的绝对值之和。 
 */  
 template <>  
 float caffe_cpu_asum<float>(const int n, const float* x) {  
   return cblas_sasum(n, x, 1);  
 }  
   
 template <>  
 double caffe_cpu_asum<double>(const int n, const double* x) {  
   return cblas_dasum(n, x, 1);  
 }  
   
 template <>  
 void caffe_cpu_scale<float>(const int n, const float alpha, const float *x,  
                             float* y) {  
   cblas_scopy(n, x, 1, y, 1);  
   cblas_sscal(n, alpha, y, 1);  
 }  
   
 template <>  
 void caffe_cpu_scale<double>(const int n, const double alpha, const double *x,  
                              double* y) {  
   cblas_dcopy(n, x, 1, y, 1);  
   cblas_dscal(n, alpha, y, 1);  
 }  
   
 }  // namespace caffe  
看一下cuda的设计技巧：
在common.hpp中用到一个宏定义CUDA_KERNEL_LOOP
#defineCUDA_KERNEL_LOOP(i,n) \
for(inti = blockIdx.x * blockDim.x + threadIdx.x;i < (n);i +=blockDim.x * gridDim.x)
先看看caffe采取的线程格和线程块的维数设计，还是从common.hpp可以看到
CAFFE_CUDA_NUM_THREADS
CAFFE_GET_BLOCKS(constintN)
明显都是一维的。
整理一下CUDA_KERNEL_LOOP格式看看，
for(inti = blockIdx.x * blockDim.x + threadIdx.x;i< (n);i+= blockDim.x * gridDim.x)
blockDim.x* gridDim.x表示的是该线程格所有线程的数量。
n表示核函数总共要处理的元素个数。
有时候，n会大于blockDim.x* gridDim.x，因此并不能一个线程处理一个元素。
由此通过上面的方法，让一个线程串行（for循环）处理几个元素。
再来看一下这个核函数的实现。
template<typename Dtype>
__global__void mul_kernel(const int n, const Dtype* a,
constDtype* b, Dtype* y)
{
CUDA_KERNEL_LOOP(index,n)
{
y[index]= a[index] * b[index];
}
}
就是算两个向量的点积了。由于向量的维数可能大于该kernel函数线程格的总线程数量。
ZhikangFu
关注
0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
梳理caffe代码math_functions(一)

转载 http://blog.csdn.net/langb2014/article/details/50986678先从caffe中使用的函数入手看看：[cpp] view plain copy #include #include #include #include "caffe/common.
复制链接

扫一扫
专栏目录