Restore ReLAPACK test folder

relapack/test/README.md

@@ -0,0 +1,48 @@
+ReLAPACK Test Suite
+===================
+This test suite compares ReLAPACK's recursive routines with LAPACK's compute
+routines in terms of accuracy:  For each test-case, we execute both ReLAPACK's
+and LAPACK's routine on the same data and consider the numerical difference
+between the two solutions.
+
+This difference is computed as the maximum error across all elements of the
+routine's outputs, where the error for each element is the minimum of the
+absolute error and the relative error (with LAPACK as the reference).  If the
+error is below the error bound configured in `config.h` (default: 1e-5 for
+single precision and 1e-14 for double precision) the test-case is considered as
+passed.
+
+For each routine the test-cases cover a variety of input argument combinations
+to ensure that ReLAPACK's routines match the functionality of LAPACK for all use
+cases.
+
+The matrix size for all experiments (default: 100) can also be specified in
+`config.h`.
+
+
+Implementation
+--------------
+`test.h` provides the framework for our tests:  It provides macros that allow to
+generalize the tests for each operation in one file covering all data-types.
+Such a file is structured as follows:
+
+ * All matrices required by the test-cases are declared globally.  For each
+   matrix, an array of two pointers is declared; one for the matrix copy passed
+   to ReLAPACK and one passed to LAPACK.
+
+ * `tests()` contains the main control flow: it allocates (and later frees) the
+   copies of the globally declared matrices.  It then defines the macro
+   `ROUTINE` to contain the name of the currently tested routine.
+   It then uses the macro `TEST` to perform the test-cases.
+   It receives the arguments of the routine, where matrices of which ReLAPACK
+   and LAPACK receive a copy are index with `i`. (Example: `TEST("L", &n, A[i],
+   &n, info);`)
+
+ * The macro `TEST` first calls `pre()`, which initializes all relevant
+   matrices, then executes the ReLAPACK algorithm on the matrices with `i` = `0`
+   and then the LAPACK counter part with `i` = `1`.  It then calls `post()`,
+   which computes the difference between the results, storing it in `error`.
+   Finally, the error is printed out and compared to the error bound.
+
+If all test-cases pass the error bound test, the program will have a `0` return
+value, otherwise it is `1`, indicating an error.

relapack/test/config.h

@@ -0,0 +1,13 @@
+#ifndef TEST_CONFIG_H
+#define TEST_CONFIG_H
+
+// error bound for single and single complex routines
+#define SINGLE_ERR_BOUND 1e-4
+
+// error bound for double an double complex routines
+#define DOUBLE_ERR_BOUND 1e-13
+
+// size of test matrices
+#define TEST_SIZE 100
+
+#endif /* TEST_CONFIG_H */

relapack/test/lapack.h

@@ -0,0 +1,64 @@
+#ifndef LAPACK_H2
+#define LAPACK_H2
+
+#include "../config.h"
+
+void LAPACK(slauum)(const char *, const int *, float *, const int *, int *);
+void LAPACK(dlauum)(const char *, const int *, double *, const int *, int *);
+void LAPACK(clauum)(const char *, const int *, float *, const int *, int *);
+void LAPACK(zlauum)(const char *, const int *, double *, const int *, int *);
+
+void LAPACK(strtri)(const char *, const char *, const int *, float *, const int *, int *);
+void LAPACK(dtrtri)(const char *, const char *, const int *, double *, const int *, int *);
+void LAPACK(ctrtri)(const char *, const char *, const int *, float *, const int *, int *);
+void LAPACK(ztrtri)(const char *, const char *, const int *, double *, const int *, int *);
+
+void LAPACK(spotrf)(const char *, const int *, float *, const int *, int *);
+void LAPACK(dpotrf)(const char *, const int *, double *, const int *, int *);
+void LAPACK(cpotrf)(const char *, const int *, float *, const int *, int *);
+void LAPACK(zpotrf)(const char *, const int *, double *, const int *, int *);
+
+void LAPACK(spbtrf)(const char *, const int *, const int *, float *, const int *, int *);
+void LAPACK(dpbtrf)(const char *, const int *, const int *, double *, const int *, int *);
+void LAPACK(cpbtrf)(const char *, const int *, const int *, float *, const int *, int *);
+void LAPACK(zpbtrf)(const char *, const int *, const int *, double *, const int *, int *);
+
+void LAPACK(ssytrf)(const char *, const int *, float *, const int *, int *, float *, const int *, int *);
+void LAPACK(dsytrf)(const char *, const int *, double *, const int *, int *, double *, const int *, int *);
+void LAPACK(csytrf)(const char *, const int *, float *, const int *, int *, float *, const int *, int *);
+void LAPACK(chetrf)(const char *, const int *, float *, const int *, int *, float *, const int *, int *);
+void LAPACK(zsytrf)(const char *, const int *, double *, const int *, int *, double *, const int *, int *);
+void LAPACK(zhetrf)(const char *, const int *, double *, const int *, int *, double *, const int *, int *);
+void LAPACK(ssytrf_rook)(const char *, const int *, float *, const int *, int *, float *, const int *, int *);
+void LAPACK(dsytrf_rook)(const char *, const int *, double *, const int *, int *, double *, const int *, int *);
+void LAPACK(csytrf_rook)(const char *, const int *, float *, const int *, int *, float *, const int *, int *);
+void LAPACK(chetrf_rook)(const char *, const int *, float *, const int *, int *, float *, const int *, int *);
+void LAPACK(zsytrf_rook)(const char *, const int *, double *, const int *, int *, double *, const int *, int *);
+void LAPACK(zhetrf_rook)(const char *, const int *, double *, const int *, int *, double *, const int *, int *);
+
+void LAPACK(sgetrf)(const int *, const int *, float *, const int *, int *, int *);
+void LAPACK(dgetrf)(const int *, const int *, double *, const int *, int *, int *);
+void LAPACK(cgetrf)(const int *, const int *, float *, const int *, int *, int *);
+void LAPACK(zgetrf)(const int *, const int *, double *, const int *, int *, int *);
+
+void LAPACK(sgbtrf)(const int *, const int *, const int *, const int *, float *, const int *, int *, int *);
+void LAPACK(dgbtrf)(const int *, const int *, const int *, const int *, double *, const int *, int *, int *);
+void LAPACK(cgbtrf)(const int *, const int *, const int *, const int *, float *, const int *, int *, int *);
+void LAPACK(zgbtrf)(const int *, const int *, const int *, const int *, double *, const int *, int *, int *);
+
+void LAPACK(ssygst)(const int *, const char *, const int *, float *, const int *, const float *, const int *, int *);
+void LAPACK(dsygst)(const int *, const char *, const int *, double *, const int *, const double *, const int *, int *);
+void LAPACK(chegst)(const int *, const char *, const int *, float *, const int *, const float *, const int *, int *);
+void LAPACK(zhegst)(const int *, const char *, const int *, double *, const int *, const double *, const int *, int *);
+
+void LAPACK(strsyl)(const char *, const char *, const int *, const int *, const int *, const float *, const int *, const float *, const int *, float *, const int *, float *, int *);
+void LAPACK(dtrsyl)(const char *, const char *, const int *, const int *, const int *, const double *, const int *, const double *, const int *, double *, const int *, double *, int *);
+void LAPACK(ctrsyl)(const char *, const char *, const int *, const int *, const int *, const float *, const int *, const float *, const int *, float *, const int *, float *, int *);
+void LAPACK(ztrsyl)(const char *, const char *, const int *, const int *, const int *, const double *, const int *, const double *, const int *, double *, const int *, double *, int *);
+
+void LAPACK(stgsyl)(const char *, const int *, const int *, const int *, const float *, const int *, const float *, const int *, float *, const int *, const float *, const int *, const float *, const int *, float *, const int *, float *, float *, float *, const int *, int *, int *);
+void LAPACK(dtgsyl)(const char *, const int *, const int *, const int *, const double *, const int *, const double *, const int *, double *, const int *, const double *, const int *, const double *, const int *, double *, const int *, double *, double *, double *, const int *, int *, int *);
+void LAPACK(ctgsyl)(const char *, const int *, const int *, const int *, const float *, const int *, const float *, const int *, float *, const int *, const float *, const int *, const float *, const int *, float *, const int *, float *, float *, float *, const int *, int *, int *);
+void LAPACK(ztgsyl)(const char *, const int *, const int *, const int *, const double *, const int *, const double *, const int *, double *, const int *, const double *, const int *, const double *, const int *, double *, const int *, double *, double *, double *, const int *, int *, int *);
+
+#endif /*  LAPACK_H2 */

relapack/test/test.h

@@ -0,0 +1,136 @@
+#ifndef TEST_H
+#define TEST_H
+
+#include "../config.h"
+#include "config.h"
+
+#if BLAS_UNDERSCORE
+#define BLAS(routine) routine ## _
+#else
+#define BLAS(routine) routine
+#endif
+
+#if LAPACK_UNDERSCORE
+#define LAPACK(routine) routine ## _
+#else
+#define LAPACK(routine) routine
+#endif
+
+#include "../inc/relapack.h"
+#include "lapack.h"
+#include "util.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+// some name mangling macros
+#define CAT(A, B) A ## B
+#define XCAT(A, B) CAT(A, B)
+#define XLAPACK(X) LAPACK(X)
+#define XRELAPACK(X) XCAT(RELAPACK_, X)
+#define STR(X) #X
+#define XSTR(X) STR(X)
+
+// default setup and error computation names: pre() and post()
+#define PRE pre
+#define POST post
+
+// TEST macro:
+// run setup (pre()), ReLAPACK routine (i = 0), LAPACK routine (i = 1), compute
+// error (post()), check error bound, and  print setup and error
+#define TEST(...) \
+    PRE(); \
+    i = 0; \
+    XRELAPACK(ROUTINE)(__VA_ARGS__); \
+    i = 1; \
+    XLAPACK(ROUTINE)(__VA_ARGS__); \
+    POST(); \
+    fail |= error > ERR_BOUND; \
+    printf("%s(%s)\t%g\n", XSTR(ROUTINE), #__VA_ARGS__, error);
+
+// generalized datatype treatment: DT_PREFIX determines the type s, d, c, or z
+#define XPREF(A) XCAT(DT_PREFIX, A)
+
+// matrix generation and error computation routines
+#define x2matgen XPREF(2matgen)
+#define x2vecerr XPREF(2vecerr)
+
+// error bounds
+#define ERR_BOUND XPREF(ERR_BOUND_)
+#define sERR_BOUND_ SINGLE_ERR_BOUND
+#define dERR_BOUND_ DOUBLE_ERR_BOUND
+#define cERR_BOUND_ SINGLE_ERR_BOUND
+#define zERR_BOUND_ DOUBLE_ERR_BOUND
+
+// C datatypes
+#define datatype XPREF(datatype_)
+#define sdatatype_ float
+#define ddatatype_ double
+#define cdatatype_ float
+#define zdatatype_ double
+
+// number of C datatype elements per element
+#define x1 XPREF(DT_MULT)
+#define sDT_MULT 1
+#define dDT_MULT 1
+#define cDT_MULT 2
+#define zDT_MULT 2
+
+// typed allocations
+#define xmalloc XPREF(malloc)
+#define imalloc(S) malloc((S) * sizeof(int))
+#define smalloc(S) malloc((S) * sizeof(float))
+#define dmalloc(S) malloc((S) * sizeof(double))
+#define cmalloc(S) malloc((S) * 2 * sizeof(float))
+#define zmalloc(S) malloc((S) * 2 * sizeof(double))
+
+// transpositions
+#define xCTRANS XPREF(CTRANS)
+#define sCTRANS "T"
+#define dCTRANS "T"
+#define cCTRANS "C"
+#define zCTRANS "C"
+
+// some constants
+#define MONE XPREF(MONE)
+const float  sMONE[] = { -1. };
+const double dMONE[] = { -1. };
+const float  cMONE[] = { -1., 0. };
+const double zMONE[] = { -1., 0. };
+
+#define ZERO XPREF(ZERO)
+const float  sZERO[] = { 0. };
+const double dZERO[] = { 0. };
+const float  cZERO[] = { 0., 0. };
+const double zZERO[] = { 0., 0. };
+
+#define ONE  XPREF(ONE)
+const float  sONE[]  = { 1. };
+const double dONE[]  = { 1. };
+const float  cONE[]  = { 1., 0. };
+const double zONE[]  = { 1., 0. };
+
+const int iMONE[]  = { -1 };
+const int iZERO[]  = { 0 };
+const int iONE[]   = { 1 };
+const int iTWO[]   = { 2 };
+const int iTHREE[] = { 3 };
+const int iFOUR[]  = { 4 };
+
+void tests();
+
+// global variables (used in tests(), pre(), and post())
+int i, n, n2, fail;
+double error;
+
+int main(int argc, char* argv[]) {
+    n = TEST_SIZE;
+    n2 = (3 * n) / 4;
+    fail = 0;
+
+    tests();
+
+    return fail;
+}
+
+#endif /* TEST_H */

relapack/test/util.c

@@ -0,0 +1,116 @@
+#include "util.h"
+#include <stdlib.h>
+#include <time.h>
+#include <math.h>
+
+#define MAX(a, b) ((a) > (b) ? (a) : (b))
+#define MIN(a, b) ((a) < (b) ? (a) : (b))
+
+///////////////////////
+// matrix generation //
+///////////////////////
+// Each routine x2matgen is passed the size (m, n) of the desired matrix and
+// geneartes two copies of such a matrix in in its output arguments A and B.
+// The generated matrices is filled with random entries in [0, 1[ (+i*[0, 1[ in
+// the complex case).  Then m is added to the diagonal; this is numerically
+// favorable for routines working with triangular and symmetric matrices.  For
+// the same reason the imaginary part of the diagonal is set to 0.
+
+void s2matgen(const int m, const int n, float *A, float *B) {
+    srand(time(NULL) + (size_t) A);
+    int i, j;
+    for (i = 0; i < m; i++)
+        for (j = 0; j < n; j++)
+            A[i + m * j] = B[i + m * j] = (float) rand() / RAND_MAX + m * (i == j);
+}
+
+void d2matgen(const int m, const int n, double *A, double *B) {
+    srand(time(NULL) + (size_t) A);
+    int i, j;
+    for (i = 0; i < m; i++)
+        for (j = 0; j < n; j++)
+            A[i + m * j] = B[i + m * j] = (double) rand() / RAND_MAX + m * (i == j);
+}
+
+void c2matgen(const int m, const int n, float *A, float *B) {
+    srand(time(NULL) + (size_t) A);
+    int i, j;
+    for (i = 0; i < m; i++)
+        for (j = 0; j < n; j++) {
+            A[2* (i + m * j)]     = B[2 * (i + m * j)]     = (float) rand() / RAND_MAX + m * (i == j);
+            A[2* (i + m * j) + 1] = B[2 * (i + m * j) + 1] = ((float) rand() / RAND_MAX) * (i != j);
+        }
+}
+
+void z2matgen(const int m, const int n, double *A, double *B) {
+    srand(time(NULL) + (size_t) A);
+    int i, j;
+    for (i = 0; i < m; i++)
+        for (j = 0; j < n; j++) {
+            A[2* (i + m * j)]     = B[2 * (i + m * j)]     = (double) rand() / RAND_MAX + m * (i == j);
+            A[2* (i + m * j) + 1] = B[2 * (i + m * j) + 1] = ((double) rand() / RAND_MAX) * (i != j);
+        }
+}
+
+////////////////////////
+// error computations //
+////////////////////////
+// Each routine x2vecerrr is passed a vector lengh n and two vectors x and y.
+// It returns the maximum of the element-wise error between these two vectors.
+// This error is the minimum of the absolute difference and the relative
+// differene with respect to y.
+
+double i2vecerr(const int n, const int *x, const int *y) {
+    double error = 0;
+    int i;
+    for (i = 0; i < n; i++) {
+        double nom = abs(x[i] - y[i]);
+        double den = abs(y[i]);
+        error = MAX(error, (den > 0) ? MIN(nom, nom / den) : nom);
+    }
+    return error;
+}
+
+double s2vecerr(const int n, const float *x, const float *y) {
+    float error = 0;
+    int i;
+    for (i = 0; i < n; i++) {
+        double nom = fabs((double) x[i] - y[i]);
+        double den = fabs(y[i]);
+        error = MAX(error, (den > 0) ? MIN(nom, nom / den) : nom);
+    }
+    return error;
+}
+
+double d2vecerr(const int n, const double *x, const double *y) {
+    double error = 0;
+    int i;
+    for (i = 0; i < n; i++) {
+        double nom = fabs(x[i] - y[i]);
+        double den = fabs(y[i]);
+        error = MAX(error, (den > 0) ? MIN(nom, nom / den) : nom);
+    }
+    return error;
+}
+
+double c2vecerr(const int n, const float *x, const float *y) {
+    double error = 0;
+    int i;
+    for (i = 0; i < n; i++) {
+        double nom = sqrt(((double) x[2 * i] - y[2 * i]) * ((double) x[2 * i] - y[2 * i]) + ((double) x[2 * i + 1] - y[2 * i + 1]) * ((double) x[2 * i + 1] - y[2 * i + 1]));
+        double den = sqrt((double) y[2 * i] * y[2 * i] + (double) y[2 * i + 1] * y[2 * i + 1]);
+        error = MAX(error, (den > 0) ? MIN(nom, nom / den) : nom);
+    }
+    return error;
+}
+
+double z2vecerr(const int n, const double *x, const double *y) {
+    double error = 0;
+    int i;
+    for (i = 0; i < n; i++) {
+        double nom = sqrt((x[2 * i] - y[2 * i]) * (x[2 * i] - y[2 * i]) + (x[2 * i + 1] - y[2 * i + 1]) * (x[2 * i + 1] - y[2 * i + 1]));
+        double den = sqrt(y[2 * i] * y[2 * i] + y[2 * i + 1] * y[2 * i + 1]);
+        error = MAX(error, (den > 0) ? MIN(nom, nom / den) : nom);
+    }
+    return error;
+}

relapack/test/util.h

@@ -0,0 +1,15 @@
+#ifndef TEST_UTIL_H
+#define TEST_UTIL_H
+
+void s2matgen(int, int, float *, float *);
+void d2matgen(int, int, double *, double *);
+void c2matgen(int, int, float *, float *);
+void z2matgen(int, int, double *, double *);
+
+double i2vecerr(int, const int *, const int *);
+double s2vecerr(int, const float *, const float *);
+double d2vecerr(int, const double *, const double *);
+double c2vecerr(int, const float *, const float *);
+double z2vecerr(int, const double *, const double *);
+
+#endif /* TEST_UTIL_H */