在数值模拟中使用 std::valarray

Question

我在 Code Review 中发布了一个简单的 n-body class，它是我用 C++ here 编写的。

我被告知要使用 std::valarray 而不是普通的 std::array ，目的是我可以重写一些目前看起来像这样的代码

void Particle::update_position() {
    for (unsigned int d = 0; d < DIM; ++d) {
        position[d] += dt*(velocity[d] + a*force_new[d]);
        force_old[d] = force_new[d];
    }
}

至此

void Particle::update_position() {
    position += 0.1*(velocity + force_new);
    force_old = force_new;
}

因为我是 C++ 的初学者，所以我写了一个使用 std::valarray 的小程序，这样我就可以学习如何使用这个数据结构。但是，编译器会抛出很多错误，我不知道为什么。我希望你能帮我解决这个问题。下面是我写的小程序：

#include <iostream>
#include <vector>
#include <array>
#include <valarray>

constexpr unsigned int DIM = 2;

struct Particle{
    std::array<std::valarray<double>, DIM> position; 
    std::array<std::valarray<double>, DIM> velocity;
    std::array<std::valarray<double>, DIM> force_new;
    std::array<std::valarray<double>, DIM> force_old;
    void update_position();
};

void Particle::update_position() {
    position += 0.1*(velocity + force_new);
    force_old = force_new;
}

void compute_position(std::vector<Particle>& particles) {
    for (auto& particle: particles) { 
        particle.update_position();
    }
}

void print_data(const std::vector<Particle>& particles) {
    for (const auto& particle: particles) {
        for (const auto& x: particle.position) std::cout << x << " ";
        for (const auto& v: particle.position) std::cout << v << " ";
        for (const auto& F: particle.position) std::cout << F << " ";
        std::cout << std::endl;
    }
}

void init_data(std::vector<Particle>& particles) {
    for (auto& particle: particles) {
        for (const auto& p: particle) {
            p.position = 1.0
            p.velocity = 2.0
            p.force_new = 3.0
            p.force_old = 4.0
        }
    }
}

int main() { 
    const unsigned int n = 10;
    std::vector<Particle> particles(n);
    init_data(particles);
    compute_position(particles);
    print_data(particles); 
    return 0;
}

当我尝试编译这段代码时，出现以下错误：

so.cpp: In member function ‘void Particle::update_position()’:
so.cpp:17:31: error: no match for ‘operator+’ (operand types are ‘std::array<std::valarray<double>, 2>’ and ‘std::array<std::valarray<double>, 2>’)
     position += 0.1*(velocity + force_new);

so.cpp: In function ‘void print_data(const std::vector<Particle>&)’:
so.cpp:29:58: error: no match for ‘operator<<’ (operand types are ‘std::ostream’ {aka ‘std::basic_ostream<char>’} and ‘const std::valarray<double>’)
         for (const auto& x: particle.position) std::cout << x << " ";


so.cpp: In function ‘void init_data(std::vector<Particle>&)’:
so.cpp:38:29: error: ‘begin’ was not declared in this scope
         for (const auto& p: particle) {
                             ^~~~~~~~
so.cpp:38:29: note: suggested alternative:
In file included from so.cpp:4:
/usr/include/c++/8/valarray:1211:5: note:   ‘std::begin’
     begin(const valarray<_Tp>& __va)
     ^~~~~
so.cpp:38:29: error: ‘end’ was not declared in this scope
         for (const auto& p: particle) {

Answer 1

您的 Particle 成员是 std::array 的 std::valarray 的 double 的 ，而不仅仅是一个简单的 std::valarray.也就是说，标准库只提供了operator+, operator*, [...] for the std::valarray，其余的你需要自己提供。

例如，以下将解决行 position += 0.1*(velocity + force_new); (See live online)

中的第一个编译器错误

template<typename ValArray, std::size_t N>
auto operator+(std::array<ValArray, N> lhs, const std::array<ValArray, N>& rhs)
{
   for (std::size_t idx{}; idx < N; ++idx)
      lhs[idx] += rhs[idx];
   return lhs;
}


template<typename ValArray, std::size_t N, typename T>
auto operator*(T constant, std::array<ValArray, N> rhs)
{
   for (std::size_t idx{}; idx < N; ++idx)
      rhs[idx] *= constant;
   return rhs;
}

或者将 std::array<std::valarray<double>, DIM> 包装成 class 并提供必要的运算符，继续进行。 (See live online)

#include <array>
#include <valarray>    

template<typename Type, std::size_t N>
class ValArray2D final
{
   std::valarray<std::valarray<Type>> _arra2D;
public:
   explicit constexpr ValArray2D(const std::valarray<Type>& valArray) noexcept
      : _arra2D{ valArray , N }
   {}

   ValArray2D& operator+=(ValArray2D rhs) noexcept
   {
      _arra2D += rhs._arra2D;
      return *this;
   }

   ValArray2D& operator*(Type constant) noexcept
   {
      for(auto& valArr: this->_arra2D)
         valArr = valArr * constant;
      return *this;
   }

   void operator=(Type constant) noexcept
   {
      for (auto& valArr : this->_arra2D)
         valArr = constant;
   }

   friend ValArray2D operator+(ValArray2D lhs, const ValArray2D& rhs) noexcept
   {
      lhs += rhs;
      return lhs;
   }

   friend std::ostream& operator<<(std::ostream& out, const ValArray2D& obj) noexcept
   {
      for (const std::valarray<Type>& valArray : obj._arra2D)
         for (const Type element : valArray)
            out << element << " ";
      return out;
   }
};

template<typename Type, std::size_t N>
ValArray2D<Type, N> operator*(Type constant, ValArray2D<Type, N> lhs)
{
   return lhs = lhs * constant;
}

Answer 2

首先，当您编写或更改代码时，始终从第一个工作版本开始，并确保代码在每个步骤之间进行编译。这将使隔离错误编译代码变得更加容易。

我为什么要告诉你这个？这是因为您的部分代码从未正确编译过。无论是在引入valarray之前还是之后。

例如，这个：

for (auto& particle : particles) {
    for (const auto& p: particle) {
        p.position = 1.0
        p.velocity = 2.0
        p.force_new = 3.0
        p.force_old = 4.0
    }
}

单粒子不是可迭代类型，行尾没有分号

一步一步来，确保代码在每个步骤之间编译。

其次，我认为 valarray 不是您要找的东西。除非你希望每个粒子都具有每个属性的动态维数，否则这将是非常令人惊讶的。

我建议您引入一个 vec 类型，它包含您进行简化所需的组件。这样的 vec 类型可以在库中找到，例如 glm 提供了一个 vec2 class，它具有 +-/* 等运算符。

即使没有库，也可以创建简单的向量类型。

这是一个使用向量（在数学意义上）而不是 std::valarray 的代码示例：

struct Particle{
    glm::vec2 position; 
    glm::vec2 velocity;
    glm::vec2 force_new;
    glm::vec2 force_old;
    void update_position();
};

void Particle::update_position() {
    position += 0.1*(velocity + force_new);
    force_old = force_new;
}

void compute_position(std::vector<Particle>& particles) {
    for (auto& particle: particles) { 
        particle.update_position();
    }
}

void print_data(const std::vector<Particle>& particles) {
    for (const auto& particle : particles) {
        std::cout << particle.position.x << ", " << particle.position.y << " ";
        std::cout << particle.velocity.x << ", " << particle.velocity.y << " ";
        std::cout << particle.force_new.x << ", " << particle.force_new.y << " ";
        std::cout << std::endl;
    }
}

void init_data(std::vector<Particle>& particles) {
    for (auto& particle : particles) {
        particle.position = {1, 2};
        particle.velocity = {2, 24};
        particle.force_old = {1, 5};
        particle.force_new = {-4, 2};
    }
}

int main() { 
    const unsigned int n = 10;
    std::vector<Particle> particles(n);
    init_data(particles);
    compute_position(particles);
    print_data(particles); 
    return 0;
}

Live example with custom (imcomplete) vec2 type

Answer 3

好的，我们来解释一些概念，以便您了解错误的原因。

std::valarray vs 普通 C++ 数组 vs std::array

当我们谈论纯 C++ 数组时，我们指的是类似下面的内容

double myArray[DIM];

这只是 C++ 中数组的原生形式。然后我们有 std::array

std::array<double, DIM> myStdArry;

这只是一个简单的 c++ 数组的包装器 class，不同之处在于您可以访问一些实用函数以便更轻松地操作数据，即

myStdArry.size() //get the number of elements
myStdArry.begin() & myStdArry.end() //iterators to the begin and end of the array
myStdArry.fill(1.4) //replace all the values of the array to 1.4

对于普通数组和标准数组，您必须使用下标运算符 ([]) 才能访问数组的每个元素，即

for (size_t i = 0; i < DIM /*or myStdArry.size()*/; ++i;) {
  myArray[i] = ...;
} 

//Or with range-based for loop

for (auto& element : myArray) {
  element = ...;
}

因此，您不能直接在容器（数组）上使用算术运算符，因为它们没有超载。这是 valarray 出现的时候，因为它是专门为这种操作而设计的。 Valarray 只是另一种类型的数组，它通过将算术运算符应用于数组的每个元素来重载它们。

即假设我们想要对数组的每个元素进行平方。使用普通数组和标准数组，我们可以通过以下方式实现它：

for (auto& element : myStdArray) {
  element *= element;
}

但是使用 valarray 我们可以简单地做：

myValArray *= myValArray;

我们得到了相同的结果。

其他主要区别是，虽然普通数组和标准数组都是固定大小（您必须在编译时设置其大小），但 val 数组的大小可以动态增长，这就是您不必指定其大小的原因在申报时但直到施工或之后

myValArray{DIM} //Or
myValArray.resize(DIM)

在数值模拟中使用 std::valarray

Using std::valarray in numerical simulation

c++

c++11

data-structures

stdarray

valarray

std::valarray vs 普通 C++ 数组 vs std::array