我找到了计算boost:：ublas矩阵的行列式的函数：
template<typename ValType>
ValType det_fast(const ublas::matrix<ValType>& matrix)
{
// create a working copy of the input
ublas::matrix<ValType> mLu(matrix);
ublas::permutation_matrix<std::size_t> pivots(matrix.size1());
auto isSingular = ublas::lu_factorize(mLu, pivots);
if (isSingular)
return static_cast<ValType>(0);
ValType det = static_cast<ValType>(1);
for (std::size_t i = 0; i < pivots.size(); ++i)
{
if (pivots(i) != i)
det *= static_cast<ValType>(-1);
det *= mLu(i, i);
}
return det;
}
此函数运行良好，但仅适用于非整数类型(它适用于浮点和双精度)。当我试图传递相同的矩阵，但类型为int时，我收到编译错误：
在文件c:\boost\boost_1_58_0\boost\numeric\ublas\lu.hpp的第167行检查失败：单数！=0|| detail:：expression_type_check(prod(triangular_adaptor(m)，triangular_adaptor(m))，cm)未知位置(0)："BaseTest"中的致命错误：std:：logic_error:内部逻辑
是boost的错误还是我的函数错误？我可以做些什么来避免这个错误？
当然这不是一个bug，因为像这样的检查遍布uBLAS。我想这是因为它的大多数操作对非浮点类型没有意义。
您可以使用禁用类型检查
#define BOOST_UBLAS_TYPE_CHECK 0
包含之前。但要三思！看看的例子
#include <iostream>
#define BOOST_UBLAS_TYPE_CHECK 0
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
namespace ublas = boost::numeric::ublas;
template<typename ValType>
ValType det_fast(const ublas::matrix<ValType>& matrix)
{
// create a working copy of the input
ublas::matrix<ValType> mLu(matrix);
ublas::permutation_matrix<std::size_t> pivots(matrix.size1());
auto isSingular = ublas::lu_factorize(mLu, pivots);
if (isSingular)
return static_cast<ValType>(0);
ValType det = static_cast<ValType>(1);
for (std::size_t i = 0; i < pivots.size(); ++i)
{
if (pivots(i) != i)
det *= static_cast<ValType>(-1);
det *= mLu(i, i);
}
return det;
}
int main()
{
ublas::matrix<int> m (3, 3);
for (unsigned i = 0; i < m.size1 (); ++ i)
for (unsigned j = 0; j < m.size2 (); ++ j)
m (i, j) = 3 * i + j;
std::cout << det_fast(m) << std::endl;
}
很明显，矩阵m是奇异的，如果将类型从int更改为double，函数将返回零。但是对于int，它返回-48。
编辑#1
此外，您可以在没有任何错误的情况下创建ublas::matrix<int>，并将其分配给ublas::matrix<float>。因此，正确计算行列式的一种方法是重载det_fast并删除define语句。
int det_fast(const ublas::matrix<int>& matrix)
{
return (int)det_fast(ublas::matrix<float>(matrix));
}
编辑#2
看一下表，其中average time是5个程序运行的全行列式计算时间(对于带有复制操作的int)。
size | average time, ms
------------------------
| int float
------------------------
100 | 9.0 9.0
200 | 46.6 46.8
300 | 156.4 159.4
400 | 337.4 331.4
500 | 590.8 609.2
600 | 928.2 1009.4
700 | 1493.8 1525.2
800 | 2162.8 2231.0
900 | 3184.2 3025.2
你可能认为int更快，但事实并非如此。平均而言，对于更多的跑步，我相信float会有轻微的加速度(我猜大约是0-15毫秒，这是时间测量误差)。此外，如果我们测量复制时间，我们会看到，对于小于3000 x 3000的矩阵大小，它接近于零(大约为0-15毫秒，这是时间测量误差)。对于较大的矩阵，复制时间会增加(例如，对于5000 x 5000矩阵，它是125毫秒)。但有一点很重要！int型矩阵的复制时间小于所有行列式计算的1.0%，并且随着大小的增长而大大减少！
所有这些措施都是针对一个使用Visual Studio 2013编译的程序制定的，该程序的发布配置为boost版本1.58。用CCD_ 15函数测量时间。CPU为Intel Xeon E5645 2.4GHz。
此外，性能主要取决于如何使用您的方法。如果你要为小矩阵多次调用这个函数，那么复制时间可能会变得很重要。