[BLAS] Matrix Multiplication Library

BLAS는 Basic Linear Algebra SubPrograms의 약자로써 기초 선형대수(행렬곱) 에 대한 표준(거의 표준) 규격 안이다.
이러한 BLAS의 종류로는 기업용으로는 AMD ACML , Intel MKL , IBM ESSL , Nvidia CUBLAS , Apple Accelearate 등이 있고 오픈 소스로는 Netlib , ATLAS , GotoBLAS , OpenBLAS 가 있다.

Intel계열 CPU를 사용한다면 MKL(2020기준 oneMKL)이 제일 성능이 좋다. 범용적으로 사용하려면 OpenBLAS를 사용하면 된다.

OpenBLAS라면 기본적으로 Fortran77 라이브러리 이며 C언어 버전을 제공한다.
C++이라면 굳이 이딴 저수준 라이브러리를 쓸 필요없이 OpenCV나 AMP등을 사용하면 된다.

BLAS는 level 1,2,3의 표준 규격이 존재하는데 level1은 vector vector 연산, level2는 vector matrix 연산 level3는 matrix matrix 연산이다.

또한, 접두사로 아래와 같은 것들이 붙는다.

s : Single Precision
d : Double Precision
c : Complex
z : Half Complex(16bit)

level3를 주로 사용하므로 level3에 있는 gemm 함수에 대해 설명한다. level3 에 있는 나머지 연산들은 사실 쓸데가 거의 없거나, 최적화를 위해 있는것이므로 gemm만 알면 된다. 어차피 행렬곱말고 다른건 자주 쓰이지 않는다.

float을 사용할 경우 sgemm(Single Precision) double을 사용할 경우 dgemm(Double Precision) , Complex 타입을 사용할 경우 cgemm(Complex) , 16bit-Complex를 사용할 경우 zgemm(16bit-Complex) 을 사용하면 된다.
OpenBLAS에서의 정확한 함수명은 cblas_sgemm, cblas_dgemm 이다. 주로 "젬" 이라고 한다.

level3

gemm

gemm의 뜻은 GEneral Matrix-Matrix Multiplication 의 약자이다.
왜 행렬곱 따위에 이런 규격이 있어야 하나 싶겠지만 이 함수는 단순 AB=C 를 계산해 주는 함수는 아니고, αAB + βC = C 의 수식을 계산한다.

여기서 βC 에서 C는 결과의 C가 아닌 입력에서의 C이다. C는 즉 입력과 동시에 출력이 되는 것이다.

α,β는 스칼라값이며 A,B,C는 각각 행렬이다. 잘은 모르겠지만 이런 형태로 여러가지를 표현할 수 있어서 그런것 같다. 예를들어 α를 1.0 으로 주고 β를 0.0 으로 준다면 AB=C를 계산할 수 있다.

BLAS에 대한 레퍼런스는 어디에서 보나 비슷할테지만 여기가 제일 잘 정리되어 있다.

아무튼 gemm 시리즈의 원형은 아래와 같다.

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언어 저수준 라이브러리라 그런지 파라매터가 많다.

행렬 곱셈의 정의는 아래와 같다.

A, B를 각각 m × k, k × n 행렬이라고 하자. A와 B의 곱 AB는 다음과 같은 항을 갖는 m × n 행렬로 정의된다.

Order : CblasRowMajor 또는 CblasColMajor 인데 C언어는 RowMajor가 기본이다.
TransA : A의 전치행렬을 쓸것인지 전치가 아니라면 CblasNoTrans 전치라면 CblasTrans.
TransB : 위와 동일

눈치 챘겠지만 4~6번째 파라매터는 위의 정의에서의 M,N,K를 그대로 쓴다.
alpha,와 beta는 위에서 설명한 그 α,β 이고
A,B는 입력행렬이고 C는 출력행렬인데 모두 1차원(C언어) 형태의 메모리블록이여야 한다.
이제 좀 짜증나는 부분이 lda,ldb,ldc이다.
각각 정의를 하자면
lda : NoTrans일 경우 K , Trans일 경우 M
ldb : NoTrans일 경우 N , Trans일 경우 K
ldc : C행렬의 첫번째 차원의 크기, A,B가 모두 NoTrans일 경우 M. Trans일 경우 계산해서 넣어라!

당연하지만 이 blas를 이용한 행렬곱이 기존의 행렬곱 보다 매우 빠르다. blas는 하드웨어에 최적화를 하기때문에 복잡한 스트라센 같은 알고리즘을 사용하지 않는다.
일반적으로 ATLAS < OpenBLAS < MKL 순으로 속도가 빠르며 일반적으로 구현하는 행렬곱 과는 비교할수 없는 속도 차이를 보여준다.
다만, 행렬의 수가 매우 적을때(100 이하) 에서는 BLAS보다 간단한 행렬곱이 빠르다.

이러한 gemm 을 간단하게 구현해본다면 아래와 같다.
Rowmajor만 동작하게 구현하였다.

void my_sgemm(int major, int transA, int transB, int M, int N, int K, float alpha, float* A, int lda, float* B, int ldb, float beta, float* C, int ldc) {
	//aAB+bC=C
	for (int m = 0; m < M; m++) {
		for (int n = 0; n < N; n++) {
			float c = C[m*ldc + n];
			C[m*ldc + n] = 0.0F;
			for (int k = 0; k < K; k++) {
				float a, b;
				a = transA == 0 ? A[k*lda + m] : A[m*lda + k];
				b = transB == 0 ? B[n*ldb + k] : B[k*ldb + n];
				C[m*ldc + n] += a*b;
			}
			C[m*ldc + n] = alpha*C[m*ldc + n] + beta*c;
		}
	}
}

Rowmajor 행렬곱을 Colmajor 행렬곱으로 바꾸기.

위의 A(3x2) , B(2x3) , C(3x3) 행렬을 보자. 이 행렬은 Rowmajor 행렬곱이고, 메모리는 실제로 아래와 같다.
A : 1 2 3 4 5 6
B : 1 2 3 4 5 6
C : 9 12 15 19 26 33 29 40 51

만일 Blas가 Rowmajor 곱을 지원하지 않는다면 어떨까(cublas).

이 행렬을 Colmajor라고 생각하고 곱셈을 해야 한다.

Rowmajor와 Colmajor를 변환시키는 방법은 단순한 전치행렬이므로 아래의 공식을 이용하면 된다.

간단히 말하면 AxB 가 아닌 BxA를 하면 된다. Colmajor 상태이므로 BxA 는 실제로는 Bt x At 일 테고 이의 결과는 (AB)t 일 것이므로 이는 사실 Rowmajor상태에서 올바른 결과값이 된다.

sbmv

Scales a Summetric band matrix, then multiplies by a vector, then adds a vector.

sbmv는 y=αAx + βy (Matrix A, Vector x,y)이고 앞으로 BLAS에서 A,B,C가 나오면 행렬, x,y가 나오면 벡터, alpha, beta 가 나오면 스칼라라고 보면 된다.

void cblas_ssbmv(const enum CBLAS_ORDER __Order, const enum CBLAS_UPLO __Uplo, const int __N, const int __K, const float __alpha, const float *__A, const int __lda, const float *__X, const int __incX, const float __beta, float *__Y, const int __incY);

sbmv는 주로 element-wise 곱셉 연산을 할 때 사용된다. 원리는 아래 그림 참조.

이 함수는 원래 벡터(A)를 대각행렬로 바꾸어 계산한다.

UPLO는 "L" 또는 "U"만 가능한데, L은 대각행렬, U는 반대각행렬이다. N과 K는 gemm과 같고, M은 무조건 1이니 입력하지 않는다.

level2

gemv

General Matrix-Vector Multiplication 의 약자이다.

위 그림을 보면 쉽게 이해가 될 것이다.
gemv는 스칼라값 α,β 와 벡터 X,Y 행렬 A가 있을때 벡터 Y는 아래의 수식으로 계산한다.
αAX + βY = Y

마찬가지로 BLAS에서는 모든 반환값이 입력값으로도 작용하니, 반환값의 초기화가 매우 중요하다.

아래는 간단한 예시 코드이다.
gemv example

float X[] = {
		1.0F
		,2.0F
		,3.0F
		,4.0F
	};
	float A[] = {
		1.0F,1.0F,1.0F,1.0F
		,2.0F,2.0F,2.0F,2.0F
		,3.0F,3.0F,3.0F,3.0F
	};
	int n = 4;
	int m = 3;
	float Y[3] = { 0 };
	cblas_sgemv(CblasRowMajor, CblasNoTrans, m, n, 1.0F, A, n, X, 1, 0.0F, Y, 1);
	for (int i = 0; i < 3; i++) {
		printf("[%f]\n", Y[i]);
	}

level1

dot

그냥 dot product 이다.
dot example

float X[] = {1.0F,2.0F,3.0F,4.0F};
float Y[] = { 1.0F,2.0F,3.0F,4.0F };
int n = 4;
float dot=cblas_sdot(n, X, 1, Y, 1);
printf("%f\n", dot);

axpy

Ax + y 의 약자이며 문자 그대로
스칼라값 A와 벡터 X,Y가 있을때 Y=AX+Y 를 계산한다.

float X[] = { 1.0F,2.0F,3.0F,4.0F };
float Y[4] = { 0 };
int n = 4;
cblas_saxpy(n, 2.0F, X, 1, Y, 1);
for (int i = 0; i < 4; i++) {
    printf("[%f]", Y[i]);
}

Build

https://github.com/springkim/WSpring 에서 OpenBLAS 라이브러리 빌드 배치 파일을 제공한다.

Reference

• Apple의 gemm 문서
• Wikipedia의 행렬곱 정의
• 행렬계산기
• openblas에 대한 간단한 설명
• level3
• level2
• level1
• gemm trick