CLAPACK 注意事项
在使用 CLAPACK 前,要先清楚四点:
首先是 Levels of Routines,即函数的层次,链接中有详细的介绍,这里不再赘述。
其次是 Naming Scheme,即命名规则,链接中有详细的介绍,这里不再赘述。
再次是 CLAPACK 的函数不接收二维数组,即只能用一维数组代替二维数组。例如,我想要个 array2 = {{1,2},{3,4}},那么正确的写法是 array[2*2] = {1,2,3,4}。
最后是行主序与列主序的问题。CLAPACK 看一维数组时会将其看成是按列存放的。
封装 CLAPACK
这里以使用 dgemm_ 函数为例,对 CLAPACK 的使用方法进行说明。
// 函数原型:
// dgemm_ 执行一个matrix-matrix的操作:C = ALPHA * op(A) * op(B) + BETA * C,其中ALPHA和BETA是双精度浮点数表示的系数,op(A)是m×k矩阵,op(B)是k×n矩阵,C是m×n矩阵
// 其它参数含义:
// 字符TRANSA:表示矩阵A的形式,若TRANSA为'N'或'n',则op(A)=A;否则op(A)=A'(矩阵A的转置)
// 字符TRANSB:表示矩阵B的形式,若TRANSB为'N'或'n',则op(B)=B;否则op(B)=B'(矩阵B的转置)
// 整数M:表示矩阵op(A)或矩阵op(C)的行数,并且要求M≥0
// 整数N:表示矩阵op(B)或矩阵op(C)的列数,并且要求N≥0
// 整数K:表示矩阵op(A)的列数或矩阵op(B)的行数,并且要求K≥0
// 整数LDA:表示矩阵A的行数,并且要求LDA≥1
// 整数LDB:表示矩阵B的行数,并且要求LDB≥1
// 整数LDC:表示矩阵C的行数,并且要求LDC≥1
int dgemm_(char* TRANSA, char* TRANSB, integer* M, integer*
N, integer* K, doublereal* ALPHA, doublereal* A, integer* LDA,
doublereal* B, integer* LDB, doublereal* BETA, doublereal* C,
integer* LDC);
从上述函数原型可以看出,即便是简单的两个矩阵相乘的函数也有 13 个参数,其它函数也类似。由于这些函数接口有部分参数是互相关联的,并且每次都传入所有参数对使用者来说也非常不方便。因此,这里用一个矩阵类封装二维数组,并对 CLAPACK 中的原型函数进行封装,使得函数接口方便使用。
封装矩阵类
定义一个类 DMatrix,将一维数组封装为按列存放的二维矩阵,如下
class DMatrix
{
public:
DMatrix(unsigned int row, unsigned int column)
{
unsigned int maxSize = (1 << 15);
if (row < 1) { row = 1; }
else if (row > maxSize) { row = maxSize; }
if (column < 1) { column = 1; }
else if (column > maxSize) { column = maxSize; }
m_row = row;
m_col = column;
unsigned int size = m_row * m_col;
m_data = new double[size];
memset(m_data, 0, sizeof(double) * size);
}
DMatrix(const DMatrix& matrix)
{
m_row = matrix.m_row;
m_col = matrix.m_col;
unsigned int size = m_row * m_col;
m_data = new double[size];
memcpy(m_data, matrix.head(), sizeof(double) * size);
}
DMatrix& operator=(const DMatrix& matrix)
{
if (this == &matrix) { return *this; } // 自赋值
delete[] m_data;
m_row = matrix.m_row;
m_col = matrix.m_col;
unsigned int size = m_row * m_col;
m_data = new double[size];
memcpy(m_data, matrix.head(), sizeof(double) * size);
return *this;
}
~DMatrix() { delete[] m_data; }
// 重载操作符
double& operator()(unsigned int i, unsigned int j)
{ return m_data[j * m_row + i]; }
double operator()(unsigned int i, unsigned int j) const
{ return m_data[j * m_row + i]; }
// 成员函数
double* head() const { return m_data; }
unsigned int GetRowSize() const { return m_row; }
unsigned int GetColSize() const { return m_col; }
private:
unsigned int m_row;
unsigned int m_col;
double* m_data;
};
封装矩阵运算函数
// 矩阵相乘: C = TRANSA * A * B + BETA * C => C
void my_dgemm(DMatrix& A, DMatrix& B, DMatrix& C, char TRANSA = 'N', char TRANSB = 'N', double ALPHA = 1.0, double BETA = 0.0)
{
integer M, N, K, LDA, LDB, LDC;
M = ('N' == TRANSA) ? A.GetRowSize() : A.GetColSize();
N = ('N' == TRANSB) ? B.GetColSize() : B.GetRowSize();
K = ('N' == TRANSA) ? A.GetColSize() : A.GetRowSize();
LDA = ('N' == TRANSA) ? M : K;
LDB = ('N' == TRANSB) ? K : N;
LDC = C.GetRowSize();
dgemm_(&TRANSA, &TRANSB, &M, &N, &K,
(doublereal*)&ALPHA,
(doublereal*)A.head(), &LDA,
(doublereal*)B.head(), &LDB,
(doublereal*)&BETA,
(doublereal*)C.head(), &LDC);
}
// 求解线性方程组: A * X = B => X
integer my_dgesv(DMatrix& A, DMatrix& B)
{
integer N, NRHS, LDA, LDB, INFO;
N = A.GetRowSize();
NRHS = B.GetColSize();
LDA = N;
LDB = N;
INFO = 0;
integer* IPIV = new integer[N];
dgesv_(&N, &NRHS, A.head(), &LDA, IPIV, B.head(), &LDB, &INFO);
delete[] IPIV;
return INFO;
}
封装后的矩阵相乘
// C = A * B => C
DMatrix matrixA(M, K), matrixB(K, N), matrixC(M, N);
// 矩阵matrixA赋值
// 矩阵matrixB赋值
// 计算矩阵相乘 C = A * B
my_dgemm(matrixA, matrixB, matrixC);
这样使用起来就方便许多,并且计算的效率也比直接按公式展开要高。
注:
- 这里仅给出矩阵相乘和求解线性方程组的封装函数接口,如果需要调用 CLAPACK 库中的其它函数接口,也可按照类似的方式进行封装。
- 完整的测试程序见附件,附件中还给出与直接按公式展开的计算方式的对比测试程序。