专门用于访客模式的循环

specialise for loop with visitor pattern?

本文关键字:循环 模式 访客 用于      更新时间:2023-10-16

假设我有一个结构,名称如下:

struct OneBodyMatrixElements {
  static const bool is_selfadjoint = 1;
  static const bool property2 = 0; 
  // ........
  static value_type val( const std::size_t i, const std::size_t j ) {
    value_type ret = 0.;
    if( i == j )  ret = 1.;
    if( i == j + 1 || i == j - 1 ) ret = -1.;
    if( i == 0 && j == 6 ) ret = -1.;
    if( i == 6 && j == 0 ) ret = -1.;
    return static_cast<value_type>(ret);
  }
};

静态变量定义了这个类的不同属性,静态函数val根据(i,j)返回值。这实际上表示一个大小为7的矩阵,其标量类型为value_type(float、double、std::complex)。

在我的代码的另一部分中,我需要循环(I,j),因为I,j=0,。。。,N(例如N=7):

for( std::size_t i = 0; i < N; i++ ) {
  for( std::size_t j = 0; j < N; j++ ) {
    if( "OneBodyMatrixElements::val( i, j ) is non zero" ) {
      //...
      OneBodyMatrixElements::val( i, j ) * func( i, j );
      //...
    }
}

我有一个if close来检查OneBodyMatrixElements是否为零,因为下面的乘法是昂贵的,并且如果前置因子为零,则不需要计算。

我的问题是:既然我知道OneBodyMatrixElements的性质,即自伴或对角线,我也知道零元素的位置,我如何优化两个for循环,使其仅在非零元素上迭代,并考虑自伴性质?

这个问题的明显答案是使用特征稀疏矩阵,但我不是我想要的,因为

  • 我不想实际存储矩阵
  • 实际上,我有另一个具有四个索引(I,j,k,l)的结构,并在这四个索引上循环。因此,我需要一些通用的东西,也适用于四个指数。在Eigen中,不支持稀疏张量

我确实看过访客模式,但我不确定这是否是需要的。对代码的任何注释也是受欢迎的,并且可以接受以稀疏形式将值实际存储在内存中的解决方案。只涉及非零元素上的迭代的解决方案也是可以接受的。非常感谢。

编辑:答案:

非常感谢你的回答。我确实实现了一个自定义迭代器,如下所述。为了完整起见,这里是我的版本:

class tridiag_iterator {
  std::size_t i, j;
public:
  struct value_type {
    std::size_t i, j;
  };
  tridiag_iterator( std::size_t i_, std::size_t j_ ) : i( i_ ), j( j_ ) { }
  value_type operator*() const {
    return { i, j };
  }
  bool operator!=( const tridiag_iterator &other ) {
    if( i >= other.i && j > other.j ) return 0;
    else return 1;
  }
  tridiag_iterator &operator++() {
    if( j == i ) j = i + 1;
    else if( j == i - 1 ) j = i;
    else if( j == i + 1 ) {
      i = i + 1;
      j = i - 1;
    }
  }
};

我在OneBodyMatrixElements中添加了两个成员:

  static tridiag_iterator begin() {
    return tridiag_iterator( 0, 0 );
  }
  static tridiag_iterator end() {
    return tridiag_iterator( N - 1, N - 1 );
  }

其中N在别处定义。用法:

  for( auto it = OneBodyMatrixElements::begin(), end = OneBodyMatrixElements::end(); it != end; ++it ) {
    OneBodyMatrixElements::val( (*it).i, (*it).j ) * func( (*it).i, (*it).j );
  }

欢迎任何评论。

您可以让OneBodyMatrixElements在非零元素上提供迭代器。类似这样的东西:

struct Iterator
{
  struct value_type
  {
    std::size_t i, j;
    OneBodyMatrixElements::value_type val;
  };
  Iterator() : i(7), j(7) {}
  Iterator(std::size_t i, std::size_t j) : i(i), j(j) {}
  value_type operator* () const
  {
    return {i, j, OneBodyMatrixElements::val(i, j)};
  }
  Iterator& operator++ ()
  {
    if (j == i || j == i - 1) j = i + 1;
    else if (j == i + 1) { i = i +1; j = i - 1; }
    else :::
  }
private:
  std::size_t i, j;
};

struct OneBodyMatrixElements
{
  static Iterator begin()
  {
    return Iterator(0, 0);
  }
  static Iterator end()
  {
    return Iterator();
  }
};
// Usage:
for (auto it = OneBodyMatrixElements::begin(), end = OneBodyMatrixElements::end(); it != end; ++it)
{
  it->val * func(it->i, it->j);
}

上面的代码旨在说明这个想法,而不是提供一个100%现成的解决方案。必须添加其他运算符定义、迭代器的typedef等。

我的第一种方法是简单地让编译器找出哪些是零。由于生成器函数相对简单(见下文),不依赖于其自身之外的任何状态,因此对其进行足够深入的分析以内联它并消除不必要的迭代应该是可行的。

我的第二种方法是提供一个游标/迭代器,它只返回非零元素的位置,可能还会返回它们的值。在任何情况下,都应该首先进行分析!

无论如何,编译可能需要一些帮助:

  • static_cast的结果ret,已经有了该类型,删除它。相反,您可能必须强制转换if子句中的赋值
  • 如果所有可能的类型都可以在布尔上下文中使用,请使用组合赋值和检查来避免额外的调用:if (auto v = val(i, j)) ...
  • 使用else来明确第一个和匹配的第一个if子句决定结果。或者,使用return而不是分配给一个临时的来明确这一点
  • 关于你的目标,另一件需要考虑的事情是不要存储矩阵。如果您只拼写或计算了一次,您将只有几个索引操作,这些操作通常可以从缓存中传递。也许内存中的矩阵甚至比计算它的代码还要小