如何将自定义Matrix类传入LAPACK子程序？

算力阿峰

2026-7-13

问题

我想借助LAPACK库实现自己写起来复杂的功能，但没法把自定义的Matrix类传入LAPACK的特征值求解函数里得到预期结果。

我的Matrix类核心实现代码如下：

#ifndef MATRIX_T_H
#define MATRIX_T_H

#include <iostream>
#include <vector>
#include <functional>
#include <fstream>
#include "Complex.hpp"
#include "abs.hpp"

#define This (*this)
#define cout std::cout
#define endl std::endl
#define pa(a) std::pair<a>
#define O(a) std::optional<a>


template <class T>
class Column;

template <class T>
class Row;

template<class T>
struct pivot {
    uint index;
    T val;
};

template<class T>
class Matrix
{
public:
//**CONSTRUCTORS AND DESTRUCTORS*********************************************

    Matrix() : m_rows(0), m_cols(0) {}
    explicit Matrix(int nrows, int ncols=0);
    explicit Matrix (const T& val, uint nrows, uint ncols=0);
    explicit Matrix(T (&fun) (int,int), uint nrows, uint ncols=0);
    explicit Matrix(std::function<T(int,int)>& fun, uint nrows, uint ncols=0)
        : m_rows(nrows), m_cols(ncols==0 ? nrows : ncols)
    {
        createMatrix();
        fill(fun);
    }
    Matrix(const std::vector<std::vector<T>>& matrix);
    Matrix(const Matrix<T>& other); //Copy constructor
    Matrix(Matrix<T>&& other);      //Move constructor
    ~Matrix();
//****************************************************************************
//*****************  methods *************************************************
    void createMatrix(uint size) {
        m_rows = size;
        m_cols = size;
        createMatrix();
    }
    void createMatrix(uint n_rows, uint n_cols) {
        m_rows = n_rows;
        m_cols = n_cols;
        createMatrix();
    }
    T **getMatrix() {return m_matrix;}
    T *getLinearMatrix() {return m_line_matrix;}


//*****************  operators *******************************************
    operator T** () {return m_matrix;}
    operator T* () {return m_line_matrix;}

// MORE CODE

//********************************************************************


protected:
    T **m_matrix = nullptr;
    T *m_line_matrix = nullptr;
    bool m_created_matrix = false;
    uint m_rows;
    uint m_cols;

protected:
    void createMatrix();
    void deleteMatrix();
};


template<class T>
Matrix<T>::Matrix(int nrows, int ncols)
    : m_rows(nrows), m_cols(ncols==0 ? nrows : ncols) {
    createMatrix();
}

template<class T>
Matrix<T>::Matrix(const T &val, uint nrows, uint ncols)
    : m_rows(nrows), m_cols(ncols==0 ? nrows : ncols) {
    createMatrix();
    for (uint i=0; i<m_rows; i++)
        m_matrix[i][i] = val;
}

template<class T>
Matrix<T>::Matrix(T (&fun)(int, int), uint nrows, uint ncols)
    : m_rows(nrows), m_cols(ncols==0 ? nrows : ncols) {
    createMatrix();
    fill(fun);
}

template<class T>
Matrix<T>::Matrix(const std::vector<std::vector<T>> &matrix) : m_rows(matrix.size()), m_cols(matrix[0].size()) {
    createMatrix();
    for (int i=0; i<m_rows; i++)
        for (int j=0; j<m_cols; j++)
            m_matrix[i][j] = matrix[i][j];
}

template<class T>
Matrix<T>::Matrix(const Matrix<T> &other) : m_rows(other.rows()), m_cols(other.cols()) {
//    cout << "Matrix::copy_constructor" << endl;
    createMatrix();
    copyMatrix(other.m_matrix);
}

template<class T>
Matrix<T>::Matrix(Matrix<T> &&other) : m_created_matrix(true), m_rows(other.m_rows), m_cols(other.m_cols){
    delete[] m_matrix;
    delete[] m_line_matrix;
    m_matrix = other.m_matrix;
    m_line_matrix = other.m_line_matrix;
    other.m_line_matrix = nullptr;
    other.m_matrix = nullptr;
    other.m_created_matrix = false;
}

template<class T>
Matrix<T>::~Matrix() {
    deleteMatrix();
}


template<class T>
T *&Matrix<T>::operator [](int index) const {
    return m_matrix[index];
}


template<class T>
void Matrix<T>::createMatrix() {
    if (m_line_matrix!=nullptr || m_matrix != nullptr)
        deleteMatrix();

    m_line_matrix = new T[m_rows*m_cols]();
    if (m_line_matrix==0) {
        cout << "Could not allocate new memory. Needed " << m_rows*m_cols*sizeof(T) << " bytes" << endl;
        abort();
    }
    m_matrix = new T*[m_rows]();

    for (uint i=0; i<m_rows; i++)
        m_matrix[i] = m_line_matrix + m_cols*i;

    m_created_matrix = true;

}

template<class T>
void Matrix<T>::deleteMatrix() {
    delete [] m_line_matrix;
    delete[] m_matrix;
    m_line_matrix = nullptr;
    m_matrix = nullptr;
    m_created_matrix = false;
}

// MORE CODE
}

#undef This
#undef cout
#undef endl
#undef pa
#undef O

#endif // MATRIX_T_H

从createMatrix()能看到，我用动态分配的连续数组m_line_matrix，再通过指针数组m_matrix模拟二维矩阵，[][]运算符能正常访问元素。

我要调用的LAPACK特征值求解函数签名是：

void LAPACK_dsygv(
    lapack_int const* itype, char const* jobz, char const* uplo,
    lapack_int const* n,
    double* A, lapack_int const* lda,
    double* B, lapack_int const* ldb,
    double* W,
    double* work, lapack_int const* lwork,
    lapack_int* info );

其中lapack_int就是int类型。我试了三种传参方式：

使用getLinearMatrix()传递连续数组指针：

dsygv_(&itype, &jobz , &uplo , &ndim, amm.getLinearMatrix(), &nnn , axx.getLinearMatrix(), &nnn, &wr[0], &aux[0], &lwork , &info);

这里amm和axx是Matrix类实例，aux和wr是std::vector实例。

取矩阵第一个元素的地址：

dsygv_(&itype, &jobz , &uplo , &ndim, &(amm[0][0]), &nnn , &(axx[0][0]), &nnn, &wr[0], &aux[0], &lwork , &info);

变量同上。

使用静态二维数组：
```
dsygv_(&itype, &jobz , &uplo , &ndim, *amm, &nnn , *axx, &nnn, wr, aux, &lwork , &info);
```
这里amm和axx是double amm[NMAX][NMAX]静态数组，wr和aux是静态数组。

三种方式里矩阵和数组的数值都一致，但只有第三种能得到正确结果。我推测问题出在LAPACK无法识别自定义Matrix类的内存布局，请问该怎么解决？

解决方案

核心问题

LAPACK的dsygv函数要求输入矩阵采用列优先存储（Fortran默认存储方式），而你的Matrix类当前是行优先存储：m_matrix[i] = m_line_matrix + m_cols*i让每行元素在内存中连续排列，和LAPACK的要求不匹配。静态二维数组在C++中默认按列优先存储，所以第三种传参方式能正常工作。

修复方案

方案1：修改Matrix类为列优先存储

这是最彻底的解决方案，让类的内存布局直接匹配LAPACK要求：

调整createMatrix()的指针映射逻辑：

template<class T>
void Matrix<T>::createMatrix() {
    if (m_line_matrix!=nullptr || m_matrix != nullptr)
        deleteMatrix();

    m_line_matrix = new T[m_rows*m_cols]();
    if (m_line_matrix==0) {
        cout << "Could not allocate new memory. Needed " << m_rows*m_cols*sizeof(T) << " bytes" << endl;
        abort();
    }
    m_matrix = new T*[m_rows]();

    // 列优先映射：每行指针指向连续内存中该行的起始位置，元素间隔为行数
    for (uint i=0; i<m_rows; i++)
        m_matrix[i] = m_line_matrix + i;

    m_created_matrix = true;
}

修正元素访问逻辑：
替换原有的operator[]为更明确的行列访问运算符，避免索引混淆：

template<class T>
T& Matrix<T>::operator()(int row, int col) const {
    // 列优先的索引计算公式：列号 * 行数 + 行号
    return m_line_matrix[col * m_rows + row];
}

传递正确的lda参数：
lda（主维度）在列优先存储中等于矩阵的行数，调用时传入m_rows而非列数：

lapack_int ndim = amm.m_rows;
lapack_int lda = ndim;
lapack_int ldb = ndim;
dsygv_(&itype, &jobz, &uplo, &ndim, amm.getLinearMatrix(), &lda, axx.getLinearMatrix(), &ldb, wr.data(), aux.data(), &lwork, &info);

方案2：不修改存储布局，临时转置矩阵

如果不想改动原有类的存储方式，可以在调用LAPACK前将行优先矩阵转置为列优先的临时数组，调用完成后再转置结果（若需要特征向量）：

// 转置行优先矩阵到列优先临时数组
std::vector<double> temp_A(amm.m_rows * amm.m_cols);
for (int i = 0; i < amm.m_rows; ++i) {
    for (int j = 0; j < amm.m_cols; ++j) {
        temp_A[j * amm.m_rows + i] = amm[i][j];
    }
}

// 调用LAPACK，lda传入矩阵行数
lapack_int ndim = amm.m_rows;
lapack_int lda = ndim;
dsygv_(&itype, &jobz, &uplo, &ndim, temp_A.data(), &lda, ...);

// 若需要特征向量，将结果转置回行优先
if (jobz == 'V') {
    for (int i = 0; i < ndim; ++i) {
        for (int j = 0; j < ndim; ++j) {
            amm[i][j] = temp_A[j * ndim + i];
        }
    }
}

验证内存布局

可以打印连续内存的元素顺序，确认是否为列优先：

Matrix<double> mat(2,2);
mat(0,0) = 1; mat(0,1) = 2;
mat(1,0) = 3; mat(1,1) = 4;

for (int i = 0; i < 4; ++i) {
    std::cout << mat.getLinearMatrix()[i] << " "; // 列优先应输出：1 3 2 4
}

内容的提问来源于stack exchange，提问作者Alessandro

如何将自定义Matrix类传入LAPACK子程序？

核心问题 #

修复方案 #

方案1：修改Matrix类为列优先存储 #

方案2：不修改存储布局，临时转置矩阵 #

验证内存布局 #

核心问题

修复方案

方案1：修改Matrix类为列优先存储

方案2：不修改存储布局，临时转置矩阵

验证内存布局