使用重载运算符的OpenMP归约

Question

我正在尝试使用OpenMP并行化下面函数中的循环。

 void CEnergymulti::forcetwobody(vector<CMolecule*>  m_mols,CPnt force0,CPnt torque0)
{
 const int nmol=m_mols.size();
 vector<CMolecule*> twomols(2);
 CPnt forcetemp,torquetemp;
 twomols.clear();
 force0.zero();
 torque0.zero();
 forcetemp.zero();
 torquetemp.zero();
 #pragma omp parallel for reduction(+:force0,torque0) private(twomols)
 for(int j=1;j<nmol;j++)
       { twomols.push_back(m_mols[0]);
         twomols.push_back(m_mols[j]);
         CMolecule::polarize_mutual(twomols,false, 1000);
         twomols[0]->computeMol_Force_and_Torque(forcetemp,torquetemp);
         force0+=forcetemp;
         torque0+=torquetemp;
         forcetemp.zero();
         torquetemp.zero();
         twomols.clear();
        }
     REAL converter=COUL_K*IKbT;
     force0*=converter;
     torque0*=converter;
     return;
     }

当我编译代码时，会出现以下提示信息：

EnergyD_multi.cpp: In static member function ‘static void
CEnergymulti::forcetwobody(std::vector<CMolecule*,
std::allocator<CMolecule*> >, CPnt, CPnt)’: EnergyD_multi.cpp:226:
error: ‘torque0’ has invalid type for ‘reduction’
EnergyD_multi.cpp:226: error: ‘force0’ has invalid type for
‘reduction’

我知道变量'force0'和'torque0'不是double或integer类型的数据，而是'CPnt'类型。'CPnt'是一个表示三维向量的类。对于类'CPnt'，操作符'+'和'-'已经被重载定义了。因此，我的问题是：OpenMP中的reduction能否处理这样的重载操作符？有没有其他方法在OpenMP中并行化这个循环，而不是对'force0'和'torque0'的每个分量都进行reduction？

非常感谢。

- user2226358

---

你搞定了吗？我已经修改了你的代码，使用OpenMP进行了重载操作。如果可以，请告诉我它是否有效。如果有错误，请查看链接以了解主要思路。 - user2088790

@raxman。实际上，即使我根据您的修改编译了代码，当我运行代码时，仍然会出现导致核心转储的分段错误。但是，如果我将twomols和forcetemp添加为私有变量，代码就可以平稳运行，尽管结果仍与串行代码不同。我的代码中可能有某个函数调用有问题，但是您认为在使用“nowait”时是否有必要拥有“private”变量？谢谢。 - user2226358

如果没有完整的代码和了解你在做什么，我很难调试你的代码。我最初发布了我的工作示例，以便你可以理解这个概念。你应该能够从链接中获得主要思路。你应该得到与串行代码相同的结果。 - user2088790

不等待只是移除了for循环的屏障，因此当线程完成时，它会直接进入关键块。 - user2088790

你确定你正确使用了 "vector<CMolecule*>" 吗？这应该是一个指针的向量，但是看起来你在推入值时使用了 "twomols.push_back(m_mols[0]);"。也许你想要使用 &m_mols[0] 或者 vector<CMolecule>？ - user2088790

Answer 1

OpenMP的归约无法处理这样的重载运算符，但是有一种替代方法。在OpenMP中重写归约的一种方式是使用nowait和atomic参数。http://bisqwit.iki.fi/story/howto/openmp/#ReductionClause。这种方法与正常方式一样快。

如果将atomic替换为critical，则可以使用更复杂的重载运算符。根据我的经验，虽然不如使用atomic快，但仍然表现良好。

我这样做是为了使用能够同时操作4个或8个浮点数（使用SEE或AVX）的运算符。reduction with OpenMP with SSE/AVX 编辑：我更改了您的代码以反映我认为会实现您想要的功能的内容。

void CEnergymulti::forcetwobody(vector<CMolecule*>  m_mols,CPnt force0,CPnt torque0)
{
    const int nmol=m_mols.size();
    force0.zero();
    torque0.zero();
    #pragma omp parallel
    {
        CPnt force0_private;
        CPnt torque0_private; 
        force0_private.clear();
        torque0_private.clear();
        #pragma omp for nowait
        for(int j=1;j<nmol;j++)
        { 
            CPnt forcetemp,torquetemp;
            forcetemp.zero();
            torquetemp.zero();
            vector<CMolecule*> twomols(2);
            twomols.clear();
            twomols.push_back(m_mols[0]);
            twomols.push_back(m_mols[j]);
            CMolecule::polarize_mutual(twomols,false, 1000);
            twomols[0]->computeMol_Force_and_Torque(forcetemp,torquetemp);
            force0_private+=forcetemp;
            torque0_private+=torquetemp;
        }
        #pragma omp critical 
        {
           force0 += force0_private;
           torque0 += torque0_private;
        }

    }
    REAL converter=COUL_K*IKbT;
    force0*=converter;
    torque0*=converter;
    return;
}

- user2088790

---

我认为第二个编译指示（即 #pragma omp parallel for nowait）应该改成 #pragma omp for nowait，即去掉 parallel（parallel 构造已在外部设置）。我的编译器报错，说不能同时使用 nowait 和 parallel for，因为在 parallel for 之后应该有一个隐式屏障。似乎 nowait 和 parallel for 是冲突的指令。 - Varun Gulshan

@VarunGulshan，你是正确的，应该是#pragma omp for nowait。我已经修复了它。 - Z boson

这样做不会失去对数复杂度吗？

最坏情况是每个线程都在等待进入临界区，然后它们按顺序添加其结果。

因此，最终您将得到 O(nmol / thread_num) + O(thread_num) 而不是 O(nmol / thread_num) + O( log (thread_num))

这可能会在强大的计算机（例如超级计算机）上产生巨大的差异。 - Jan

Answer 2

现在，你也可以使用declare reduction指令（OpenMP 4.0+）来实现；例如：

#pragma omp declare reduction(mysum:CPnt:omp_out += omp_in) initializer(omp_priv.zero())
#pragma omp parallel for reduction(mysum:force0,torque0) private(twomols)
for(int j=1;j<nmol;j++)
...

- Streets