version with functor

main.cpp

@@ -3,6 +3,7 @@
 #include <cstring>
 #include "sound_column.h"
 #include <execution>
+#include <algorithm>
 
 #define DEFAULT_DISTANCE    20'000.0f
 #define DEFAULT_DEPTH       5'000.0f
@@ -67,15 +68,15 @@ int run_simulation(float r_max, float z_max, float f_min, float f_max){
     auto nf = size_t(f_max / f_min);
     float delta_r = r_max/10'000;
 
+    std::vector<sound_column*> column(nf, nullptr);
 
-    std::vector<sound_column*> column;
-    column.reserve(nf);
+    //std::ranges::iota_view indexes((size_t) 0, nf);
 
-    std::ranges::iota_view indexes((size_t) 0, nf);
-    std::for_each(std::execution::par_unseq, indexes.begin(), indexes.end(),
-                  [&](size_t i) {column[i] = new sound_column(delta_r, (size_t) 1.10 * z_max*(4*(i+1)*f_min)/WATER_SPEED, 7, r_max, z_max, (i+1)*f_min, WATER_SPEED);
-                      column[i]->init_source(gaussian_init);});
-    std::for_each(std::execution::par_unseq, indexes.begin(), indexes.end(),
-                  [&](size_t i) {while (column[i]->step_r(media));});
+    std::for_each(std::execution::par_unseq, column.begin(), column.end(),
+                  [&](sound_column*& c) {int i = &c - &column[0];
+                  c = new sound_column(delta_r, (size_t) 1.10 * z_max*(4*(i+1)*f_min)/WATER_SPEED, 7, r_max, z_max, (i+1)*f_min, WATER_SPEED);
+                      c->init_source(gaussian_init);});
+    std::for_each(std::execution::par_unseq, column.begin(), column.end(),
+                  [&](sound_column*& c) {while (c->step_r(media));});
     return 0;
 }

sound_column.cpp

@@ -14,13 +14,13 @@ sound_column::sound_column(float dr, size_t nz, size_t np, float dm, float zm, f
     n_p=np;
     freq=f;
     c0=speed;
-    pressures.reserve(nz);
-    _u_.reserve(nz-1);
-    _l_.reserve(nz-1);
-    _d_.reserve(nz);
-    n_next.reserve(nz);
-    fapbj.reserve(np);
-    fbmaj.reserve(np);
+    pressures.resize(nz);
+    _u_.resize(nz-1);
+    _l_.resize(nz-1);
+    _d_.resize(nz);
+    n_next.resize(nz);
+    fapbj.resize(np);
+    fbmaj.resize(np);
     for(int i=0; i < np; i++) {
         faj = std::complex<float>(0, 4* M_1_PIf32 * freq / c0)  / float((2.f * n_p + 1.f) *
                                     powf(sinf(float(i+1) * M_1_PIf32 / (2.f * n_p + 1.f)), 2) * delta_r);
@@ -36,34 +36,35 @@ sound_column::~sound_column() = default;
 bool sound_column::step_r(float (*fun)(float, float)) {
     if (r >= d_max)
         return true;
-    std::ranges::iota_view z_ind((size_t) 0, n_z);
+    //std::ranges::iota_view z_ind((size_t) 0, n_z);
     r += delta_r;
-    std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
-                   [&](size_t i) {n_next[i] = (*fun)(r, z_max/n_z*i);});
+    std::for_each(std::execution::par_unseq, n_next.begin(), n_next.end(),
+                   [&](float& n) {int i = &n - &n_next[0]; n = (*fun)(r, z_max/n_z*i);});
     for(int i=0; i < n_p;i++){
         _step_r_(i);
     }
     // Hank
-    std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
-                  [&](size_t i) {pressures[i] *= 2.f / (M_1_PIf32 * M_1_PIf32 * freq / c0 * r) *
+    std::for_each(std::execution::par_unseq, pressures.begin(), pressures.end(),
+                  [&](std::complex<float>& p) {p *=  static_cast<std::complex<float>>(2.f / (M_1_PIf32 * M_1_PIf32 * freq / c0 * r)) *
                   std::exp(std::complex<float>(0, 2.f* M_1_PIf32 * freq / c0 * r - M_1_PIf32 / 4));});
     return false;
 }
 
 void sound_column::init_source(std::complex<float> (*fun)(float, float)) {
-    std::ranges::iota_view z_ind((size_t) 0, n_z);
-    std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
-                  [&](size_t i) {pressures[i] = (*fun)(z_max/n_z*i, 2*M_1_PIf32 * freq / c0);});
+//    std::ranges::iota_view z_ind((size_t) 0, n_z);
+    std::for_each(std::execution::par_unseq, pressures.begin(), pressures.end(),
+                  [&](std::complex<float>& p) {int i = &p - &pressures[0];
+                  p = (*fun)(z_max/n_z*i, 2*M_1_PIf32 * freq / c0);});
 }
 
 void sound_column::_step_r_(int ip) {
-    std::ranges::iota_view z_ind((size_t) 0, n_z);
+//    std::ranges::iota_view z_ind((size_t) 0, n_z);
     //tri diag init
     std::fill(std::execution::par_unseq, _u_.begin(), _u_.end(), fbmaj[ip]/h2);
     std::fill(std::execution::par_unseq, _l_.begin(), _l_.end(), fbmaj[ip]/h2);
 
-    std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
-                  [&](size_t i) {_d_[i] = 4* M_1_PIf32 * freq / c0 + fbmaj[ip] /
+    std::for_each(std::execution::par_unseq, _d_.begin(), _d_.end(),
+                  [&](std::complex<float>& d) {int i = &d - &_d_[0];d =  static_cast<std::complex<float>>(4* M_1_PIf32 * freq / c0) + fbmaj[ip] /
                   static_cast<std::complex<float>>(-2.f/h2 + powf(4* M_1_PIf32 * freq / c0, 2) * (powf(n_next[i], 2) - 1.f));});
 
     //  first matmul

sound_column.h

@@ -6,8 +6,9 @@
 #define PARABOLIC___SOUND_COLUMN_H
 #include <vector>
 #include <complex>
-#include <ranges>
+//#include <ranges>
 #include <execution>
+#include <algorithm>
 
 # define M_1_PIf32	(0.318309886183790671537767526745028724)
 
@@ -30,6 +31,7 @@ private:
     std::vector<std::complex<float>> _l_;
     std::vector<std::complex<float>> _d_;
     std::vector<float> n_next;
+    std::complex<float> (*fun)(float, float) source_fun;
     void _step_r_(int ip);
 
 
@@ -38,6 +40,7 @@ public:
     ~sound_column();
     bool step_r(float (*fun)(float, float));
     void init_source(std::complex<float> (*fun)(float, float));
+    void init_source();
 
 };