4
\$\begingroup\$

I'm very new to C++ and made this code using the Eigen library to port an old and slow python code which computes the three-dimensional proper orthogonal decomposition (POD) from a set of temporal point cloud data.

In this case the data is a obtained from a CFD (computational fluid dynamics) simulation as a cloud of 106x60x60 points and for each point three velocity components are written. An example can be seen here. It is stored in the input/OFout folder. The physical time associated with each cloud is retrieved from the times.txt file.

This data is concatenated in a big matrix which is used for several mathematical operations like calculating a projection matrix, obtaining the eigenvalues and eigenvectors of it and calculating the POD modes. Finally, these modes are written out to different files in the output/mode folder.

My aim is to improve both the C++ aspect (syntax, making it easier to read, more idiomatic, object oriented...) as well as the performance of the code. Tips/pointers regarding how to parallelise it would also be nice.

The code I wrote is the following:

#include <iostream> #include <string> #include <fstream> #include <iomanip> #include <Eigen/Dense> #define MSIZE 381600 // Rows of matrix (number of points) #define TSIZE 1804 // Size of time data (number of snapshots) #define VSIZE 3 // Size of the variable (1 if scalar, 3 if vector) #define NSIZE 20 // Size of output POD modes (number of modes to write) using namespace Eigen; int main() { MatrixXd m = MatrixXd::Zero(MSIZE * VSIZE, TSIZE); MatrixXd pm = MatrixXd::Zero(TSIZE, TSIZE); // READING INPUT FILES std::string t; std::string timedir = "../input/times.txt"; std::ifstream timefile(timedir); std::cout << t << std::endl; for (size_t k = 0; k < TSIZE; k++) { if (timefile.is_open()) { timefile >> t; while (t.back() == '0') // Remove trailing zeros { t.pop_back(); } } std::string dir = "../input/OFout/" + t + "/pointCloud_U.xy"; std::ifstream file(dir); if (file.is_open()) { for (size_t i = 0; i < MSIZE; i++) { for (size_t j = 0; j < VSIZE; j++) { file >> m(i + MSIZE * j, k); } } std::cout << "Finished reading file " + std::to_string(k + 1) + " of " + std::to_string(TSIZE) << std::endl; file.close(); } else { std::cout << "Unable to open file" << std::endl; } } // COMPUTING NORMALISED PROJECTION MATRIX std::cout << "Computing projection matrix..." << std::flush; for (size_t i = 0; i < TSIZE; i++) { for (size_t j = 0; j < TSIZE; j++) { pm(i, j) = (1.0 / TSIZE) * (m.col(i).dot(m.col(j).transpose())); } } std::cout << " Done" << std::endl; // COMPUTING SORTED EIGENVALUES AND EIGENVECTORS std::cout << "Computing eigenvalues..." << std::flush; SelfAdjointEigenSolver<MatrixXd> eigensolver(pm); if (eigensolver.info() != Success) abort(); VectorXd eigval = eigensolver.eigenvalues().reverse(); MatrixXd eigvec = eigensolver.eigenvectors().rowwise().reverse(); std::cout << " Done" << std::endl; // COMPUTING POD MODES std::cout << "Computing POD modes..." << std::flush; MatrixXd podx = MatrixXd::Zero(MSIZE, TSIZE); MatrixXd pody = MatrixXd::Zero(MSIZE, TSIZE); MatrixXd podz = MatrixXd::Zero(MSIZE, TSIZE); for (size_t i = 0; i < TSIZE; i++) { for (size_t j = 0; j < TSIZE; j++) { podx.col(i) = podx.col(i) + (1.0 / (eigval(i) * TSIZE)) * sqrt(eigval(i) * TSIZE) * eigvec(j, i) * m.block<MSIZE, 1>(0 * MSIZE, j); pody.col(i) = pody.col(i) + (1.0 / (eigval(i) * TSIZE)) * sqrt(eigval(i) * TSIZE) * eigvec(j, i) * m.block<MSIZE, 1>(1 * MSIZE, j); podz.col(i) = podz.col(i) + (1.0 / (eigval(i) * TSIZE)) * sqrt(eigval(i) * TSIZE) * eigvec(j, i) * m.block<MSIZE, 1>(2 * MSIZE, j); } } std::cout << " Done" << std::endl; // WRITING SORTED EIGENVALUES std::cout << "Writing eigenvalues..." << std::flush; std::ofstream writeEigval("../output/chronos/A.txt"); if (writeEigval.is_open()) { writeEigval << std::scientific << std::setprecision(10) << eigval; writeEigval.close(); } std::cout << " Done" << std::endl; // WRITING CHRONOS std::cout << "Writing chronos..." << std::flush; for (size_t i = 0; i < NSIZE; i++) { std::string chronos = "../output/chronos/chronos_" + std::to_string(i) + ".dat"; std::ofstream writeChronos(chronos); for (size_t j = 0; j < TSIZE; j++) { writeChronos << std::scientific << std::setprecision(10) << sqrt(eigval(i) * TSIZE) * eigvec(j, i) << '\n'; } writeChronos.close(); } std::cout << " Done" << std::endl; // WRITING POD MODES std::cout << "Writing POD modes..." << std::flush; std::string xyz = "xyz_"; std::string modeHead = "../output/mode/mode_U"; for (size_t j = 0; j < VSIZE; j++) { for (size_t i = 0; i < NSIZE; i++) { std::string modeTail = std::to_string(i) + ".dat"; std::ofstream writeMode(modeHead + xyz.at(j) + xyz.at(3) + modeTail); if (writeMode.is_open()) { if (j == 0) { writeMode << std::scientific << std::setprecision(6) << podx.col(i); } else if (j == 1) { writeMode << std::scientific << std::setprecision(6) << pody.col(i); } else if (j == 2) { writeMode << std::scientific << std::setprecision(6) << podz.col(i); } writeMode.close(); } } } std::cout << " Done" << std::endl; }

Thanks for your insights!

\$\endgroup\$

1 Answer 1

Reset to default
2
\$\begingroup\$
  1. Instead of #define use enum or constexpr for the constant numbers.
  2. I don't think that the number of time data or number of points should be absolute constants - rather they should be determined during runtime. Else you have to recompile your code each time you receive new data.
  3. Better declare some sort of output path instead of relying on ../, also you need to make sure that the directories exist, else it will fail to save.
  4. You don't want to store small pieces of data over huge amount of files and directories. Better make a serialization standard for yourselves and store a sizable amount of data inside each file. Also, consider binary save/load for better speed, accuracy, and consistency, though, it will be less readable for humans and portability might suffer a bit - some odd platforms use less common endianess.

  5. You can speed up the loop

    for (size_t i = 0; i < TSIZE; i++) { for (size_t j = 0; j <= i; j++) { pm(j, i) = pm(i, j) = (1.0 / TSIZE) * (m.col(i).dot(m.col(j).transpose())); } }

    Note: by dot function I expect to see scalar pruduct, so it is odd for me to see that you transpose the vector inside it. Use operator * for matrix multiplication.

  6. I am not too familiar with these Matrices classes, but most likely, it will be better if you replace

    podx.col(i) = podx.col(i) + (1.0 / (eigval(i) * TSIZE)) * sqrt(eigval(i) * TSIZE) * eigvec(j, i) * m.block<MSIZE, 1>(0 * MSIZE, j);

    with

    podx.col(i) += (eigvec(j, i) / sqrt(eigval(i) * TSIZE)) * m.block<MSIZE, 1>(0 * MSIZE, j);

    also, what to do if eigval(i) == 0. or too close to it?

  7. If you plan to grow your project you'd better eventually replace std::cout with a usage of a logger class as std::cout is not thread friendly despite being technically thread safe. Also, in this, case move the code inside a class and each relevant section move inside a function with a sensible name.
\$\endgroup\$
4
  • \$\begingroup\$ "Instead of #define use enum for the constant numbers." No, use constexpr. \$\endgroup\$ Commented Oct 12, 2019 at 10:47
  • \$\begingroup\$ @L.F. I see constexpr being better when you require computation or non trivial classes. But for simple integer constants? \$\endgroup\$ Commented Oct 12, 2019 at 14:22
  • \$\begingroup\$ Still better. For both semantics and consistency. \$\endgroup\$ Commented Oct 12, 2019 at 14:31
  • \$\begingroup\$ @L.F. edited and made recommendation for both. \$\endgroup\$ Commented Oct 12, 2019 at 14:35

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.