使用密集和稀疏矩阵

Posted 2023-02-21

【中文标题】使用密集和稀疏矩阵

【英文标题】：Working with both dense and sparse matrices

【发布时间】：2018-06-01 12:18:51

【问题描述】：

我正在编写一个 C++ 函数，它对矩阵进行操作，作为参数传递，并且希望代码能够处理不同类型的矩阵（例如 Boost 稀疏矩阵、std::vectors 的 std::vectors）。我目前的方法是为不同类型的基本矩阵操作定义重载方法，为不同类型的矩阵提供统一的接口，并将我的函数定义为仅使用这些重载方法的模板函数

#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <iostream>

typedef std::vector<double> vec;
typedef std::vector<vec> mat;
typedef boost::numeric::ublas::compressed_matrix<double, boost::numeric::ublas::row_major> spmat;

namespace matrix

    inline void set(spmat & input, u_int i, u_int j, double val)
    
        input(i, j) = val;
    

    inline void set(mat & input, u_int i, u_int j, double val)
    
        input[i][j] = val;
    

    inline u_int size1(const mat & input)
    
        return input.size();
    

    inline u_int size2(const mat & input)
    
        return input[0].size();
    

    inline u_int size1(const spmat & input)
    
        return input.size1();
    

    inline u_int size2(const spmat & input)
    
        return input.size2();
    

    inline double get(const spmat & input, u_int i, u_int j)
    
        return input(i, j);
    

    inline double get(const mat & input, u_int i, u_int j)
    
        return input[i][j];
    

对于简单的任务，这种方法似乎很有效。然而，目前，我正在尝试编写一个需要迭代所有条目的函数，在密集矩阵的情况下，或者在稀疏矩阵的情况下仅非零条目。我知道如何为每种情况单独执行此操作，但希望只有一个在这两种情况下都有效的实现。实现这一目标的标准方法是什么？

【问题讨论】：

你能有一个模板方法 IndexIsEmpty() 对于非稀疏存储总是返回 true 吗？然后两次迭代都会调用 IndexIsEmpty()。希望优化器会丢弃所有非稀疏冗余检查。

惯用的 C++ 方法将定义 dense_iterator 和 sparse_iterator 并将操作编写为模板函数

@Caleth，类似的东西正是我想要的。您对迭代器有很好的参考吗？我发现的大多数来源看起来都非常令人生畏。

【参考方案1】：

我也遇到了同样的问题。由于我懒得写迭代器（这很快就会遇到限制），我决定使用 lambda 方法，为每个矩阵定义专门的函数，每个元素调用一个 lambda：

template<class Func>
void forMatrixEntries(const VecOfVecMatrix& mat, Func&& func)

    for (auto& row : mat.getData())
        for (auto& elem : row)
            func(elem); // You could track and pass the indices if you want.


template<class Func>
void forMatrixEntries(const CompressedSparseMatrix& mat, Func&& func)

    for (auto& elem : mat.getElements())
        func(elem); // You could track and pass the indices if you want.

（这些也可以是成员函数，因此它们可以更轻松地访问内部 - 您的选择）。 然后您可以统一使用它们：

template<class Mat>
void scale(const Mat& mat, double factor)

    forMatrixEntries(mat, [factor](double& elem) 
        elem *= factor;
    );

唯一的缺点是矩阵专用函数当然需要在标题中（因为模板）。但我认为这种方法不仅优雅而且非常有表现力（你可以给“迭代矩阵条目”命名而不是复杂的循环语法，但循环体保持不变）。

【讨论】：

谢谢，很优雅！我忘记了 lambda 函数是 C++ 中的一个东西。矩阵专用函数仍然可以放在 cpp 文件中，但是您需要显式地实例化它们。 :)

scale-like 函数是的，那些采用Func 的函数没有。不确定你指的是哪一个。