梳理caffe代碼math_functions(一)

標籤：

先從caffe中使用的函數入手看看：

#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*Cconst 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=0const 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+YN：為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 加上常數 alphatemplate <>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*XN： 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());}templatefloat caffe_nextafter(const float b);templatedouble 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();  }}templatevoid caffe_rng_uniform<float>(const int n, const float a, const float b,                              float* r);templatevoid 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();  }}templatevoid caffe_rng_gaussian<float>(const int n, const float mu,                               const float sigma, float* r);templatevoid 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();  }}templatevoid caffe_rng_bernoulli<double>(const int n, const double p, int* r);templatevoid 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());  }}templatevoid caffe_rng_bernoulli<double>(const int n, const double p, unsigned int* r);templatevoid 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);}templatefloat caffe_cpu_dot<float>(const int n, const float* x, const float* y);templatedouble 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函數線程格的匯流排程數量。

梳理caffe代碼math_functions(一)