From e640588cffeedb724310d2ccabd2921cb37883d5 Mon Sep 17 00:00:00 2001 From: Theodore Chang Date: Tue, 28 Nov 2023 23:35:39 +0100 Subject: [PATCH] Reformat --- src/Gauss.hpp | 30 +++++++++++++----------- src/VPMR.cpp | 65 ++++++++++++++++----------------------------------- 2 files changed, 36 insertions(+), 59 deletions(-) diff --git a/src/Gauss.hpp b/src/Gauss.hpp index 3b42e11..03b142b 100644 --- a/src/Gauss.hpp +++ b/src/Gauss.hpp @@ -19,7 +19,6 @@ #include #include -#include #include #include #include @@ -28,11 +27,11 @@ using namespace mpfr; -mpreal MP_PI = [] { +inline mpreal MP_PI = [] { mpreal::set_default_prec(512); return const_pi(); }(); -mpreal MP_PI_HALF{MP_PI / 2}; +inline mpreal MP_PI_HALF{MP_PI / 2}; extern int DIGIT; class Evaluation { @@ -42,7 +41,10 @@ class Evaluation { public: explicit Evaluation(const mpreal& X, const size_t D) - : degree(D), _x(X), _v(1, DIGIT), _d(0, DIGIT) { this->evaluate(X); } + : degree(D) + , _x(X) + , _v(1, DIGIT) + , _d(0, DIGIT) { this->evaluate(X); } void evaluate(const mpreal& x) { this->_x = x; @@ -75,8 +77,10 @@ class LegendrePolynomial { public: explicit LegendrePolynomial(const size_t D) - : degree(D > 2 ? D : 2), _r(std::make_unique(degree)), _w(std::make_unique(degree)) { - tbb::parallel_for(size_t(0), degree / 2 + 1, [&](const size_t i) { + : degree(D > 2 ? D : 2) + , _r(std::make_unique(degree)) + , _w(std::make_unique(degree)) { + tbb::parallel_for(static_cast(0), degree / 2 + 1, [&](const size_t i) { mpreal dr{1, DIGIT}; Evaluation eval(cos(MP_PI * mpreal(4 * i + 3, DIGIT) / mpreal(4 * degree + 2, DIGIT)), degree); @@ -95,12 +99,12 @@ class LegendrePolynomial { this->_w[i] = this->_w[degree - i - 1]; }); - tbb::parallel_for(size_t(0), degree, [&](const size_t i) { this->_r[i] = MP_PI_HALF * this->_r[i] + MP_PI_HALF; }); + tbb::parallel_for(static_cast(0), degree, [&](const size_t i) { this->_r[i] = MP_PI_HALF * this->_r[i] + MP_PI_HALF; }); } - [[nodiscard]] mpreal root(int i) const { return this->_r[i]; } + [[nodiscard]] mpreal root(const int i) const { return this->_r[i]; } - [[nodiscard]] mpreal weight(int i) const { return this->_w[i]; } + [[nodiscard]] mpreal weight(const int i) const { return this->_w[i]; } }; class Quadrature { @@ -110,17 +114,15 @@ class Quadrature { explicit Quadrature(const size_t D) : poly(D) {} - template - mpreal integrate(Function&& f) const { + template mpreal integrate(Function&& f) const { // mpreal sum{0, DIGIT}; // for (int i = 0; i < int(poly.degree); ++i) // sum += poly.weight(i) * f(i, poly.root(i)); // return sum; return tbb::parallel_deterministic_reduce( - tbb::blocked_range(0, int(poly.degree)), mpreal(0, DIGIT), + tbb::blocked_range(0, static_cast(poly.degree)), mpreal(0, DIGIT), [&](const tbb::blocked_range& r, mpreal running_total) { - for(auto i = r.begin(); i < r.end(); ++i) - running_total += poly.weight(i) * f(i, poly.root(i)); + for(auto i = r.begin(); i < r.end(); ++i) running_total += poly.weight(i) * f(i, poly.root(i)); return running_total; }, std::plus<>()); diff --git a/src/VPMR.cpp b/src/VPMR.cpp index 396b099..e196310 100644 --- a/src/VPMR.cpp +++ b/src/VPMR.cpp @@ -82,9 +82,10 @@ class Expression { } }; +Expression kernel{}; + // integrand in eq. 2.2 --- K(t)cos(jt) class Integrand { - static Expression* expression; static tbb::concurrent_unordered_map value; const int j; @@ -93,16 +94,12 @@ class Integrand { explicit Integrand(const int J) : j(J) {} - static void set_expression(Expression* E) { expression = E; } - mpreal operator()(const int i, const mpreal& r) const { - if(value.find(i) == value.end()) - value.insert(std::make_pair(i, expression->value(-NC * log((mpreal(1, DIGIT) + cos(r)) >> 1)))); + if(value.find(i) == value.end()) value.insert(std::make_pair(i, kernel.value(-NC * log((mpreal(1, DIGIT) + cos(r)) >> 1)))); return value[i] * cos(j * r); } }; -Expression* Integrand::expression = nullptr; tbb::concurrent_unordered_map Integrand::value(QUAD_ORDER); using mat = Eigen::Matrix; @@ -111,9 +108,7 @@ using cx_vec = Eigen::Matrix, Eigen::Dynamic, 1>; // step 2 mat lyap(const vec& A, mat&& C) { - tbb::parallel_for(0, static_cast(A.size()), [&](const int i) { - for(auto j = 0; j < A.size(); ++j) C(i, j) /= (A(i) + A(j)); - }); + tbb::parallel_for(0, static_cast(A.size()), [&](const int i) { for(auto j = 0; j < A.size(); ++j) C(i, j) /= (A(i) + A(j)); }); return std::move(C); } @@ -164,14 +159,12 @@ std::tuple model_reduction(const vec& A, const vec& B, const vec // step 7 const auto P = pos(SIG); std::cout << "[4/6] Transforming (P=" << P << ")...\n"; - if(P == SIG.size()) - std::cout << "WARNING: No singular value is smaller than the given tolerance.\n"; + if(P == SIG.size()) std::cout << "WARNING: No singular value is smaller than the given tolerance.\n"; // step 5 auto SS = SIG; for(auto I = 0; I < SS.size(); ++I) - if(SS(I) > std::numeric_limits::epsilon()) - SS(I) = pow(SS(I), -.5); + if(SS(I) > std::numeric_limits::epsilon()) SS(I) = pow(SS(I), -.5); else { std::cout << "WARNING: Need to increase digits.\n"; SS(I) = pow(std::numeric_limits::epsilon(), -.5); @@ -266,17 +259,15 @@ std::tuple vpmr() { auto [M, S] = model_reduction(A, B, C); auto ID = sort_index(M); - for(int I = int(ID.size()) - 1; I >= 0; --I) - if(norm(M(ID[I])) < TOL) - ID.erase(ID.begin() + I); - else - break; + for(int I = static_cast(ID.size()) - 1; I >= 0; --I) + if(norm(M(ID[I])) < TOL) ID.erase(ID.begin() + I); + else break; std::cout << "[6/6] Done.\n\n"; if(abs(W(0)) < TOL) return std::make_tuple(M(ID), -S(ID)); - cx_vec MM = cx_vec::Zero(long(ID.size()) + 1), SS = cx_vec::Zero(long(ID.size()) + 1); + cx_vec MM = cx_vec::Zero(static_cast(ID.size()) + 1), SS = cx_vec::Zero(static_cast(ID.size()) + 1); MM(0) = W(0); SS(0) = mpreal(0, DIGIT); MM.tail(ID.size()) = M(ID); @@ -313,28 +304,18 @@ int main(const int argc, const char** argv) { bool has_digit = false; for(auto I = 1; I < argc; ++I) { - const auto token = std::string(argv[I]); - - if(token == "-nc") - NC = std::stoi(argv[++I]); - else if(token == "-n") - N = std::stoi(argv[++I]); + if(const auto token = std::string(argv[I]); token == "-nc") NC = std::stoi(argv[++I]); + else if(token == "-n") N = std::stoi(argv[++I]); else if(token == "-d") { DIGIT = std::stoi(argv[++I]); has_digit = true; } - else if(token == "-q") - QUAD_ORDER = std::stoi(argv[++I]); - else if(token == "-m") - SCALE = std::stoi(argv[++I]); - else if(token == "-e") - TOL = mpreal(argv[++I]); - else if(token == "-w") - OUTPUT_W = true; - else if(token == "-s") - OUTPUT_EIG = true; - else if(token == "-h") - return print_helper(); + else if(token == "-q") QUAD_ORDER = std::stoi(argv[++I]); + else if(token == "-m") SCALE = std::stoi(argv[++I]); + else if(token == "-e") TOL = mpreal(argv[++I]); + else if(token == "-w") OUTPUT_W = true; + else if(token == "-s") OUTPUT_EIG = true; + else if(token == "-h") return print_helper(); else if(token == "-k") { const std::string file_name(argv[++I]); std::ifstream file(file_name); @@ -359,8 +340,7 @@ int main(const int argc, const char** argv) { while((comb_max /= 2) > 0) ++comb_digit; comb_digit = std::max(5 * N, SCALE * comb_digit); - if(!has_digit) - DIGIT = comb_digit; + if(!has_digit) DIGIT = comb_digit; else if(comb_digit >= DIGIT) { std::cout << "WARNING: Too few digits to hold combinatorial number, resetting digits to " << comb_digit << ".\n"; DIGIT = comb_digit; @@ -376,12 +356,10 @@ int main(const int argc, const char** argv) { TOL /= 2; - Expression kernel; if(!kernel.compile()) { std::cerr << "Cannot compile kernel function: " << KERNEL << ".\n"; return 1; } - Integrand::set_expression(&kernel); std::cout << std::scientific << std::setprecision(4); @@ -435,8 +413,7 @@ std::tuple>, std::vector>> while((comb_max /= 2) > 0) ++comb_digit; comb_digit = std::max(5 * N, SCALE * comb_digit); - if(0 == d) - DIGIT = comb_digit; + if(0 == d) DIGIT = comb_digit; else if(comb_digit >= DIGIT) { std::cout << "WARNING: Too few digits to hold combinatorial number, resetting digits to " << comb_digit << ".\n"; DIGIT = comb_digit; @@ -454,12 +431,10 @@ std::tuple>, std::vector>> std::vector> mm, ss; - Expression kernel; if(!kernel.compile()) { std::cerr << "Cannot compile kernel function: " << KERNEL << ".\n"; return {mm, ss}; } - Integrand::set_expression(&kernel); const auto [M, S] = vpmr();