将自定义 Matrix 类传递给 Lapack 子例程

问题描述 投票:0回答:1

我想要 lapack 包来实现一些非常有用的功能,但我不喜欢/无法自己很好地实现这些功能。问题是我无法将自定义 Matrix 类传递给 lapack 函数,我没有得到想要的行为。这是我的矩阵类

#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()
中看到的那样,我动态创建了一个数组,并通过将地址保存在指向指针的指针中来模拟矩阵结构。因此,从技术上讲,使用双运算符 [][] 给出了矩阵所需的元素。

问题是当我尝试将其传递给我感兴趣的 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 值。 我尝试以多种不同的方式传递我的矩阵,i.e.

  • dsygv_(&itype, &jobz , &uplo , &ndim, *amm.getLinearMatrix(), &nnn , *axx.getLinearMatrix(), &nnn, &wr[0], &aux[0], &lwork , &info);
    ,其中 ammaxx 是 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);
    ,现在我有
    double amm[NMAX][NMAX]
    和 axx 相同,而我有
    int wr[NMAX], aux[4*NMAX]

在所有情况下,矩阵和数组的值都是相同的,但仅在最后一种情况下才有效。我想这与该函数无法将矩阵[][]识别为双精度**有关。

我该如何解决这个问题?

正如我所说,我尝试改变传递参数的方式,但无法直接访问 lapack 库代码(而且我不想手动重新编译所有内容),我不能只使用 gdb 或其他工具来检查已发生的内容传递给函数。

c++ matrix parameter-passing lapack double-pointer
1个回答
0
投票
  • 要修复第一个:
    dsygv_(&itype, &jobz , &uplo , &ndim, amm.getLinearMatrix(), &nnn , axx.getLinearMatrix(), &nnn, &wr[0], &aux[0], &lwork , &info)
    ,因为
    *amm.getLinearMatrix()
    *axx.getLinearMatrix()
    的类型为
    double
    ,而不是
    double *
    。这里应该是编译错误吧?
  • 修复第二个
    dsygv_(&itype, &jobz , &uplo , &ndim, amm[0], &nnn , axx[0], &nnn, &wr[0], &aux[0], &lwork , &info);
    ,因为
    &(amm[0][0])
    是堆栈上临时
    double
    对象
    amm[0][0]
    的地址,因此它的值不等于
    amm.getLinearMatrix()
© www.soinside.com 2019 - 2024. All rights reserved.