小编Ale*_*exD的帖子

我正在尝试使用 OpenMP 缩减来并行化以下循环；

#define EIGEN_DONT_PARALLELIZE
#include <iostream>
#include <cmath>
#include <string>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Eigenvalues>
#include <omp.h>

using namespace Eigen;
using namespace std;

VectorXd integrand(double E)
{
    VectorXd answer(500000);
    double f = 5.*E + 32.*E*E*E*E;
    for (int j = 0; j !=50; j++)
        answer[j] =j*f;
    return answer;
}

int main()
{
    omp_set_num_threads(4);
    double start = 0.;
    double end = 1.;
    int n = 100;
    double h = (end - start)/(2.*n);

    VectorXd result(500000);
    result.fill(0.);
    double E = start;
    result = integrand(E); …

我正在编写一些计算成本高但高度可并行化的代码。并行化后，我打算在 HPC 上运行它，但是为了将运行时间缩短到一周内，问题需要随着处理器数量的增加而很好地扩展。

下面是我试图实现的一个简单而荒谬的例子，它足够简洁，可以编译和演示我的问题；

#include <iostream>
#include <ctime>
#include "mpi.h"

using namespace std;

double int_theta(double E){
    double result = 0;
    for (int k = 0; k < 20000; k++)
        result += E*k;
    return result;
}

int main() 
{
    int n = 3500000;
    int counter = 0;
    time_t timer;
    int start_time = time(&timer);
    int myid, numprocs;
    int k;
    double integrate, result;
    double end = 0.5;
    double start = -2.;
    double E;
    double factor = (end - start)/(n*1.);
    integrate = 0;
    MPI_Init(NULL,NULL); …

我有一个函数，它要求返回类型是一个容器。问题是我需要尽可能有效地集成容器的内容，并希望使用自适应 Gauss-Kronrod 集成或同样有效（或更好）的方法。

我希望使用 GNU Scientific Library 正交例程 qags，但它返回的结果是 double 类型。就目前情况而言，我认为我可能不得不重写 GSL 正交例程的部分以将返回类型转换为 std 向量，但这将是一个相当冗长且可能容易出错的弯路。我希望有人可以推荐更好的解决方案！

下面是我试图整合的那种函数的例子，目前使用基本的梯形规则，但表明我更喜欢在哪里实现 GSL 例程 gaq。虽然它比我的实际问题简单得多，但它表明放入向量中的每个元素都是根据先前的结果计算的，因此需要容器。

#include <iostream>
#include <vector>
#include <gsl/gsl_integration.h>
using namespace std;

vector<double> f(double E, int N) {
    vector<double> result;
    result.reserve(N);
    double x = E;
    for (int it=0; it < N; ++it){
        result.push_back(x);
        x *= x;
    }
    return result;
}

vector<double> f_qag(double E, void * params) {
    int N = *(int *) params;
    vector<double> result;
    result.reserve(N);
    double x = E;
    for (int it=0; it …