Added C++ lang implementation of the split operator

[rajdakin]

This pull request is one of the two new pull requests requested in #409 and is made to integrate the C++ lang with the split operator method.

[leios]

A few comments:

1. The OOP nature of this makes it difficult to read. I would argue that most people implementing this algorithm will be coming from a Scientific computing background, meaning that this style of programming might not comfortable for them. I am up for debate here, though.
2. The energy calculation can also be put into the quantum systems chapter. Would you be willing to put your energy calculation code there as well?

C++ is a sticky point for the Algorithm Archive in that a lot of people tend to want to write different styles of C++ code. We tend to try to take the middle ground with this project and avoid overly OOP and Functional programming and instead use what makes sense at the time for readability.

[berquist]

+1 for not using classes when you don't have any methods

[nic-hartley]

A couple of big things:

• Please, please use std::. I know it's not required on most compilers, but either it or using is by the standard, and using is a bad habit to give peope.
• If you're going to use C-style string formatting and IO, do it right, especially in teaching code.

[rajdakin]

@nic-hartley Where do I implicitly use std:: instead of explicitly use it?

[nic-hartley]

@rajdakin Basically everywhere you were using a C function, before the changes. The standard mandates that they be put in std, and allows them to be put in the global namespace as well, but it explicitly says it's undefined (not even implementation-defined) whether or not they are.

Or, in the standard's words:

S20.5.1.1.2-3

All library entities except operator new and operator delete are defined within the namespace std or namespaces nested within namespace std. It is unspecified whether names declared in a specific namespace are declared directly in that namespace or in an inline namespace inside that namespace. Whenever a name x defined in the standard library is mentioned, the name x is assumed to be fully qualified as ::std::x, unless explicitly described otherwise. For example, if the Effects: section for library function F is described as calling library function G, the function ::std::G is meant.

S20.1.5.2.4

Except as noted in Clauses 20 through 33 and Annex D, the contents of each header cname is the same as that of the corresponding header name.h as specified in the C standard library (Clause 2). In the C++ standard library, however, the declarations (except for names which are defined as macros in C) are within namespace scope (6.3.6) of the namespace std. It is unspecified whether these names (including any overloads added in Clauses 21 through 33 and Annex D) are first declared within the global namespace scope and are then injected into namespace std by explicit using-declarations (10.3.3).

[ Note: The names defined as macros in C include the following: assert, offsetof, setjmp, va_arg, va_end, and va_start. — end note ]

(formatting omitted for ease of copy/pasting)

latest draft, p. 456-457

Just to be explicitly clear, though: You are not implicitly using it anymore.

[nic-hartley]

Looks good to me now!

[berquist]

I am getting compilation errors:

$ g++ -lfftw3 split_op.cpp
split_op.cpp: In constructor ‘Operators::Operators(Params&, double, double)’:
split_op.cpp:61:84: error: missing template arguments before ‘(’ token
                 ke.emplace_back(exp(-0.5 * par.dt * pow(par.k[i], 2) * std::complex(0.0, 1.0)));
                                                                                    ^
split_op.cpp:62:72: error: missing template arguments before ‘(’ token
                 pe.emplace_back(exp(-0.5 * par.dt * v[i] * std::complex(0.0, 1.0)));
                                                                        ^
split_op.cpp: In function ‘double calculate_energy(Params, Operators)’:
split_op.cpp:150:58: error: missing template arguments before ‘(’ token
     std::vector<std::complex<double>> wfc_r = std::vector(opr.wfc);
                                                          ^
split_op.cpp:151:58: error: missing template arguments before ‘(’ token
     std::vector<std::complex<double>> wfc_k = std::vector(opr.wfc);
                                                          ^
split_op.cpp:163:50: error: missing template arguments before ‘(’ token
         energy_k[i] = wfc_k[i] * pow(std::complex(par.k[i], 0.0), 2);
                                                  ^

[berquist]

diff --git a/contents/split-operator_method/code/c++/split_op.cpp b/contents/split-operator_method/code/c++/split_op.cpp
index 814ef2a..525d66f 100644
--- a/contents/split-operator_method/code/c++/split_op.cpp
+++ b/contents/split-operator_method/code/c++/split_op.cpp
@@ -58,8 +58,8 @@ public:
                 ke.emplace_back(exp(-0.5 * par.dt * pow(par.k[i], 2)));
                 pe.emplace_back(exp(-0.5 * par.dt * v[i]));
             } else {
-                ke.emplace_back(exp(-0.5 * par.dt * pow(par.k[i], 2) * std::complex(0.0, 1.0)));
-                pe.emplace_back(exp(-0.5 * par.dt * v[i] * std::complex(0.0, 1.0)));
+                ke.emplace_back(exp(-0.5 * par.dt * pow(par.k[i], 2) * std::complex<double>(0.0, 1.0)));
+                pe.emplace_back(exp(-0.5 * par.dt * v[i] * std::complex<double>(0.0, 1.0)));
             }
         }
     }
@@ -147,8 +147,8 @@ void split_op(Params &par, Operators &opr) {
 }
 
 double calculate_energy(Params par, Operators opr) {
-    std::vector<std::complex<double>> wfc_r = std::vector(opr.wfc);
-    std::vector<std::complex<double>> wfc_k = std::vector(opr.wfc);
+  std::vector<std::complex<double>> wfc_r(opr.wfc);
+  std::vector<std::complex<double>> wfc_k(opr.wfc);
     std::vector<std::complex<double>> wfc_c;
     fft(wfc_k, opr.size, false);
 
@@ -160,7 +160,7 @@ double calculate_energy(Params par, Operators opr) {
     std::vector<std::complex<double>> energy_r;
 
     for (size_t i = 0; i < opr.size; ++i) {
-        energy_k[i] = wfc_k[i] * pow(std::complex(par.k[i], 0.0), 2);
+        energy_k[i] = wfc_k[i] * pow(std::complex<double>(par.k[i], 0.0), 2);
     }
 
     fft(energy_k, opr.size, true);

There is also a segfault that probably originates from accessing wfc_c before initialization.

[rajdakin]

The problem was not due to the fact that I accessed wfc_c before initialization (as it was initialized), but rather to the fact that I assigned values to an address that was not allocated (wfc_c.reserve(...); fixes this), and there were also more of these problems.
I also needed to use push_back as it didn't want to assign values the way I tried to.

[leios]

It looks like this problem has been fixed? I just compiled and it worked just fine.

[leios]

Would it be possible to put the energy calculation in the quantum systems chapter as well? It's a small modification to what you currently have.

[nic-hartley]

The problem was not due to the fact that I accessed wfc_c before initialization (as it was initialized), but rather to the fact that I assigned values to an address that was not allocated (wfc_c.reserve(...); fixes this), and there were also more of these problems.
I also needed to use push_back as it didn't want to assign values the way I tried to.

Uh, no. .reserve doesn't fix it; it's still UB to access beyond the bounds of the vector. You want to resize, or better yet, just call the ctor with the desired size. All reserve does is keep you from having to reallocate multiple times while pushing back.

[rajdakin]

@berquist I fixed the problems, it works now.
@leios OK, I'm going to put the code there too.
@nic-hartley What do you mean by UB? I tried assign the vector's value using vector[i] = ...; and it "worked" (it only crashed later when doing fftw_execute as the vectors were empty).

[nic-hartley]

@rajdakin You can also assign and read outside the bounds of a normal array without crashing; that's how Heartbleed worked. It's still UB. All reserve does is increase the capacity, which is the number of elements the vector can hold without needing to reallocate; you should still push_back into a reserved vector. Or you can call the constructor, set its initial size, and assign in. Either one will work.

[leios]

Hey, what is the status on this PR? Once this is merged, we can attack #476 . I would like to point out that the code seems to output that the energy is 0 right now, which shouldn't be correct. I can look into this more, so let me know if you want me to do this.

As a note: I am not sure what the best solution to the Undefined Behaviour (UB) is, but I can also take a swing at that, if desired.

[berquist]

There's also an incredible amount of whitespace in the file that needs to be removed.

[rajdakin]

@berquist I also removed all the unused whitespaces.

[berquist]

One last thing.

[leios]

It looks like this review is mostly done. After this, I'll double check that the energy calculation is consistent in #476 and attack that next.

[berquist]

All the code looks good. One last book thing.

[leios]

Hey @berquist sorry to bug you, but the changes also seem to be made. Is this one good to go?