Initial commit for accuracy lapack.

9 files changed, 431 insertions, 6 deletions

diff --git a/accuracy/TestAccuracy.hpp b/accuracy/TestAccuracy.hpp
new file mode 100644
index 0000000..b0a6cc1
--- /dev/null
+++ b/accuracy/TestAccuracy.hpp
@@ -0,0 +1,75 @@
+#ifndef TESTACCURACY_HPP
+#define TESTACCURACY_HPP
+
+#include <utility>
+#include <cmath>
+#include <vector>
+#include <iostream>
+
+#include "utilities/Timer.hpp"
+
+std::vector<int> getLogSizes(
+        const int& minsize = 4,
+        const int& maxsize = 3000,
+        const unsigned& N = 100u
+    )
+{
+    std::vector<int> sizes(N);
+    const double amin = std::log(minsize), amax = std::log(maxsize);
+    const double k = (amax-amin) / (N-1u);
+    double a;
+
+    for (unsigned i = 0; i < N; ++i) {
+        a = amin + i*k;
+        a = std::floor(std::exp(a)+.5);
+        sizes[i] = a;
+    }
+
+    sizes.front() = minsize;
+    sizes.back() = maxsize;
+
+    return sizes;
+}
+
+
+template<
+    template<class> class Action,
+    typename value_t
+>
+std::pair<value_t, value_t> testAccuracy(
+        const int& minsize = 4,
+        const int& maxsize = 1000,
+        const unsigned& N = 100u)
+{
+    std::vector<int> sizes = getLogSizes(minsize, maxsize, N);
+    Timer t;
+
+    for (std::vector<int>::const_iterator i = sizes.begin(), e = sizes.end();
+         i != e; ++i) {
+
+        int size = *i;
+        std::cout << " -- Size: " << size << "  -  " << std::flush;
+
+        t.start();
+        int N = 0;
+        double e, emean = 0., e2mean = 0.;
+        do {
+            Action<value_t> action(size, N);
+            action.compute();
+            e = action.getResidual();
+            emean += e;
+            e2mean += e*e;
+        } while(++N < 100 && t.elapsed() < 1. || N < 4);
+
+        std::cout << N << " samples  -  " << t.elapsed() << " seconds  -  ";
+
+        emean /= N;
+        e2mean /= N;
+        e2mean = std::sqrt(e2mean - emean*emean);
+
+        std::cout << emean << " +/- " << e2mean << std::endl;
+
+    }
+}
+
+#endif // TESTACCURACY_HPP
diff --git a/accuracy/actions/ActionGESV.hpp b/accuracy/actions/ActionGESV.hpp
new file mode 100644
index 0000000..454e5f6
--- /dev/null
+++ b/accuracy/actions/ActionGESV.hpp
@@ -0,0 +1,53 @@
+#ifndef ACTIONGESV_HPP
+#define ACTIONGESV_HPP
+
+#include <vector>
+#include <iostream>
+#include "LinearCongruential.hpp"
+
+extern "C" {
+    void dgesv_(const int*, const int*, double*, const int*, double*,
+                double*, const int*, int*);
+
+    void dgemv_(const char*, const int*, const int*, const double*,
+                const double*, const int*, const double*, const int*,
+                const double*, double*, const int*);
+
+    double dnrm2_(const int*, const double*, const int*);
+}
+
+template<typename value_t = double>
+class ActionGESV {
+    typedef std::vector<value_t> storage_t;
+    typedef LinearCongruential<> rangen_t;
+
+public:
+    ActionGESV(const int& size, const unsigned& seed=0)
+        : size(size), rg(seed), A(rg.fillVector<>(size*size, 0.)), Acopy(A),
+          x(rg.fillVector(size, 0.)), b(x), bcopy(b)
+    { }
+
+    void compute() {
+        const int ONE = 1;
+        int info;
+        std::vector<value_t> ipiv(size);
+        dgesv_(&size, &ONE, &A[0], &size, &ipiv[0], &x[0], &size, &info);
+        if (info != 0)
+            std::cerr << "Info: " << info << "\n";
+    }
+
+    double getResidual() {
+        const double alpha = -1., beta = 1.;
+        const int ONE = 1;
+        dgemv_("N", &size, &size, &alpha, &Acopy[0], &size, &x[0], &ONE, &beta,
+                &b[0], &ONE);
+        return dnrm2_(&size, &b[0], &ONE)/dnrm2_(&size, &bcopy[0], &ONE);
+    }
+
+//private:
+    const int size;
+    rangen_t rg;
+    storage_t A, Acopy, x, b, bcopy;
+};
+
+#endif //  ACTIONGESV_HPP
diff --git a/accuracy/libraries/LAPACK/main.cpp b/accuracy/libraries/LAPACK/main.cpp
new file mode 100644
index 0000000..4dafab4
--- /dev/null
+++ b/accuracy/libraries/LAPACK/main.cpp
@@ -0,0 +1,23 @@
+#include <iostream>
+#include <string>
+
+#include "TestAccuracy.hpp"
+
+// Actions
+#include "ActionGESV.hpp"
+
+using namespace std;
+
+int main(int argc, char **argv)
+{
+    bool do_gesv = false;
+
+    // Parse input
+    for (int i = 1; i < argc; ++i) {
+        std::string arg = argv[i];
+        if (arg == "general_solve") do_gesv = true;
+    }
+
+    if (do_gesv)
+        testAccuracy<ActionGESV, double>(4, 3000, 20);
+}
diff --git a/accuracy/utilities/LinearCongruential.hpp b/accuracy/utilities/LinearCongruential.hpp
new file mode 100644
index 0000000..a735f77
--- /dev/null
+++ b/accuracy/utilities/LinearCongruential.hpp
@@ -0,0 +1,77 @@
+#ifndef LINEARCONGRUENTIAL_HPP
+#define LINEARCONGRUENTIAL_HPP
+
+#include "stdint.h"
+#include <vector>
+#include <limits>
+
+template<
+    typename value_t = uint32_t,
+    value_t a = 1664525u,
+    value_t c = 1013904223u
+>
+class LinearCongruential {
+    typedef typename std::vector<value_t> buf_t;
+    typedef typename buf_t::iterator iterator_t;
+    typedef typename buf_t::const_iterator constiterator_t;
+
+public:
+    LinearCongruential(const value_t& seed, const int bufsize = 4*1024)
+        : maxvalue(std::numeric_limits<value_t>::max())
+    {
+        buffer.resize(bufsize);
+        fillBuffer(seed);
+    }
+
+    value_t operator()()
+    {
+        if (j == buffer.end())
+            fillBuffer(buffer.back());
+        return *(j++);
+    }
+
+    template<typename result_t>
+    void fillVector (
+            result_t *vector,
+            const int& size,
+            const result_t& min = 0.,
+            const result_t& max = 1.
+        )
+    {
+        const result_t maxvalue = this->maxvalue;
+        const result_t diff = max-min;
+
+        for (const result_t *const endp = vector+size; vector!=endp; ++vector)
+            *vector = diff * (this->operator()() / maxvalue) + min;
+    }
+
+
+    template<typename result_t>
+    std::vector<result_t> fillVector (
+            const int& size,
+            const result_t& min = 0.,
+            const result_t& max = 1.
+        )
+    {
+        std::vector<result_t> result(size);
+        fillVector(&result[0], size, min, max);
+        return result;
+    }
+
+private:
+    buf_t buffer;
+    constiterator_t j;
+
+    const value_t maxvalue;
+
+    void fillBuffer(int seed)
+    {
+        for (iterator_t i = buffer.begin(), e = buffer.end(); i!=e; ++i) {
+            *i = ((seed *= a) += c);
+        }
+
+        j = buffer.begin();
+    }
+};
+
+#endif // LINEARCONGRUENTIAL_HPP
diff --git a/accuracy/utilities/Timer.hpp b/accuracy/utilities/Timer.hpp
new file mode 100644
index 0000000..6cfeac1
--- /dev/null
+++ b/accuracy/utilities/Timer.hpp
@@ -0,0 +1,43 @@
+#include <sys/time.h>
+
+class Timer
+{
+public:
+    Timer() : ZERO(0.) { }
+
+    void start()
+    {
+        t = walltime(&ZERO);
+    }
+
+    double elapsed()
+    {
+        return walltime(&t);
+    }
+
+private:
+    double t;
+
+    const double ZERO;
+
+    static double walltime(const double *t0) {
+        double mic;
+        double time;
+        double mega = 0.000001;
+        struct timeval tp;
+        struct timezone tzp;
+        static long base_sec = 0;
+        static long base_usec = 0;
+
+        (void) gettimeofday(&tp, &tzp);
+        if (base_sec == 0){
+            base_sec  = tp.tv_sec;
+            base_usec = tp.tv_usec;
+        }
+
+        time = (double)(tp.tv_sec - base_sec);
+        mic  = (double)(tp.tv_usec - base_usec);
+        time = (time + mic * mega) - *t0;
+        return(time);
+    }
+};
diff --git a/numbench/benchconfig.py b/numbench/benchconfig.py
index 56ef304..bb58651 100644
--- a/numbench/benchconfig.py
+++ b/numbench/benchconfig.py
@@ -29,6 +29,7 @@ imageformat = None
 curdir = None
 scriptdir = None
 btldir = None
+accudir = None
 libdir = None
 
 # Storage directories
@@ -88,17 +89,16 @@ def parseArguments():
 
 
 def setDirs():
-    global curdir, scriptdir, btldir, libdir, basedir
+    global curdir, scriptdir, btldir, accudir, libdir, basedir
     global testsdir, rootsdir, pkgsdir, reportdir, logdir
 
 
     # Script directories
     curdir = os.path.abspath('.')
     scriptdir = os.path.dirname(os.path.realpath(__file__))
-    if os.environ.has_key('BTLDIR'):
-        btldir = os.environ['BTLDIR']
-    else:
-        btldir = '/usr/include/numbench/btl'
+
+    btldir = os.environ.get('BTLDIR', '/usr/include/numbench/btl')
+    accudir = os.environ.get('ACCUDIR', '/usr/include/numbench/accuracy')
 
     # Library directory (lib vs. lib32 vs. lib64)
     libdir = sp.Popen \
diff --git a/numbench/modules/__init__.py b/numbench/modules/__init__.py
index 9696252..cb61b6a 100644
--- a/numbench/modules/__init__.py
+++ b/numbench/modules/__init__.py
@@ -25,7 +25,8 @@ modnames = [
   'lapack',
   'lapacke',
   'scalapack',
-  'fftw'
+  'fftw',
+  'lapack_accuracy'
 ]
 
 
diff --git a/numbench/modules/internal/accuracyBase.py b/numbench/modules/internal/accuracyBase.py
new file mode 100644
index 0000000..e6beeb1
--- /dev/null
+++ b/numbench/modules/internal/accuracyBase.py
@@ -0,0 +1,125 @@
+from os.path import dirname, join as pjoin
+import os, subprocess as sp
+
+from numbench import benchchildren, benchconfig as cfg
+from numbench.utils import alternatives as alt, benchutils as bu
+from numbench.benchprint import Print
+
+def compileExe(test, libname, implementation):
+    src = pjoin(cfg.accudir, 'libraries', libname.upper(), 'main.cpp')
+    exe = pjoin(test['testdir'], implementation, 'test')
+
+    # Compiler switches
+    includes = ['-I' + pjoin(cfg.accudir, i) for i in
+                ('', 'actions', 'utilities')]
+    libraries = ['-lm']
+
+    # Flags (also update compile and run environments)
+    flags = alt.getFlags(test, libname, implementation)
+    libenv = ""
+    for i in [j[2:] for j in flags if j.startswith('-L')]:
+        libenv += i + ":"
+
+    libenvc = libenv + test['compileenv'].get('LIBRARY_PATH', '')
+    test['compileenv']['LIBRARY_PATH'] = libenvc
+
+    libenvr = libenv + test['runenv'].get('LD_LIBRARY_PATH', '')
+    print "\n\nLIBENV: ", libenvr, "\n\n"
+    test['runenv']['LD_LIBRARY_PATH'] = libenvr
+
+    # Set compile-time environment
+    compileenv = test['compileenv'].copy()
+
+
+    # Set C++ compiler
+    cxx = compileenv.get('CXX', '') or \
+          bu.run_cmd(['portageq', 'envvar', 'CXX']).strip() or '/usr/bin/g++'
+
+
+    # Form command-line arguments
+    args = [cxx, src, '-o', exe] + includes + libraries + flags
+
+    # Open logfile and write environment
+    logdir = pjoin(test['logdir'], implementation)
+    bu.mkdir(logdir)
+    logfile = pjoin(logdir, "accuracyCompile.log")
+    logfs = file(logfile, 'w')
+    logfs.write('\n'.join([n + '=' + v for n, v in compileenv.items()]))
+    logfs.write(3 * '\n' + ' '.join(args) + 3 * '\n')
+    logfs.flush()
+
+    # Execute compilation
+    bu.mkdir(dirname(exe))
+    proc = sp.Popen(args, stdout=logfs, stderr=sp.STDOUT, env=compileenv)
+    proc.wait()
+    logfs.flush()
+    retcode = proc.returncode
+    if retcode == 0:
+        logfs.write("\n\n<<< Compilation terminated successfully >>>")
+    else:
+        logfs.write("\n\n<<< Compilation failed >>>")
+
+    # Close, return
+    logfs.close()
+    return retcode, exe
+
+
+def runExe(test, implementation, exe, args):
+    logdir = pjoin(test['logdir'], implementation)
+
+    # Check linking
+    logfs = file(pjoin(logdir, 'accuracyLinking.log'), 'w')
+    sp.Popen(['ldd', '-v', exe], stdout=logfs, env=test['runenv']).wait()
+    logfs.close()
+
+    # Prepare arguments
+    args = (exe,) + tuple(args)
+
+    # Open log
+    logfs = file(pjoin(logdir, "accuracyRun.log"), 'w')
+    logfs.write('\n'.join([n + '=' + v for n, v in test['runenv'].items()]))
+    logfs.write(3 * '\n' + ' '.join(args) + 3 * '\n')
+    logfs.flush()
+
+    # Error log
+    errfname = pjoin(logdir, "accuracyRun.err")
+    errfs = file(errfname, 'w')
+
+    # Open pipe
+    try:
+        testdir = pjoin(test['testdir'], implementation)
+        proc = sp.Popen(args, bufsize=1, stdout=sp.PIPE, stderr=errfs, \
+                        env=test['runenv'], cwd=testdir)
+        benchchildren.append(proc)
+    except OSError:
+        Print('Execution failed to start')
+        Print('Command line: ' + ' '.join(args))
+        return -1, None
+
+    result = {}
+
+    # Interpret output
+    Print('Begin execution')
+    # TODO: interpret output, store results,...
+    proc.wait()
+    Print("Execution finished with return code " + str(proc.returncode))
+
+    # Close logs
+    logfs.close()
+    errp = errfs.tell()
+    errfs.close()
+    if errp == 0:
+        os.unlink(errfname)
+
+
+    # Close, return
+    logfs.close()
+    return proc.returncode, result
+
+
+def runTest(self, test, implementation):
+    retcode, exe = compileExe(test, self.libname, implementation)
+    runExe(test, implementation, exe, self.tests)
+
+def reportConf(*args):
+    return {'type':'plot', 'xlabel':'size', 'ylabel':'Error'}
diff --git a/numbench/modules/lapack_accuracy.py b/numbench/modules/lapack_accuracy.py
new file mode 100644
index 0000000..53f82cc
--- /dev/null
+++ b/numbench/modules/lapack_accuracy.py
@@ -0,0 +1,28 @@
+#=====================================================
+# Copyright (C) 2012 Andrea Arteaga <andyspiros@gmail.com>
+#=====================================================
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+#
+from internal import accuracyBase, lapackBase
+
+class Module:
+    libname = 'lapack'
+    descr = 'Accuracy test module for LAPACK implementations'
+
+    __init__ = lapackBase.init
+    getImplementations = lapackBase.getImplementations
+    runTest = accuracyBase.runTest
+    getTests = lapackBase.getTests
+    reportConf = accuracyBase.reportConf