CBLAS編譯安裝與使用舉例

在Github上看到有人用BLAS library優化自己的源碼，對此產生了強烈興趣。

准備自己動手實踐一下，網上搜索了一大堆編譯安裝BLAS教程的資料，沒一個靠譜的，編譯過程中遇到一堆的問題。因為自己沒有root權限，所以只能在home目錄中本地編譯使用cblas，然後本地鏈接編譯得到的庫文件到應用程序。

最後自己憑著直覺連蒙帶猜，終於把BLAS與CBLAS裝上，並投入到實例中優化運行應用程序。填補了很多Linux知識。 

首先要解釋一下BLAS，CBLAS與LAPAXK之間的區別與聯系。

• BLAS（Basic Linear Algebra Subprograms）庫，是用Fortran語言實現的向量和矩陣運算庫,是許多數值計算軟件庫的核心, 但也有一些其它的包裝, 如cblas是C語言, 也有C++的包裝, boost/ublas 是C++ template class的實現； 另外還有一些特別的實現, 如intel MKL, AMD core math library blas就是做向量、矩陣的基本運算，如加、減、乘等操作。
• CBLAS是BLAS的C語言接口。
• LAPACK（Linear Algebra PACKage）庫，是用Fortran語言編寫的線性代數計算庫，包含線性方程組求解（AX=b）、矩陣分解、矩陣求逆、求矩陣特征值、奇異值等。該庫用BLAS庫做底層運算，許多高層的數學庫都用BLAS和LAPACK做底層。 

CBLAS只是BLAS的C語言版本，所以CBLAS安裝需要先裝BLAS

安裝步驟

• 確保機器上安裝了gfortran編譯器
• 下載blas, cblas源代碼，都可以在官網http://www.netlib.org/ 上找到
• 解壓後得到兩個文件夾
• 編譯

1. 編譯blas，進入BLAS目錄執行下面的命令

    gfortran -c  -O3    *.f                # 編譯所有的 .f 文件，生成 .o文件 
    ar rv libblas.a      *.o                # 鏈接所有的 .o文件，生成 .a 文件 

2. 編譯cblas，進入CBLAS目錄，首先根據自己的操作系統平台，將某個Makefiel.XXX復制為Makefile.in，XXX表示操作系統。如果是Linux，那麼就將Makefile.LINUX 復制為 Makefile.in。

    cp ../BLAS/libblas.a  testing  # 將上一步編譯成功的 libblas.a 復制到 CBLAS目錄下的testing子目錄 
    make                                            # 編譯所有的目錄 

此時會在CBLAS安裝目錄下的lib目錄中產生一個靜態鏈接庫文件cblas_LINUX.a，這個庫文件和上面得到的libblas.a文件就是我們所需要的。另外還需要的就是CBLAS/include中的cblas.h頭文件。將三個文件全部拷貝到，你需調用的應用程序源碼目錄中。

到此BLAS和CBLAS的安裝任務完成，可以看出，這裡安裝的實際含義是編譯得到兩個庫文件和一個頭文件，再將這三個文件放置到gcc的搜索路徑中去（例如可以在拷貝到/esr/local/lib，或在/usr/local/lib下做一個快捷鏈接，也可直接像我上面那樣復制的）。

cd /usr/local/lib
ln -s  ./CBLAS/lib/cblas_LINUX.a  ./libcblas.a

簡單運用

CBLAS/BLAS分為3個level，level1是用於向量的計算，level2是用於向量和矩陣之間的計算，level3是矩陣之間的計算。比如計算矩陣的乘法就是屬於level3，這裡就用矩陣乘法來學習使用CBLAS。

計算矩陣乘法的函數之一是 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 enum CBLAS_TRANSPOSE TransA和 const enum CBLAS_TRANSPOSE TransB，這兩個參數影響的是op( A )和op( B)，可選參數為CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113，其中TransA = CblasNoTrans, op( A ) = A，TransA = CblasTrans, op( A ) = A'，TransA = CblasConjTrans, op( A ) = A'。 TransB類似。

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的數據

計算兩個簡單矩陣的乘法。
A:
1,2,3
4,5,6
7,8,9
8,7,6
B:
5,4
3,2
1,0

// <strong>因為程序是C++，而CBLAS是C語言寫的，所以在此處用extern關鍵字</strong>
extern"C"
{
    #include"cblas.h"      // <strong>由於cblas.h文件已經拷貝到工作目錄中，只需用雙引號 </strong>
}
#include<iostream>
using namespace std;
int main(void) {
    const enum CBLAS_ORDER Order=CblasRowMajor;
    const enum CBLAS_TRANSPOSE TransA=CblasNoTrans;
    const enum CBLAS_TRANSPOSE TransB=CblasNoTrans;
    const int M=4;//A的行數，C的行數
    const int N=2;//B的列數，C的列數
    const int K=3;//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的列
    const float A[K*M]={1,2,3,4,5,6,7,8,9,8,7,6};
    const float B[K*N]={5,4,3,2,1,0};
    float C[M*N];
 
    cblas_sgemm(Order, TransA, TransB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   
    for(int i=0;i<M;i++)
    {
      for(int j=0;j<N;j++)
      {
          cout<<C[i*N+j]<<"/t";
      }
      cout<<endl;
    }
 
    return EXIT_SUCCESS;
}

在編譯的時候需要帶上cblas_LINUX.a和libblas.a 

比如：g++ main.cpp cblas_LINUX.a libblas.a -o main

如果是寫成Makefile文件，則修改變量OBJS = main.o cblas_LINUX.a libblas.a

當然，這裡假定是這兩個.a文件是放在可以直接訪問的位置，或者寫全路徑也可以。

應用優化

這個有時間再補充，現在這個時間比較尴尬。不過確實對應用起到了優化作用，運行時間從160多s降到了115s，只是簡單將串行代碼用函數clabs_caxpy替代。