c – boost :: ublas如何获得int矩阵的行列式？

2023年6月16日 238次阅读

我找到了计算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;
}

这个函数工作正常,但只有非整数类型(它适用于float和double).当我尝试传递相同的矩阵但类型为int时,我收到了编译错误：

Check failed in file c:\boost\boost_1_58_0\boost\numeric\ublas\lu.hpp at line 167:
singular != 0 || detail::expression_type_check (prod (triangular_adaptor (m), triangular_adaptor (m)), cm)
unknown location(0): fatal error in “BaseTest”: std::logic_error: internal logic

这是提升的错误还是我的功能错了？
我可以做些什么改变来避免这个错误？

最佳答案 当然这不是一个错误,因为这样的检查完全是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

看一下表,其中平均时间是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它更快,但事实并非如此.平均来说,对于更多的运行,我确信你会得到轻微的浮动加速度(我猜大约0-15毫秒,这是时间测量误差).另外,如果我们测量复制时间,我们将看到,对于小于3000 x 3000的矩阵大小,它接近于零(它也是大约0-15 ms,这是时间测量误差).对于较大的矩阵,复制时间增加(例如,对于5000 x 5000矩阵,它是125 ms).但有一个重要的注意事项！复制时间小于int型矩阵的所有行列式计算的1.0％,并随着大小的增长而大大减少！

所有这些措施都是针对使用Visual Studio 2013在发布配置中使用boost版本1.58编译的程序而制定的.时间用时钟功能测量. CPU是Intel Xeon E5645 2.4GHz.

此外,性能主要取决于您的方法的使用方式.如果您要为小矩阵多次调用此函数,则复制时间可能会变得很大.

点赞