markcmiller86 · frobnitzem · Jun 10, 2021 · Jun 10, 2021 · Jun 10, 2021 · Jun 10, 2021
diff --git a/Number.H b/Number.H
@@ -1,4 +1,5 @@
 #include <ostream>
+#include <vector>
 
 #include "Half.H"
 
@@ -8,109 +9,43 @@
 
 // 0=half, 1=float, 2=double, 3=quad
 #if   FPTYPE == 0
-typedef half_float::half fpnumber;
+typedef half_float::half Number;
+template <typename T> Number castNum(T x) {
+    return half_float::half_cast<half_float::half>(x);
+}
 #elif FPTYPE == 1
-typedef float fpnumber;
+typedef float Number;
+template <typename T> Number castNum(T x) { return (Number) x; }
 #elif FPTYPE == 2
-typedef double fpnumber;
+typedef double Number;
+template <typename T> Number castNum(T x) { return (Number) x; }
 #elif FPTYPE == 3
-typedef long double fpnumber;
+typedef long double Number;
+template <typename T> Number castNum(T x) { return (Number) x; }
 #else
 #error UNRECOGNIZED FPTYPE
 #endif
 
-// Encapsulate a floating point number to easily track various metrics.
-// Note: Doing this well is complicated by the mix of types in numerical
-// statements (e.g. int*double) made explicit by modern compilers.
-class Number {
-  public:
-    static void *operator new(std::size_t sz) { Number::nbytes += sz; return ::operator new(sz); };
-    static void *operator new[](std::size_t sz) { Number::nbytes += sz; return ::operator new(sz); };
-    static char const *counts_string() { static char counts[256];
-        sprintf(counts, "Adds:%d, Mults:%d, Divs:%d, Bytes:%d", Number::nadds, Number::nmults, Number::ndivs, (int) Number::nbytes);
-        return counts; };
-    inline Number() : x(0) {};
-    inline Number(fpnumber _x) : x(_x) {};
-#if FPTYPE!=2
-    inline Number(double _x) : x((fpnumber) _x) {};
-#endif
-    inline Number(int _x) : x((fpnumber) _x) {};
+class Vector {
+    std::vector<Number> x;
+    const int n;
 
-    inline Number &operator=(const Number& rhs) { this->x = rhs.x; return *this; };
-    inline operator fpnumber() const { return x; };
-    inline operator int() const { return static_cast<int>(x); };
-#if FPTYPE!=2
-    inline operator double() const { return static_cast<double>(x); };
-#endif
-  private:
-    static int nadds;
-    static int nmults;
-    static int ndivs;
-    static std::size_t nbytes;
-    fpnumber x;
+    public:
+        Vector(int n_) : x(n_), n(n_) { }
 
-    // Various operators on Number w/mix of int
-    friend Number operator+(const Number& lhs, const Number& rhs) {  Number::nadds++; return lhs.x + rhs.x; };
-    friend Number operator+(const int& lhs, const Number& rhs) {  Number::nadds++; return lhs + rhs.x; };
-    friend Number operator+(const Number& lhs, const int& rhs) {  Number::nadds++; return lhs.x + rhs; };
-    friend Number operator+(const double& lhs, const Number& rhs) {  Number::nadds++; return lhs + rhs.x; };
-    friend Number operator+(const Number& lhs, const double& rhs) {  Number::nadds++; return lhs.x + rhs; };
-    friend Number operator+=(Number& lhs, const Number& rhs) {  Number::nadds++; return lhs.x += rhs.x; };
-    friend Number operator+=(Number& lhs, const double& rhs) {  Number::nadds++; return lhs.x += rhs; };
-    friend Number operator+=(double& lhs, const Number& rhs) {  Number::nadds++; return lhs += rhs.x; };
-    friend Number operator-(const Number& lhs, const Number& rhs) {  Number::nadds++; return lhs.x - rhs.x; };
-    friend Number operator-(const int& lhs, const Number& rhs) {  Number::nadds++; return lhs - rhs.x; };
-    friend Number operator-(const Number& lhs, const int& rhs) {  Number::nadds++; return lhs.x - rhs; };
-    friend Number operator-(const double& lhs, const Number& rhs) {  Number::nadds++; return lhs - rhs.x; };
-    friend Number operator-(const Number& lhs, const double& rhs) {  Number::nadds++; return lhs.x - rhs; };
-    friend Number operator-(const Number& rhs) {  Number::nadds++; return -rhs.x; };
-    friend Number operator-=(Number& lhs, const Number& rhs) {  Number::nadds++; return lhs.x -= rhs.x; };
-    friend Number operator*(const Number& lhs, const Number& rhs) { Number::nmults++; return lhs.x * rhs.x; };
-    friend Number operator*(const int& lhs, const Number& rhs) { Number::nmults++; return lhs * rhs.x; };
-    friend Number operator*(const Number& lhs, const int& rhs) { Number::nmults++; return lhs.x * rhs; };
-    friend Number operator*(const double& lhs, const Number& rhs) { Number::nmults++; return lhs * rhs.x; };
-    friend Number operator*(const Number& lhs, const double& rhs) { Number::nmults++; return lhs.x * rhs; };
-    friend Number operator*=(Number& lhs, const Number& rhs) { Number::nmults++; return lhs.x *= rhs.x; };
-    friend Number operator/(const Number& lhs, const Number& rhs) { Number::ndivs++; return lhs.x / rhs.x; };
-    friend Number operator/(const int& lhs, const Number& rhs) { Number::ndivs++; return lhs / rhs.x; };
-    friend Number operator/(const Number& lhs, const int& rhs) { Number::ndivs++; return lhs.x / rhs; };
-    friend Number operator/(const double& lhs, const Number& rhs) { Number::ndivs++; return lhs / rhs.x; };
-    friend Number operator/(const Number& lhs, const double& rhs) { Number::ndivs++; return lhs.x / rhs; };
-    friend Number operator/=(Number& lhs, const Number& rhs) { Number::ndivs++; return lhs.x /= rhs.x; };
-    friend bool operator< (const Number& lhs, const Number& rhs){ return lhs.x < rhs.x; }
-    friend bool operator< (const int& lhs, const Number& rhs){ return lhs < rhs.x; }
-    friend bool operator< (const Number& lhs, const int& rhs){ return lhs.x < rhs; }
-    friend bool operator< (const double& lhs, const Number& rhs){ return lhs < rhs.x; }
-    friend bool operator< (const Number& lhs, const double& rhs){ return lhs.x < rhs; }
-    friend bool operator> (const Number& lhs, const Number& rhs){ return rhs < lhs; }
-    friend bool operator> (const int& lhs, const Number& rhs){ return rhs < lhs; }
-    friend bool operator> (const Number& lhs, const int& rhs){ return rhs < lhs; }
-    friend bool operator> (const double& lhs, const Number& rhs){ return rhs < lhs; }
-    friend bool operator> (const Number& lhs, const double& rhs){ return rhs < lhs; }
-    friend bool operator<=(const Number& lhs, const Number& rhs){ return !(lhs > rhs); }
-    friend bool operator<=(const int& lhs, const Number& rhs){ return !(lhs > rhs); }
-    friend bool operator<=(const Number& lhs, const int& rhs){ return !(lhs > rhs); }
-    friend bool operator<=(const double& lhs, const Number& rhs){ return !(lhs > rhs); }
-    friend bool operator<=(const Number& lhs, const double& rhs){ return !(lhs > rhs); }
-    friend bool operator>=(const Number& lhs, const Number& rhs){ return !(lhs < rhs); }
-    friend bool operator>=(const int& lhs, const Number& rhs){ return !(lhs < rhs); }
-    friend bool operator>=(const Number& lhs, const int& rhs){ return !(lhs < rhs); }
-    friend bool operator>=(const double& lhs, const Number& rhs){ return !(lhs < rhs); }
-    friend bool operator>=(const Number& lhs, const double& rhs){ return !(lhs < rhs); }
-    friend bool operator==(const Number& lhs, const Number& rhs){ return lhs.x == rhs.x; }
-    friend bool operator==(const int& lhs, const Number& rhs){ return lhs == rhs.x; }
-    friend bool operator==(const Number& lhs, const int& rhs){ return lhs.x == rhs; }
-    friend bool operator==(const double& lhs, const Number& rhs){ return lhs == rhs.x; }
-    friend bool operator==(const Number& lhs, const double& rhs){ return lhs.x == rhs; }
-    friend bool operator!=(const Number& lhs, const Number& rhs){ return !(lhs == rhs); }
-    friend bool operator!=(const int& lhs, const Number& rhs){ return !(lhs == rhs); }
-    friend bool operator!=(const Number& lhs, const int& rhs){ return !(lhs == rhs); }
-    friend bool operator!=(const double& lhs, const Number& rhs){ return !(lhs == rhs); }
-    friend bool operator!=(const Number& lhs, const double& rhs){ return !(lhs == rhs); }
-    friend std::ostream& operator<<(std::ostream& os, const Number& rhs)  { os << rhs.x; return os; }
+        Number& operator[](int i) {
+            return x[i];
+        }
+        const Number operator[](int i) const {
+            return x[i];
+        }
+        int size() const {
+            return n;
+        }
+        const Number *data() const {
+            return x.data();
+        }
+        Number *data() {
+            return x.data();
+        }
 };
-
-#define TSTART -1
-#define TFINAL -2
-#define RESIDUAL -3
-#define ERROR -4
diff --git a/args.C b/args.C
@@ -6,7 +6,7 @@ static char clargs[2048];
 { \
     char const *style = #STYLE; \
     char const *q = style[1]=='s'?"\"":""; \
-    void *valp = (void*) &VAR; \
+    void *valp = (void*) &ret.VAR; \
     int const len = strlen(#VAR)+1; \
     std::stringstream strmvar; \
     for (i = 1; i < argc; i++) \
@@ -20,12 +20,12 @@ static char clargs[2048];
             if      (style[1] == 'd') /* int */ \
                 *((int*) valp) = (int) strtol(argv[i]+len,0,10); \
             else if (style[1] == 'g') /* double */ \
-                *((Number*) valp) = (fpnumber) strtod(argv[i]+len,0); \
+                *((Number*) valp) = (Number) strtod(argv[i]+len,0); \
             else if (style[1] == 's') /* char* */ \
                 *((char**) valp) = (char*) strdup(argv[i]+len); \
         }\
     }\
-    strmvar << VAR; \
+    strmvar << ret.VAR; \
     if (help) \
     {\
         char tmp[256]; \
@@ -46,22 +46,7 @@ static char clargs[2048];
     } \
 }
 
-extern Number alpha;
-extern Number lenx;
-extern Number dx;
-extern Number dt;
-extern Number maxt;
-extern Number bc0;
-extern Number bc1;
-extern Number min_change;
-extern char const *runame;
-extern char const *ic;
-extern char const *alg;
-extern int savi;
-extern int save;
-extern int outi;
-extern int noout;
-int const prec = FPTYPE;
+//int const prec = FPTYPE;
 
 static void handle_help(char const *argv0)
 {
@@ -86,9 +71,8 @@ static void handle_ic_help()
         "\nconditions such that they *combine* smoothly with the initial conditions.\n");
 }
 
-void
-process_args(int argc, char **argv)
-{
+Args process_args(int argc, char **argv) {
+    Args ret{};
     int i;
     int help = 0;
 
@@ -100,26 +84,26 @@ process_args(int argc, char **argv)
     if (help)
         fprintf(stderr, "Usage: %s <arg>=<value> <arg>=<value>...\n", argv[0]);
 
-    HANDLE_ARG(alpha, fpnumber, %g, material thermal diffusivity (sq-meters/second));
-    HANDLE_ARG(lenx, fpnumber, %g, material length (meters));
-    HANDLE_ARG(dx, fpnumber, %g, x-incriment. Best if lenx/dx==int. (meters));
-    HANDLE_ARG(dt, fpnumber, %g, t-incriment (seconds));
-    HANDLE_ARG(maxt, fpnumber, %g, >0:max sim time (seconds) | <0:min l2 change in soln);
-    HANDLE_ARG(bc0, fpnumber, %g, boundary condition @ x=0: u(0,t) (Kelvin));
-    HANDLE_ARG(bc1, fpnumber, %g, boundary condition @ x=lenx: u(lenx,t) (Kelvin));
+    HANDLE_ARG(alpha, Number, %g, material thermal diffusivity (sq-meters/second));
+    HANDLE_ARG(lenx, Number, %g, material length (meters));
+    HANDLE_ARG(dx, Number, %g, x-incriment. Best if lenx/dx==int. (meters));
+    HANDLE_ARG(dt, Number, %g, t-incriment (seconds));
+    HANDLE_ARG(maxt, Number, %g, >0:max sim time (seconds) | <0:min l2 change in soln);
+    HANDLE_ARG(bc0, Number, %g, boundary condition @ x=0: u(0,t) (Kelvin));
+    HANDLE_ARG(bc1, Number, %g, boundary condition @ x=lenx: u(lenx,t) (Kelvin));
     HANDLE_ARG(runame, char*, %s, name to give run and results dir);
     HANDLE_ARG(ic, char*, %s, initial condition @ t=0: u(x,0) (Kelvin));
     HANDLE_ARG(alg, char*, %s, algorithm ftcs|dufrank|crankn);
     HANDLE_ARG(savi, int, %d, save every i-th solution step);
     HANDLE_ARG(save, int, %d, save error in every saved solution);
     HANDLE_ARG(outi, int, %d, output progress every i-th solution step);
     HANDLE_ARG(noout, int, %d, disable all file outputs);
-    HANDLE_ARG(prec, const int, %d, precision 0=half/1=float/2=double/3=long double)
+    //HANDLE_ARG(prec, const int, %d, precision 0=half/1=float/2=double/3=long double)
 
     if (help)
         handle_help(argv[0]);
 
-    if (strcasestr(ic, "help"))
+    if (strcasestr(ret.ic, "help"))
     {
         handle_ic_help();
         exit(1);
@@ -129,24 +113,31 @@ process_args(int argc, char **argv)
         exit(1);
 
     // Handle possible termination by change threshold criterion
-    if (maxt < 0)
+    if (ret.maxt < 0)
     {
-        min_change = -maxt * -maxt;
-        maxt = INT_MAX;
+        ret.min_change = -ret.maxt * -ret.maxt;
+        ret.maxt = INT_MAX;
     }
 
     // Handle output results dir creation and save of command-line
-    if (access(runame, F_OK) == 0)
+    if (access(ret.runame, F_OK) == 0)
     {
-        fprintf(stderr, "An entry \"%s\" already exists\n", runame);
+        fprintf(stderr, "An entry \"%s\" already exists\n", ret.runame);
         exit(1);
     } 
 
     // Make the output dir and save clargs there too
-    mkdir(runame, S_IRWXU|S_IRWXG|S_IROTH|S_IXOTH);
+    mkdir(ret.runame, S_IRWXU|S_IRWXG|S_IROTH|S_IXOTH);
     char fname[128];
-    sprintf(fname, "%s/clargs.out", runame);
+    sprintf(fname, "%s/clargs.out", ret.runame);
     FILE *outf = fopen(fname, "w");
     fprintf(outf, "%s", clargs);
     fclose(outf);
+
+    // final setup on args
+    ret.Nx = (int) round((double)(ret.lenx/ret.dx))+1;
+    ret.Nt = (int) (ret.maxt/ret.dt);
+    ret.dx = ret.lenx/(ret.Nx-1);
+
+    return ret;
 }
diff --git a/args.H b/args.H
@@ -0,0 +1,44 @@
+#pragma once
+
+// Command-line argument variables
+struct Args {
+    // Number of points in space, x, and time, t.
+    int Nx;
+    int Nt;
+
+    int noout;
+    int savi;
+    int outi;
+    int save;
+    char const *runame;
+    char const *alg;
+    char const *ic;
+    Number lenx;
+    Number alpha;
+    Number dt;
+    Number dx;
+    Number bc0;
+    Number bc1;
+    Number maxt;
+    Number min_change;
+
+    Args() :
+        Nx(0), Nt(0),
+        noout        ( 0 ),
+        savi         ( 0 ),
+        outi         ( 100 ),
+        save         ( 0 ),
+        runame ( "heat_results" ),
+        alg  ( "ftcs" ),
+        ic   ( "const(1)" ),
+        lenx      ( Number(1.0) ),
+        alpha     ( Number(0.2) ),
+        dt        ( Number(0.004) ),
+        dx        ( Number(0.1) ),
+        bc0       ( Number(0.0) ),
+        bc1       ( Number(1.0) ),
+        maxt      ( Number(2.0) ),
+        min_change ( Number(1e-8*1e-8) ) {}
+};
+
+Args process_args(int argc, char **argv);
diff --git a/crankn.C b/crankn.C
@@ -4,8 +4,9 @@
 // Modified: 30 May 2009 Author: John Burkardt
 // Modified by Mark C. Miller, July 23, 2017
 static void
-r83_np_fa(int n, Number *a)
+r83_np_fa(Vector &a)
 {
+    int n = a.size()/3;
     int i;
 
     for ( i = 1; i <= n-1; i++ )
@@ -22,16 +23,15 @@ r83_np_fa(int n, Number *a)
     assert( a[1+(n-1)*3] != 0.0 );
 }
 
-void
+Vector
 initialize_crankn(int n,
-    Number alpha, Number dx, Number dt,
-    Number **_cn_Amat)
+    Number alpha, Number dx, Number dt)
 {
     int i;
     Number const w = alpha * dt / dx / dx;
 
     // Build a tri-diagonal matrix
-    Number *cn_Amat = new Number[3*n]();
+    Vector cn_Amat(3*n);
 
     cn_Amat[0+0*3] = 0.0;
     cn_Amat[1+0*3] = 1.0;
@@ -49,10 +49,9 @@ initialize_crankn(int n,
     cn_Amat[2+(n-1)*3] = 0.0;
 
     // Factor the matrix.
-    r83_np_fa(n, cn_Amat);
+    r83_np_fa(cn_Amat);
 
-    // Return the generated matrix
-    *_cn_Amat = cn_Amat;
+    return cn_Amat;
 }
 
 // Licensing: This code is distributed under the GNU LGPL license. 

diff --git a/dufrank.C b/dufrank.C
@@ -11,7 +11,7 @@ update_solution_dufrank(
     Number bc0, Number bc1) // boundary conditions @ x=0 & x=Lx
 {
     Number r = alpha * dt / (dx * dx);
-    Number q = 1 / (1+r);
+    Number q = castNum(1.0) / (castNum(1.0)+r);
 
     // FTCS update algorithm
     for (int i = 1; i < n-1; i++)