高速缓存与矩阵乘法(三)

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);

这样使用起来就方便许多，并且计算的效率也比直接按公式展开要高。

注：

1. 这里仅给出矩阵相乘和求解线性方程组的封装函数接口，如果需要调用 CLAPACK 库中的其它函数接口，也可按照类似的方式进行封装。
2. 完整的测试程序见附件，附件中还给出与直接按公式展开的计算方式的对比测试程序。