矩阵运算库-Cblas
CBLAS
CBlas是BLAS(Basic Linear Algebra Subprograms)的C接口，提供基本的矢量和矩阵计算。其分为三个级别的计算：
- Level 1：标量、矢量、矢量-矢量计算
- Level 2：矩阵-矢量计算
- Level 3：矩阵-矩阵计算
1 实数矩阵-矢量乘法(Level 2)
双精度带状矩阵-矢量乘法
计算公式为：
\[\begin{aligned}
y:=& \alpha \bold{A} \times \bold{x} + \beta \bold{y} \\
y:=&\alpha \bold{A}^{T} \times \bold{x} + \beta \bold{y}
\end{aligned}
\]
此处 \(\alpha\) 和 \(\beta\) 为标量，\(x\) 和 \(y\) 为 \(n\) 元素的矢量，\(A\) 为 $m \times n $ 矩阵，且具有 \(kl\) 个下对角，\(ku\) 个上对角。
若\(A\) 的形式\(n\times n\)的带状方阵：
\[A = \left(\begin{matrix}
a_{11} & a_{12} & a_{13} & a_{14} & 0 & 0 & 0 & 0 \\
a_{21} & a_{22} & a_{23} & a_{24} & a_{25} & 0 & 0 & 0 \\
a_{31} & a_{32} & a_{33} & a_{34} & a_{35} & a_{36} & 0 & 0 \\
0 & a_{42} & a_{43} & a_{44} & a_{45} & a_{46} & a_{47} & 0 \\
0 & 0 & a_{53} & a_{54} & a_{55} & a_{56} & a_{57} & a_{58} \\
0 & 0 & 0 & a_{64} & a_{65} & a_{66} & a_{67} & a_{68} \\
0 & 0 & 0 & 0 & a_{75} & a_{76} & a_{77} & a_{78} \\
0 & 0 & 0 & 0 & 0 & a_{86} & a_{87} & a_{88} \\
\end{matrix}\right)
\]
则 \(A\) 的上带宽，也就是上对角个数为3；下带宽，也就是下对角个数为2。将其转化为带状矩阵的存储，用 \(A.band\) 表示，转化为\((kl + ku + 1)\times n\) 的矩阵，如下：
\[A.band = \left(\begin{matrix}
0 & 0 & 0 & a_{14} & a_{25} & a_{36} & a_{47} & a_{58} \\
0 & 0 & a_{13} & a_{24} & a_{35} & a_{46} & a_{57} & a_{68} \\
0 & a_{12} & a_{23} & a_{34} & a_{45} & a_{56} & a_{67} & a_{78} \\
a_{11} & a_{22} & a_{33} & a_{44} & a_{55} & a_{66} & a_{77} & a_{88} \\
a_{21} & a_{32} & a_{43} & a_{54} & a_{65} & a_{76} & a_{87} & 0 \\
a_{31} & a_{42} & a_{53} & a_{64} & a_{75} & a_{86} & 0 & 0 \\
\end{matrix}\right)
\]
在所调用的函数中，我们需要输入的是 \(A.band\) 矩阵的转置，即输入
\[(A.band)^{T} = \left(\begin{matrix}
0 & 0 & 0 & a_{11} & a_{21} & a_{31} \\
0 & 0 & a_{12} & a_{22} & a_{32} & a_{42} \\
0 & a_{13} & a_{23} & a_{33} & a_{43} & a_{53} \\
a_{14} & a_{24} & a_{34} & a_{44} & a_{54} & a_{64} \\
a_{25} & a_{35} & a_{45} & a_{54} & a_{65} & a_{75} \\
a_{36} & a_{46} & a_{56} & a_{66} & a_{76} & a_{86} \\
a_{47} & a_{57} & a_{67} & a_{77} & a_{87} & 0 \\
a_{58} & a_{68} & a_{78} & a_{88} & 0 & 0 \\
\end{matrix}\right)
\]
注意到，对角元在 \(ku+1\) 列，\(ku\)、\(ku-1\)、\(\cdots\) 分别为第一个上对角、第二个上对角... ，依次类推；而 \(ku+2\)、\(ku+3\) 对应到第一个下对角、第二个下对角...，依次类推。
调用函数为：clbas_dgbmv()
void cblas_dgbmv(const enum CBLAS_ORDER order,
const enum CBLAS_TRANSPOSE TransA,
const int M,
const int N,
const int KL,
const int KU,
const double alpha,
const double *A,
const int lda,
const double *X,
const int incX,
const double beta,
double *Y,
const int incY
);
example
计算\(4\times 4\)矩阵与矢量的运算：
\[A \times x =
\left(\begin{matrix}
1 & 2 & 0 & 0 \\
3 & 4 & 5 & 0 \\
0 & 6 & 7 & 8 \\
0 & 0 & 8 & 9
\end{matrix}\right)
\left(\begin{matrix}
1 \\ 2 \\ 3 \\ 4
\end{matrix}\right)
=
\left(\begin{matrix}
5 \\ 26 \\ 65 \\ 60
\end{matrix}\right)
\]
矩阵 \(A\) 转换 \((A.band)^{T}\) 为作为函数输入
\[(A.band)^{T} = \left(\begin{matrix}
0 & 1 & 3 \\
2 & 4 & 6 \\
5 & 7 & 9 \\
8 & 10 & 0
\end{matrix}\right)
\]
具体程序如下：
编译命令：g++ -o test test.cpp -I$HOME/local/include -L$HOME/local/lib -lcblas -lblas
#include <iostream>
#include <complex>
#include <cblas.h>
using std::complex;
using std::cout;
using std::cin;
using std::endl;
int main() {
double a[4][3];
a[0][0] = 0., a[0][1] = 1., a[0][2] = 3.0;
a[1][0] = 2., a[1][1] = 4., a[1][2] = 6.0;
a[2][0] = 5., a[2][1] = 7., a[2][2] = 9.0;
a[3][0] = 8., a[3][1] = 10., a[3][2] = 0.0;
double x[4] = {1.0, 2.0, 3.0, 4.0};
double y[4] = {0.};
cblas_dgbmv(CblasRowMajor,CblasNoTrans, 4, 4, 1, 1, 1, *a, 3, x, 1, 1.0, y, 1);
for(int i = 0; i < 4; i++)
{
cout << y[i] << endl;
}
return 0;
}
2. 复数矩阵-矢量乘法(Level 2)
在Cblas中，对于复数的存储，奇数位存储实数部分，偶数位存储虚数部分。
一般双精复数类型矩阵-矢量
计算公式为：
\[y:= \alpha \bold{A} \times \bold{x} + \beta \bold{y} \\
y:=\alpha \bold{A}^{T} \times \bold{x} + \beta \bold{y} \\
y:=\alpha \bold{A}^{H} \times \bold{x} + \beta \bold{y}
\]
调用的函数为clbas_zgemv()
void cblas_zgemv(const enum CBLAS_ORDER order, // CblasColMajor 列主序; CblasRawMajor 行主序
const enum CBLAS_TRANSPOSE TransA, //
const int M, // 矩阵A的行数，至少为0
const int N, // 矩阵A的列数，至少为0
const void *alpha, // 标量，复数形式，公式中的alpha
const void *A, // 复数矩阵，维度为（lda，N）
const int lda, // A的第一个维度，取MAX(1, m)
const void *X, // 增量矩阵，
// 矩阵A未转置：维度至少为：(1+(n-1))*abs(incX)
// 矩阵A转置后，维度至少为：(1+(m-1))*abs(incX)
const int incX, // X的元素增量，不能为0
const void *beta, // 标量，复数形式，公式中的beta，若为0，无需设置Y
void *Y, // 输入量及输出量，Y为最后的累加形式
// 矩阵A未转置：维度至少为：(1+(m-1))*abs(incY)
// 矩阵A转置后，维度至少为：(1+(n-1))*abs(incY)
const int incY // Y的元素增量，不可为0
);
双精复数厄密共轭矩阵-矢量
计算公式为：
\[y:= \alpha \bold{A} \cdot \bold{x} + \beta \bold{y}
\]
此处 \(\alpha\) 和 \(\beta\) 为标量，\(x\) 和 \(y\) 为 \(n\) 元素的矢量，\(A\) 为 $n \times n $ 矩阵。
调用函数cblas_zhemv：
void cblas_zhemv(const enum CBLAS_ORDER order, // 列主序或行主序
const enum CBLAS_UPLO Uplo, // A矩阵的上三角或下三角部分
// CblasUpper：引用上三角部分
// CblasLower：引用下三角部分
const int N, // 矩阵A的位数，至少为0
const void *alpha,
const void *A, // 若不设置对角虚部，则视为0
const int lda, // A的第一个维度，取MAX(1, n)
const void *X, // 维度至少为(1+(n-1))*abs(incX)
const int incX, // X的元素增量,不可为0，
const void *beta, //
void *Y, // 维度至少为(1+(n-1))*abs(incY)
const int incY // Y的元素增量，不可为0，意思是每隔几个元素乘一个beta
);
example
编译命令：g++ -o test test.cpp -I$HOME/local/include -L$HOME/local/lib -lcblas -lblas
// 计算矩阵—矢量乘法 y = a * x
#include <iostream>
#include <complex>
#include <cblas.h>
int main() {
using namespace std;
typedef complex<double> Comp; // 定义复数类型
int n = 2;
Comp a[n][n] = {{Comp(1, 0.), Comp(1, 2)},
{Comp(1, -2), Comp(2, 0.)}};
Comp x[n] = {Comp(1, 2), Comp(3, 4)};
Comp y[n] = {Comp(0, 0), Comp(0, 0)};
Comp alpha(1., 0.), beta(0, 0);
cblas_zhemv(CblasRowMajor, CblasUpper, n, &alpha, &a, 2, &x, 1., &beta, &y, 1.);
cout << "matrix a: " << endl
<< a[0][0] << "\t" << a[0][1] << endl
<< a[1][0] << "\t" << a[1][1] << endl;
cout << "vector x: " << endl
<< x[0] << endl
<< x[1] << endl;
cout << "alpha = " << alpha << endl;
cout << "beta = " << beta << endl;
cout << "vector y: " << endl
<< y[0] << endl
<< y[1] << endl;
}
带状复数矩阵-矢量乘法
计算公式为：
\[y:= \alpha \bold{A} \cdot \bold{x} + \beta \bold{y} \\
y:= \alpha \bold{A}^{T} \cdot \bold{x} + \beta \bold{y} \\
y:= \alpha \bold{A^{H}} \cdot \bold{x} + \beta \bold{y}
\]
此处 \(\alpha\) 和 \(\beta\) 为标量，\(x\) 和 \(y\) 为 \(n\) 元素的矢量，\(A\) 为包含 \(kl\) 下对角和 \(ku\) 下对角的 $m \times n $ 矩阵。
调用函数cblas_zgbmv：
void cblas_zgbmv(const enum CBLAS_ORDER order,
const enum CBLAS_TRANSPOSE TransA,
const int M, // 矩阵A的行数
const int N, // 矩阵A的列数
const int KL, // 矩阵A的下对角数，即
const int KU,
const void *alpha,
const void *A, //
const int lda,
const void *X,
const int incX,
const void *beta,
void *Y,
const int incY
);
3. 矩阵-矩阵乘法
复数对称矩阵乘法
计算公式为：
\[\bold{C}:= \alpha \bold{A} \cdot \bold{B} + \beta \bold{C}\\
\bold{C}:= \alpha \bold{B} \cdot \bold{A} + \beta \bold{C}
\]
此处 \(\alpha\) 和 \(\beta\) 为标量，\(\bold{A}\) 为 $n \times n $ 对称矩阵， \(\bold{B}\) 和 \(\bold{C}\) 为 \(m\times n\)矩阵。
void cblas_zsymm(const enum CBLAS_ORDER Order, // CblasColMajor 列主序; CblasRawMajor 行主序
const enum CBLAS_SIDE Side, // 标定A在左边或右边
const enum CBLAS_UPLO Uplo, // A 的上三角或下三角
const int M, // C 的行数
const int N, // C 的列数
const void *alpha, //
const void *A, // 维度为(lda, ka)，A在左时 ka=M，在右时ka=n
const int lda, // A的主维
const void *B,
const int ldb, // B的主维
const void *beta,
void *C,
const int ldc); // C的主维
浙公网安备 33010602011771号