[OpenBLAS] Matrix Multiplication Library

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

이런 라이브러리를 선택할때 무었을 선택할지 고민된다면 그냥 앞에 Open 붙은걸 고르면 된다.
OpenCV, OpenSSL, OpenGL 등등 Open이 붙어있으면 내가 알기론 대부분 사용하기 쉽고 라이선스 문제에서도 자유롭다. 또한 Open이 붙어있을경우 8할 이상은 Cross-Platform이다.

기본적으로 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를 계산할 수 있다.

OpenBLAS는 오픈소스라고 대놓고 문서를 다른데서 보라고 하는데, 그곳의 문서가 썩 맘에 들지 않는다.

아무튼 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상태에서 올바른 결과값이 된다.

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

kimbom

Read more posts by this author.

Seoul