MKL庫中基本線性代數子程序,BLAS(Basic Linear Algebra Subprograms)庫,是一個API標淮,用以規範發佈基礎線性代數操作的數值庫(如向量或矩陣乘法)。其中CBLAS是BLAS的C語言接口。
庫中前綴用來區分所支持處理的數據類型。
|
前綴
|
描述
|
函數名系列
|
描述
|
|
s-
|
實數、單精度
|
ge...
|
一般矩陣
|
|
c-
|
複數、單精度
|
sy...
|
對稱矩陣
|
|
d-
|
實數、雙精度
|
he...
|
Hermitian矩陣
|
|
z-
|
複數、雙精度
|
tr...
|
三角矩陣
|
基本矩陣、向量操作
|
函數(採用常規的前綴為d的接口)
|
描述
|
|
cblas_dasum
|
向量元素值模的總和
|
|
cblas_daxpy
|
縮放向量
|
|
cblas_dcopy
|
複製向量
|
|
cblas_ddot
|
向量點積
|
|
cblas_dswap
|
交換兩向量
|
|
cblas_dgemv
|
常規矩陣×向量
|
重點介紹矩陣的乘法運算。
此示例是利用Intel 的MKL庫函數完成計算矩陣(乘法)運算,計算式為:
其通過調整\(A、B、C\)矩陣及其係數,同樣可以完成矩陣的加減;如若只需矩陣\(A\)與\(B\)的乘積,設置\(\alpha=1,\beta=0\)即可。
其中\(A\)為\(m\times k\)維矩陣,\(B\)為\(k\times n\)維矩陣,\(C\)為\(m\times n\)維矩陣。
使用的函數為cblas_?gemm(gemm表示GEneric Matrix Multiplication),完成一般的矩陣乘法。
根據輸入/輸出數據的類型可以分為cblas_dgemm,cblas_sgemm,cblas_cgemm,cblas_zgemm,具體類型參見上文,不再贅述,以下以cblas_dgemm為例介紹其用法。
1 cblas_dgemm參數詳解
fun cblas_dgemm(Layout, //指定行優先(CblasRowMajor,C)或列優先(CblasColMajor,Fortran)數據排序
TransA, //指定是否轉置矩陣A,可以為CblasNoTrans或CblasTrans
TransB, //指定是否轉置矩陣B,可以為CblasNoTrans或CblasTrans
M, //矩陣A和C的行數
N, //矩陣B和C的列數
K, //矩陣A的列,B的行
alpha, //矩陣A和B乘積的比例因子
A, //A矩陣
lda, //矩陣A的第一維的大小
B, //B矩陣
ldb, //矩陣B的第一維的大小
beta, //矩陣C的比例因子
C, //(輸入/輸出) 矩陣C的地址
ldc //矩陣C的第一維的大小
)
2 定義待處理矩陣
#include <stdio.h>
#include <stdlib.h>
#include "mkl.h" // 調用mkl頭文件
#define min(x,y) (((x) < (y)) ? (x) : (y))
double* A, * B, * C; //聲明三個矩陣變量,並分配內存
int m, n, k, i, j; //聲明矩陣的維度,其中
double alpha, beta;
m = 2000, k = 200, n = 1000;
alpha = 1.0; beta = 0.0;
A = (double*)mkl_malloc(m * k * sizeof(double), 64); //按照矩陣維度分配內存
B = (double*)mkl_malloc(k * n * sizeof(double), 64); //mkl_malloc用法與malloc相似,64表示64位
C = (double*)mkl_malloc(m * n * sizeof(double), 64);
if (A == NULL || B == NULL || C == NULL) { //判空
mkl_free(A);
mkl_free(B);
mkl_free(C);
return 1;
}
for (i = 0; i < (m * k); i++) { //賦值
A[i] = (double)(i + 1);
}
for (i = 0; i < (k * n); i++) {
B[i] = (double)(-i - 1);
}
for (i = 0; i < (m * n); i++) {
C[i] = 0.0;
}
其中\(A\)和\(B\)矩陣設置為:
\(C\)矩陣為全0。
3 執行矩陣乘法
回到例子中,對照上面的參數,將C矩陣用A與B的矩陣乘法表示:
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, k, B, n, beta, C, n);
執行後的得到結果如下:
完整代碼
#include <stdio.h>
#include <stdlib.h>
#include "mkl.h"
#define min(x,y) (((x) < (y)) ? (x) : (y))
int main()
{
double* A, * B, * C;
int m, n, k, i, j;
double alpha, beta;
m = 2000, k = 200, n = 1000;
alpha = 1.0; beta = 0.0;
A = (double*)mkl_malloc(m * k * sizeof(double), 64);
B = (double*)mkl_malloc(k * n * sizeof(double), 64);
C = (double*)mkl_malloc(m * n * sizeof(double), 64);
if (A == NULL || B == NULL || C == NULL) {
mkl_free(A);
mkl_free(B);
mkl_free(C);
return 1;
}
for (i = 0; i < (m * k); i++) {
A[i] = (double)(i + 1);
}
for (i = 0; i < (k * n); i++) {
B[i] = (double)(-i - 1);
}
for (i = 0; i < (m * n); i++) {
C[i] = 0.0;
}
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
m, n, k, alpha, A, k, B, n, beta, C, n);
for (i = 0; i < min(m, 6); i++) {
for (j = 0; j < min(k, 6); j++) {
printf("%12.0f", A[j + i * k]);
}
printf("\n");
}
for (i = 0; i < min(k, 6); i++) {
for (j = 0; j < min(n, 6); j++) {
printf("%12.0f", B[j + i * n]);
}
printf("\n");
}
for (i = 0; i < min(m, 6); i++) {
for (j = 0; j < min(n, 6); j++) {
printf("%12.5G", C[j + i * n]);
}
printf("\n");
}
mkl_free(A);
mkl_free(B);
mkl_free(C);
return 0;
}