三维真正交分解: C++实现

Question

我对C++非常陌生，并使用本征库来移植一个古老而缓慢的python代码，该代码从一组时态点云数据中计算出三维适当的正交分解。

在这种情况下，数据是从CFD (计算流体力学)模拟获得的云106x60x60点，每一点，三个速度分量被写入。一个例子可以看到这里。它存储在input/OFout文件夹中。从times.txt文件中检索与每个云相关联的物理时间。

该数据连接在一个大矩阵中，用于计算投影矩阵、获取特征值和特征向量以及计算POD模式等数学运算。最后，将这些模式写入输出/模式文件夹中的不同文件。

我的目标是改善两个C++方面(语法，使它更容易阅读，更惯用，面向对象.)以及代码的性能。关于如何并行化的提示/提示也会很好。

我编写的代码如下：

#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;
}

谢谢你的洞察力！

Answer 1

1. 用#define代替enum或constexpr来表示常量。
2. 我不认为时间数据或点的数量应该是绝对常量，而是应该在运行时确定它们。否则，每次收到新数据时，都必须重新编译代码。
3. 最好声明某种类型的输出路径，而不是依赖../，还需要确保目录的存在，否则它将无法保存。
4. 您不希望在大量的文件和目录上存储小数据。最好为自己制定一个序列化标准，并在每个文件中存储相当数量的数据。此外，考虑二进制保存/加载以获得更好的速度、准确性和一致性，但是，对于人类来说，它的可读性会更低，而且可移植性可能会受到一些影响--一些奇怪的平台使用较少的公共安全性。
5. 您可以加快循环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()；}}注意:通过dot函数，我希望看到标量公式，所以我看到你在里面传递向量是奇怪的。使用运算符*进行矩阵乘法。
6. 我对这些矩阵类不太熟悉，但最有可能的是，如果您替换podx.col(i) = podx.col(i) + (1.0 / (eigval(i) * TSIZE)) * sqrt(eigval(i) * TSIZE) * eigvec(j，i) *eigvec(j，i) * m.block(0 * MSIZE，j)；用podx.col(i) += (eigvec(j，i)/ sqrt(eigval(i) *TSIZE)* m.block(0 * MSIZE，j)；另外，如果eigval(i) == 0.太近，该怎么办？
7. 如果您计划扩展您的项目，那么最好最终使用一个记录器类来替换std::cout，因为std::cout并不是线程友好的，尽管在技术上是线程安全的。同样，在这种情况下，将代码移到类中，而每个相关部分都会在一个函数中移动一个合理的名称。