CBLAS 使ってみた

CBLAS とは

• BLAS を C 言語から呼び出すためのライブラリ
• ベクトルと行列に関する基本的な演算ができる
• 行優先か列優先か選択できる

ビルド

• LAPACK/BLAS の設定メモ （mac）
• LAPACK/BLAS の設定メモ （Ubuntu）

サンプルコード

内積

$x^T y$

dot.c

#include <stdio.h>
#include <math.h>
#include <cblas.h>

int main(void){
    double X[] = {1.0,2.0,3.0,4.0,5.0};
    double Y[] = {6.0,7.0,8.0,9.0,0.0};
    int N = 5;

    printf("inner product: %f\n", cblas_ddot(N,X,1,Y,1));

    return 0;
}

出力結果

$ gcc -o dot dot.c -lcblas -lblas -lm
$ ./dot
inner product: 80.000000

ノルム

$l_1$ノルム： $||x||_1$

$l_2$ノルム： $||x||_2$

nrm.c

#include <stdio.h>
#include <math.h>
#include <cblas.h>

int main(void){
    double x[] = {1.0,2.0,3.0,4.0,5.0};
    int N = 5;

    printf("1-norm,2-norm: %f,%f\n",
    cblas_dasum(N,x,1),cblas_dnrm2(N,x,1));

    return 0;
}

出力結果

$ gcc -o nrm nrm.c -lcblas -lblas -lm
$ ./nrm
1-norm,2-norm: 15.000000,7.416198

ベクトルと行列の演算

$y := \alpha Ax + \beta y$

gemv.c

#include <stdio.h>
#include <math.h>
#include <cblas.h>

void print_vec(double *v,int v_dim){
    for (int i = 0; i < v_dim; ++i){
        if (i == v_dim-1)
            printf("%.3f\n",v[i]);
        else
            printf("%.3f ",v[i]);
    }
}

int main(void){
    double A[] = {1.0,2.0,3.0,4.0,  5.0,6.0,7.0,8.0,  9.0,10.0,11.0,12.0};
    double x[] = {10.0,20.0,30.0,40.0};
    double y[] = {50.0,60.0,70.0};
    int M = 3, N = 4, lda = N;
    double alpha = 2.0, beta = 3.0;

    cblas_dgemv(CblasRowMajor, CblasNoTrans, M, N, alpha, A, N, 
    x, 1, beta, y, 1);

    print_vec(y, M);

    return 0;
}

出力結果

$ gcc -o gemv gemv.c -lcblas -lblas -lm
$ ./gemv
750.000 1580.000 2410.000

２つの行列の演算

$C := \alpha AB + \beta C$

gemm.c

#include <stdio.h>
#include <math.h>
#include <cblas.h>

void print_matrix(double *data,int row,int col){
    int i,j;

    for (i = 0; i < row; ++i){
        for (j = 0; j < col; ++j){
            printf("%.3f ", data[i*col + j]);
        }
        printf("\n");
    }
}

int main(void){
    double A[6] = {1.0,2.0,  3.0,4.0,  5.0,6.0};
    double B[8] = {1.0,2.0,3.0,4.0,  5.0,6.0,7.0,8.0};
    double C[12];
    int M = 3, N = 4, K = 2;
    double alpha = 1.0, beta = 0.0;

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, M, N, K,
    alpha,A,K,B,N,beta,C,N);

    print_matrix(C,M,N);

    return 0;
}

出力結果

$ gcc -o gemm gemm.c -lcblas -lblas -lm
$ ./gemm
11.000 14.000 17.000 20.000
23.000 30.000 37.000 44.000
35.000 46.000 57.000 68.000