Revise configuration, use structured settings

TLCFEM · Sep 21, 2024 · 863fa65 · 863fa65
1 parent b0efd1b
commit 863fa65
Show file tree

Hide file tree

Showing 4 changed files with 153 additions and 140 deletions.
diff --git a/README.md b/README.md
@@ -159,16 +159,16 @@ Usage: vpmr [options]
 
 Options:
 
-   -n <int>     number of terms (default: 10)
-   -c <int>     maximum exponent (default: 4)
-   -d <int>     number of precision bits (default: 512)
-   -q <int>     quadrature order (default: 500)
-   -m <float>   precision multiplier (default, minimum: 1.5)
-   -e <float>   tolerance (default: 1E-8)
-   -k <string>  file name of kernel function (default: exp(-t^2/4))
-   -s           print singular values
-   -w           print weights
-   -h           print this help message
+    -n, --max-terms             <int>     number of terms (default: 10)
+    -c, --max-exponent          <int>     maximum exponent (default: 4)
+    -d, --precision-bits        <int>     number of precision bits (default: 512)
+    -q, --quadrature-order      <int>     quadrature order (default: 500)
+    -m, --precision-multiplier  <float>   precision multiplier (default: 1.5)
+    -e, --tolerance             <float>   tolerance (default: 1E-8)
+    -k, --kernel                <string>  file name of kernel function (default uses: exp(-t^2/4))
+    -s, --singular-values                 print singular values
+    -w, --weights                         print weights
+    -h, --help                            print this help message
 ```
 
 The minimum required precision can be estimated by the parameter $n$.
@@ -188,12 +188,13 @@ The output is:
 
 ```text
 Using the following parameters:
-     terms = 30.
-     order = 500.
-  exponent = 4.
- precision = 336.
- tolerance = 1.0000e-08.
-    kernel = exp(-t*t/4).
+       terms = 30.
+    exponent = 4.
+   precision = 355.
+ quad. order = 500.
+  multiplier = 1.5000e+00.
+   tolerance = 1.0000e-08.
+      kernel = exp(-t*t/4).
 
 [1/6] Computing weights... [60/60]
 [2/6] Solving Lyapunov equation...
@@ -203,7 +204,7 @@ Using the following parameters:
 [6/6] Done.
 
 M = 
-+1.1745193571738943e+01-1.4261645574068720e-100j
++1.1745193571738943e+01+6.4089561283054790e-107j
 -5.5143304351134397e+00+5.7204056791636839e+00j
 -5.5143304351134397e+00-5.7204056791636839e+00j
 -1.6161617424833762e-02+2.3459542440459513e+00j
@@ -223,7 +224,7 @@ S =
 +1.7655956664692953e+00-2.7555720406099038e+00j
 +1.7655956664692953e+00+2.7555720406099038e+00j
 
-Running time: 3 s.
+Running time: 3112 ms.
 ```
 
 ![exp(-t^2/4)](resource/example.png)

diff --git a/src/Gauss.hpp b/src/Gauss.hpp
@@ -18,22 +18,19 @@
 #pragma once
 
 #include <cmath>
-#include <iomanip>
 #include <memory>
 #include <tbb/blocked_range.h>
 #include <tbb/parallel_for.h>
 #include <tbb/parallel_reduce.h>
-#include "mpreal.h"
-
-using namespace mpfr;
 
 inline mpreal MP_PI = const_pi();
 inline mpreal MP_PI_HALF{MP_PI / 2};
-extern int DIGIT;
+
+extern Config config;
 
 class Evaluation {
-    const mpreal ONE = mpreal(1, DIGIT);
-    const mpreal TWO = mpreal(2, DIGIT);
+    const mpreal ONE = mpreal(1, config.precision_bits);
+    const mpreal TWO = mpreal(2, config.precision_bits);
 
     const size_t degree;
 
@@ -43,13 +40,13 @@ class Evaluation {
     explicit Evaluation(const mpreal& X, const size_t D)
         : degree(D)
         , _x(X)
-        , _v(1, DIGIT)
-        , _d(0, DIGIT) { this->evaluate(X); }
+        , _v(1, config.precision_bits)
+        , _d(0, config.precision_bits) { this->evaluate(X); }
 
     void evaluate(const mpreal& x) {
         this->_x = x;
 
-        mpreal left{x}, right{1, DIGIT};
+        mpreal left{x}, right{1, config.precision_bits};
 
         for(size_t i = 2; i <= degree; ++i) {
             this->_v = ((TWO * i - ONE) * x * left - (i - ONE) * right) / i;
@@ -69,8 +66,8 @@ class Evaluation {
 };
 
 class LegendrePolynomial {
-    const mpreal ONE = mpreal(1, DIGIT);
-    const mpreal TWO = mpreal(2, DIGIT);
+    const mpreal ONE = mpreal(1, config.precision_bits);
+    const mpreal TWO = mpreal(2, config.precision_bits);
 
 public:
     const size_t degree;
@@ -84,9 +81,9 @@ class LegendrePolynomial {
         , _r(std::make_unique<mpreal[]>(degree))
         , _w(std::make_unique<mpreal[]>(degree)) {
         tbb::parallel_for(static_cast<size_t>(0), degree / 2 + 1, [&](const size_t i) {
-            mpreal dr{1, DIGIT};
+            mpreal dr{1, config.precision_bits};
 
-            Evaluation eval(cos(MP_PI * mpreal(4 * i + 3, DIGIT) / mpreal(4 * degree + 2, DIGIT)), degree);
+            Evaluation eval(cos(MP_PI * mpreal(4 * i + 3, config.precision_bits) / mpreal(4 * degree + 2, config.precision_bits)), degree);
             do {
                 dr = eval.v() / eval.d();
                 eval.evaluate(eval.x() - dr);
@@ -118,12 +115,12 @@ class Quadrature {
         : poly(D) {}
 
     template<typename Function> mpreal integrate(Function&& f) const {
-        // mpreal sum{0, DIGIT};
+        // mpreal sum{0, config.precision_bits};
         // 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<int>(0, static_cast<int>(poly.degree)), mpreal(0, DIGIT),
+            tbb::blocked_range<int>(0, static_cast<int>(poly.degree)), mpreal(0, config.precision_bits),
             [&](const tbb::blocked_range<int>& r, mpreal running_total) {
                 for(auto i = r.begin(); i < r.end(); ++i) running_total += poly.weight(i) * f(i, poly.root(i));
                 return running_total;