处理密集矩阵和稀疏矩阵

Question

处理密集矩阵和稀疏矩阵

3

我正在编写一个C++函数，该函数操作作为参数传递的矩阵，并希望代码能够处理不同类型的矩阵（例如Boost稀疏矩阵、std::vectors of 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];
    }
}

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

- Rastapopoulos

你能否创建一个名为IndexIsEmpty()的模板方法，对于非稀疏存储，它始终返回true？这样两个迭代都将调用IndexIsEmpty()。希望优化器会丢弃所有非稀疏冗余检查。 - Gem Taylor

1

惯用的C++方法是定义一个dense_iterator和一个sparse_iterator，并将操作编写为模板函数。 - Caleth

@Caleth，这正是我正在寻找的东西。你有关于迭代器的好参考资料吗？我发现大多数来源看起来都很令人生畏。 - Rastapopoulos

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Max Langhof · Accepted Answer

我曾经遇到过同样的问题。因为我懒得编写迭代器（那会很快遇到限制），所以我决定采用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;
    });
}

唯一的缺点是矩阵专用函数必须在头文件中（因为模板）。但我认为这种方法不仅优雅而且非常表达（您可以给“迭代矩阵条目”的名称，而不是复杂的循环语法，但循环体保持不变）。