【发布时间】：2021-01-16 01:10:00
【问题描述】：
我正在为一个大学项目编写自己的库，其中包含模板类：Vector 和 Matrix。除了这些模板类之外，还有向量和矩阵的相关模板函数。教授明确告诉我们将矩阵定义为一维数组，其中元素按列排序（效率/优化的原因）。 “矩阵”模板类有3个模板参数：矩阵允许的数据类型、行数、列数。
template <class T, unsigned int M, unsigned int N>
class Matrix
话虽如此，我立即解决了问题。我正在编写一个函数，该函数使用列的 LaPlace 规则（使用第一列）计算任何维度 > 4 的矩阵的行列式。 我还编写了一个用于二维矩阵的函数（称为 D2MatrixDet）和一个用于测试和工作的三维矩阵的函数（称为 D3MatrixDet）：
template <class T>
double D2MatrixDet(const Matrix<T, 2, 2>& _m)
template <class T>
double D3MatrixDet(const Matrix<T, 3, 3>& _m)
我要写的模板函数有两个模板参数：输入矩阵的数据类型，矩阵的维数（因为行列式是为方阵计算的，只有一个维就够了）。它是一个递归函数；变量“结果”是在每一步将行列式保存在内存中的变量。下面是我写的代码。
template <class T, unsigned int D>
void DNMatrixDet(Matrix<T, D, D> _m, double result) //LaPlace Rule respect to the first column
{
const unsigned int new_D = D - 1;
Matrix<T, new_D, new_D> temp;
if (D > 3)
{
for (unsigned int i = 0; i < _m.row; ++i)
//Indicate the element to multiply
{
for (unsigned int j = _m.row, l = 0; j < _m.row * _m.column && l < pow(new_D, 2); ++j)
//Manage the element to be inserted in temp
{
bool invalid_row = false;
for (unsigned int k = 1; k < _m.row && invalid_row == false; ++k) //Slide over row
{
if (j == (i + k * _m.row))
{
invalid_row = true;
}
}
if (invalid_row == false)
{
temp.components[l] = _m.components[j];
++l;
}
}
DNMatrixDet(temp, result);
result += pow((-1), i) * _m.components[i] * result;
}
}
else if (D == 3)
{
result += D3MatrixDet(_m);
}
}
在 main 中，我使用 5 x 5 矩阵测试函数。 当我尝试编译时，出现了几个错误，都非常相似，并且与矩阵的大小有关，每一步都会减少一个。这是初始矩阵大小为 5 时（LA 是库的名称，Test.cpp 是包含 main 的文件）：
LA.h: In instantiation of 'void LA::DNMatrixDet(LA::Matrix<T, M, M>, double) [with T = double;
unsigned int D = 5]':
Test.cpp:437:33: required from here
LA.h:668:34: error: no matching function for call to 'D3MatrixDet(LA::Matrix<double, 5, 5>&)'
result += D3MatrixDet(_m);
~~~~~~~~~~~^~~~
In file included from Test.cpp:1:
LA.h:619:12: note: candidate: 'template<class T> double LA::D3MatrixDet(const LA::Matrix<T, 3, 3>&)'
double D3MatrixDet(const Matrix<T, 3, 3>& _m)
^~~~~~~~~~~
LA.h:619:12: note: template argument deduction/substitution failed:
In file included from Test.cpp:1:
LA.h:668:34: note: template argument '5' does not match '3'
result += D3MatrixDet(_m);
~~~~~~~~~~~^~~~
这是大小变成4的时候：
LA.h: In instantiation of 'void LA::DNMatrixDet(LA::Matrix<T, M, M>, double) [with T = double;
unsigned int D = 4]':
LA.h:662:28: required from 'void LA::DNMatrixDet(LA::Matrix<T, M, M>, double) [with T = double;
unsigned int D = 5]'
Test.cpp:437:33: required from here
LA.h:668:34: error: no matching function for call to 'D3MatrixDet(LA::Matrix<double, 4, 4>&)'
In file included from Test.cpp:1:
LA.h:619:12: note: candidate: 'template<class T> double LA::D3MatrixDet(const LA::Matrix<T, 3, 3>&)'
double D3MatrixDet(const Matrix<T, 3, 3>& _m)
^~~~~~~~~~~
LA.h:619:12: note: template argument deduction/substitution failed:
In file included from Test.cpp:1:
LA.h:668:34: note: template argument '4' does not match '3'
result += D3MatrixDet(_m);
~~~~~~~~~~~^~~~
等等。它一直在下降，直到从 4294967295 重新开始（我发现这是 32 位“无符号整数”的上限）并继续下降，直到达到模板实例的最大数量（= 900）。
在每次迭代中，编译器总是检查用于计算 3 x 3 行列式的函数，即使该函数仅在输入矩阵为 3 x 3 时才执行。那么它为什么要检查在理论不应该发生？
我多次仔细检查了我所写内容的数学逻辑，即使借助写在纸上的矩阵并缓慢执行第一步也是如此。我相信并希望这是对的。我很确定问题与使用模板和递归函数有关。
对于这个很长的问题，我深表歉意，我试图以最好的方式解释它。我希望我已经很好地解释了这个问题。
编辑： 通过在 DNMatrixDet 函数的开头定义“if constexpr”来解决问题。编译成功。我只需要修复算法，但这超出了帖子的范围。以下是所做更改的代表：
template <class T, unsigned int M, unsigned int N>
class Matrix
{
public:
T components[M * N];
unsigned int row = M;
unsigned int column = N;
Matrix()
{
for (unsigned int i = 0; i < M * N; ++i)
{
components[i] = 1;
}
}
Matrix(T* _c)
{
for (unsigned int i = 0; i < M * N; ++i, ++_c)
{
components[i] = *_c;
}
}
friend std::ostream& operator<<(std::ostream& output, const Matrix& _m)
{
output << _m.row << " x " << _m.column << " matrix:" << std::endl;
for (unsigned int i = 0; i < _m.row; ++i)
{
for (unsigned int j = 0; j < _m.column; ++j)
{
if (j == _m.column -1)
{
output << _m.components[i + j*_m.row];
}
else
{
output << _m.components[i + j*_m.row] << "\t";
}
}
output << std::endl;
}
return output;
}
};
template <class T>
double D3MatrixDet(const Matrix<T, 3, 3>& _m)
{
double result = _m.components[0] * _m.components[4] * _m.components[8] +
_m.components[3] * _m.components[7] * _m.components[2] +
_m.components[6] * _m.components[1] * _m.components[5] -
(_m.components[6] * _m.components[4] * _m.components[2] +
_m.components[3] * _m.components[1] * _m.components[8] +
_m.components[0] * _m.components[7] * _m.components[5]);
return result;
}
template <class T, unsigned int D>
void DNMatrixDet(Matrix<T, D, D> _m, double result)
{
Matrix<T, D - 1, D - 1> temp;
if constexpr (D > 3)
{
for (unsigned int i = 0; i < D; ++i)
{
for (unsigned int j = D, l = 0; j < D * D && l < (D - 1) * (D - 1); ++j)
{
bool invalid_row = false;
for (unsigned int k = 1; k < D && invalid_row == false; ++k)
{
if (j == (i + k * D))
{
invalid_row = true;
}
}
if (invalid_row == false)
{
temp.components[l] = _m.components[j];
++l;
}
}
DNMatrixDet(temp, result);
result += i & 1 ? -1 : 1 * _m.components[i] * result;
}
}
else if (D == 3)
{
result += D3MatrixDet(_m);
}
}
int main()
{
double m_start[25] = {4, 9, 3, 20, 7, 10, 9, 50, 81, 7, 20, 1, 36, 98, 4, 20, 1, 8, 5, 93, 47, 21, 49, 36, 92};
Matrix<double, 5, 5> m = Matrix<double, 5, 5> (m_start);
double m_det = 0;
DNMatrixDet(m, m_det);
std::cout << "m is " << m << std::endl;
std::cout << "Det of m is " << m_det << std::endl;
return 0;
}
【问题讨论】：
-
我尝试重构您的代码，但仅模板（删除了一些内容）不会触发该错误 (godbolt.org/z/q6czGa)
-
在我看来，您使用带有整数参数的
pow()的所有地方都非常可疑。pow(-1, i)可以写成例如i & 1 ? -1 : 1甚至是1 - 2 * (i & 1)。pow(i, 2)更好地实现为i * i。请注意，std::pow() 仅用于浮点数，可能会引入整数类型的奇怪舍入问题。 -
如果你需要 constexpr。
标签： c++ templates recursion matrix linear-algebra