提问人:Alessandro 提问时间:8/16/2023 最后编辑:Alessandro 更新时间:8/17/2023 访问量:37
将自定义 Matrix 类传递给 Lapack 子例程
Passing a custom Matrix class to a Lapack subroutine
问:
我想要一些非常有用的功能的 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()
问题是当我尝试将其传递给我感兴趣的 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 值。 我试图以许多不同的方式传递我的矩阵,即
dsygv_(&itype, &jobz , &uplo , &ndim, amm.getLinearMatrix(), &nnn , axx.getLinearMatrix(), &nnn, &wr[0], &aux[0], &lwork , &info);
,其中 amm 和 axx 是我的 Matrix 类的实例,并且是 std::vector 的实例;aux
wr
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);
,现在我有,axx 也一样,而我有double amm[NMAX][NMAX]
int wr[NMAX], aux[4*NMAX]
在所有情况下,矩阵和数组的值都是相同的,但只有在最后一种情况下它才有效。我想这与函数不将矩阵[][]识别为双精度**这一事实有关。
我该如何解决这个问题?
正如我所说,我尝试更改传递参数的方式,但无法直接访问 lapack 库代码(而且我不想手动重新编译所有内容)我不能只使用 gdb 或其他工具来检查已传递给函数的内容。
答:
0赞
Misha T
8/16/2023
#1
- 要修复第一个:,因为 ,有类型,而不是 。这里应该是编译错误,不是吗?
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()
评论
0赞
Alessandro
8/17/2023
关于第一个修复,我实际上已经这样做了。我的帖子中有一个错别字(现已修复)。所以我知道这是行不通的。您提出的第二个修复程序没有解决它,它仍然在信息输出中给我同样的错误。
评论
double[NMAX][NMAX]