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 */

relapack/test/xgbtrf.c

@@ -0,0 +1,43 @@
+#include "test.h"
+
+datatype *A[2];
+int *ipiv[2], info;
+int kl, ku, ld;
+
+void pre() {
+    int i;
+    x2matgen(ld, n, A[0], A[1]);
+    for (i = 0; i < n; i++) {
+        // set diagonal
+        A[0][x1 * (i + ld * i)] =
+        A[1][x1 * (i + ld * i)] = (datatype) rand() / RAND_MAX;
+    }
+    memset(ipiv[0], 0, n * sizeof(int));
+    memset(ipiv[1], 0, n * sizeof(int));
+}
+
+void post() {
+    error = x2vecerr(ld * n, A[0], A[1]) + i2vecerr(n, ipiv[0], ipiv[1]);
+}
+
+void tests() {
+    kl = n - 10;
+    ku = n;
+    ld = 2 * kl + ku + 1;
+
+    A[0] = xmalloc(ld * n);
+    A[1] = xmalloc(ld * n);
+    ipiv[0] = imalloc(n);
+    ipiv[1] = imalloc(n);
+
+    #define ROUTINE XPREF(gbtrf)
+
+    TEST(&n, &n, &kl, &ku, A[i], &ld, ipiv[i], &info);
+    TEST(&n, &n2, &kl, &ku, A[i], &ld, ipiv[i], &info);
+    TEST(&n2, &n, &kl, &ku, A[i], &ld, ipiv[i], &info);
+
+    free(A[0]);
+    free(A[1]);
+    free(ipiv[0]);
+    free(ipiv[1]);
+}

relapack/test/xgemmt.c

@@ -0,0 +1,65 @@
+#include "test.h"
+
+datatype *A[2], *B[2], *C[2], *Ctmp;
+int info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+    x2matgen(n, n, B[0], B[1]);
+    x2matgen(n, n, C[0], C[1]);
+}
+
+void post() {
+    error = x2vecerr(n * n, C[0], C[1]);
+}
+
+#define ROUTINE XPREF(gemmt)
+
+#define xlacpy XPREF(LAPACK(lacpy))
+#define xgemm XPREF(BLAS(gemm))
+
+extern void xlacpy(const char *, const int *, const int *, const datatype *, const int *, datatype *, const int *);
+extern void xgemm(const char *, const char *, const int *, const int *, const int *, const datatype *, const datatype *, const int *, const datatype *, const int *, const datatype *, const datatype *, const int*);
+
+void XLAPACK(ROUTINE)(
+    const char *uplo, const char *transA, const char *transB,
+    const int *n, const int *k,
+    const datatype *alpha, const datatype *A, const int *ldA,
+    const datatype *B, const int *ldB,
+    const datatype *beta, datatype *C, const int *ldC
+) {
+    xlacpy(uplo, n, n, C, ldC, Ctmp, n);
+    xgemm(transA, transB, n, n, k, alpha, A, ldA, B, ldB, beta, Ctmp, n);
+    xlacpy(uplo, n, n, Ctmp, ldC, C, n);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+    B[0] = xmalloc(n * n);
+    B[1] = xmalloc(n * n);
+    C[0] = xmalloc(n * n);
+    C[1] = xmalloc(n * n);
+    Ctmp = xmalloc(n * n);
+
+    TEST("L", "N", "N", &n, &n, ONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("L", "N", "N", &n, &n, ONE, A[i], &n, B[i], &n, MONE, C[i], &n);
+    TEST("L", "N", "N", &n, &n, MONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("L", "N", "T", &n, &n, ONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("L", "T", "N", &n, &n, ONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("L", "N", "N", &n, &n2, ONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("U", "N", "N", &n, &n, ONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("U", "N", "N", &n, &n, ONE, A[i], &n, B[i], &n, MONE, C[i], &n);
+    TEST("U", "N", "N", &n, &n, MONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("U", "N", "T", &n, &n, ONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("U", "T", "N", &n, &n, ONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+    TEST("U", "N", "N", &n, &n2, ONE, A[i], &n, B[i], &n, ONE, C[i], &n);
+
+    free(A[0]);
+    free(A[1]);
+    free(B[0]);
+    free(B[1]);
+    free(C[0]);
+    free(C[1]);
+    free(Ctmp);
+}

relapack/test/xgetrf.c

@@ -0,0 +1,32 @@
+#include "test.h"
+
+datatype *A[2];
+int *ipiv[2], info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+    memset(ipiv[0], 0, n * sizeof(int));
+    memset(ipiv[1], 0, n * sizeof(int));
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]) + i2vecerr(n, ipiv[0], ipiv[1]);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+    ipiv[0] = imalloc(n);
+    ipiv[1] = imalloc(n);
+
+    #define ROUTINE XPREF(getrf)
+
+    TEST(&n, &n, A[i], &n, ipiv[i], &info);
+    TEST(&n, &n2, A[i], &n, ipiv[i], &info);
+    TEST(&n2, &n, A[i], &n, ipiv[i], &info);
+
+    free(A[0]);
+    free(A[1]);
+    free(ipiv[0]);
+    free(ipiv[1]);
+}

relapack/test/xhegst.c

@@ -0,0 +1,32 @@
+#include "test.h"
+
+datatype *A[2], *B[2];
+int info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+    x2matgen(n, n, B[0], B[1]);
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+    B[0] = xmalloc(n * n);
+    B[1] = xmalloc(n * n);
+
+    #define ROUTINE XPREF(hegst)
+
+    TEST(iONE, "L", &n, A[i], &n, B[i], &n, &info);
+    TEST(iONE, "U", &n, A[i], &n, B[i], &n, &info);
+    TEST(iTWO, "L", &n, A[i], &n, B[i], &n, &info);
+    TEST(iTWO, "U", &n, A[i], &n, B[i], &n, &info);
+
+    free(A[0]);
+    free(A[1]);
+    free(B[0]);
+    free(B[1]);
+}

relapack/test/xhetrf.c

@@ -0,0 +1,40 @@
+#include "test.h"
+
+datatype *A[2], *Work;
+int *ipiv[2], info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+    memset(ipiv[0], 0, n * sizeof(int));
+    memset(ipiv[1], 0, n * sizeof(int));
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]) + i2vecerr(n, ipiv[0], ipiv[1]);
+}
+
+void tests() {
+    const int lWork = n * n;
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+    ipiv[0] = imalloc(n);
+    ipiv[1] = imalloc(n);
+    Work = xmalloc(lWork);
+
+    #define ROUTINE XPREF(hetrf)
+
+    TEST("L", &n, A[i], &n, ipiv[i], Work, &lWork, &info);
+    TEST("U", &n, A[i], &n, ipiv[i], Work, &lWork, &info);
+
+    #undef ROUTINE
+    #define ROUTINE XPREF(hetrf_rook)
+
+    TEST("L", &n, A[i], &n, ipiv[i], Work, &lWork, &info);
+    TEST("U", &n, A[i], &n, ipiv[i], Work, &lWork, &info);
+
+    free(A[0]);
+    free(A[1]);
+    free(ipiv[0]);
+    free(ipiv[1]);
+    free(Work);
+}

relapack/test/xlauum.c

@@ -0,0 +1,25 @@
+#include "test.h"
+
+datatype *A[2];
+int info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+
+    #define ROUTINE XPREF(lauum)
+
+    TEST("L", &n, A[i], &n, &info);
+    TEST("U", &n, A[i], &n, &info);
+
+    free(A[0]);
+    free(A[1]);
+}

relapack/test/xpbtrf.c

@@ -0,0 +1,40 @@
+#include "test.h"
+
+datatype *A[2];
+int info[2];
+int n;
+
+void pre() {
+    int i;
+    x2matgen(n, n, A[0], A[1]);
+    for (i = 0; i < n; i++) {
+        // set diagonal
+        A[0][x1 * (i + n * i)] =
+        A[1][x1 * (i + n * i)] = (datatype) rand() / RAND_MAX;
+        // set first row
+        A[0][x1 * (n * i)] =
+        A[1][x1 * (n * i)] = (datatype) rand() / RAND_MAX + n;
+    }
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+
+    #define ROUTINE XPREF(pbtrf)
+
+    const int
+        kd1 = n / 4,
+        kd2 = n * 3 / 4;
+    TEST("L", &n, &kd1, A[i], &n, &info[i]);
+    TEST("L", &n, &kd2, A[i], &n, &info[i]);
+    TEST("U", &n, &kd1, A[i] - x1 * kd1, &n, &info[i]);
+    TEST("U", &n, &kd2, A[i] - x1 * kd2, &n, &info[i]);
+
+    free(A[0]);
+    free(A[1]);
+}

relapack/test/xpotrf.c

@@ -0,0 +1,25 @@
+#include "test.h"
+
+datatype *A[2];
+int info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+
+    #define ROUTINE XPREF(potrf)
+
+    TEST("L", &n, A[i], &n, &info);
+    TEST("U", &n, A[i], &n, &info);
+
+    free(A[0]);
+    free(A[1]);
+}

relapack/test/xsygst.c

@@ -0,0 +1,32 @@
+#include "test.h"
+
+datatype *A[2], *B[2];
+int info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+    x2matgen(n, n, B[0], B[1]);
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+    B[0] = xmalloc(n * n);
+    B[1] = xmalloc(n * n);
+
+    #define ROUTINE XPREF(sygst)
+
+    TEST(iONE, "L", &n, A[i], &n, B[i], &n, &info);
+    TEST(iONE, "U", &n, A[i], &n, B[i], &n, &info);
+    TEST(iTWO, "L", &n, A[i], &n, B[i], &n, &info);
+    TEST(iTWO, "U", &n, A[i], &n, B[i], &n, &info);
+
+    free(A[0]);
+    free(A[1]);
+    free(B[0]);
+    free(B[1]);
+}

relapack/test/xsytrf.c

@@ -0,0 +1,40 @@
+#include "test.h"
+
+datatype *A[2], *Work;
+int *ipiv[2], info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+    memset(ipiv[0], 0, n * sizeof(int));
+    memset(ipiv[1], 0, n * sizeof(int));
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]) + i2vecerr(n, ipiv[0], ipiv[1]);
+}
+
+void tests() {
+    const int lWork = n * n;
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+    ipiv[0] = imalloc(n);
+    ipiv[1] = imalloc(n);
+    Work = xmalloc(lWork);
+
+    #define ROUTINE XPREF(sytrf)
+
+    TEST("L", &n, A[i], &n, ipiv[i], Work, &lWork, &info);
+    TEST("U", &n, A[i], &n, ipiv[i], Work, &lWork, &info);
+
+    #undef ROUTINE
+    #define ROUTINE XPREF(sytrf_rook)
+
+    TEST("L", &n, A[i], &n, ipiv[i], Work, &lWork, &info);
+    TEST("U", &n, A[i], &n, ipiv[i], Work, &lWork, &info);
+
+    free(A[0]);
+    free(A[1]);
+    free(ipiv[0]);
+    free(ipiv[1]);
+    free(Work);
+}

relapack/test/xtgsyl.c

@@ -0,0 +1,94 @@
+#include "test.h"
+
+datatype *A[2], *B[2], *C[2], *D[2], *E[2], *F[2], *Work, scale[2], dif[2];
+int *iWork, lWork, info;
+
+#define xlascl XPREF(LAPACK(lascl))
+void xlascl(const char *, const int *, const int *, const datatype *, const
+        datatype *, const int *, const int *, datatype *, const int *, int *);
+
+#define xscal XPREF(LAPACK(scal))
+void xscal(const int *, const datatype *, datatype *, const int *);
+
+void pre() {
+    int i;
+
+    x2matgen(n, n, A[0], A[1]);
+    x2matgen(n, n, B[0], B[1]);
+    x2matgen(n, n, C[0], C[1]);
+    x2matgen(n, n, D[0], D[1]);
+    x2matgen(n, n, E[0], E[1]);
+    x2matgen(n, n, F[0], F[1]);
+
+    for (i = 0; i < n; i++) {
+        // set diagonal
+        A[0][x1 * (i + n * i)] =
+        A[1][x1 * (i + n * i)] = (datatype) rand() / RAND_MAX;
+        E[0][x1 * (i + n * i)] =
+        E[1][x1 * (i + n * i)] = (datatype) rand() / RAND_MAX;
+        // clear first subdiagonal
+        A[0][x1 * (i + 1 + n * i)] =
+        A[1][x1 * (i + 1 + n * i)] =
+        B[0][x1 * (i + 1 + n * i)] =
+        B[1][x1 * (i + 1 + n * i)] =
+        A[0][x1 * (i + 1 + n * i) + x1 - 1] =
+        A[1][x1 * (i + 1 + n * i) + x1 - 1] =
+        B[0][x1 * (i + 1 + n * i) + x1 - 1] =
+        B[1][x1 * (i + 1 + n * i) + x1 - 1] = 0;
+    }
+}
+
+
+void post() {
+    if (scale[0] != 1 || scale[0] != 1)
+        printf("scale[RELAPACK] = %12g\tscale[LAPACK] = %12g\n", scale[0], scale[1]);
+    if (scale[0]) {
+        xlascl("G", iZERO, iZERO, &scale[0], &scale[1], &n, &n, C[0], &n, &info);
+        xlascl("G", iZERO, iZERO, &scale[0], &scale[1], &n, &n, F[0], &n, &info);
+    }
+    error = x2vecerr(n * n, C[0], C[1]) + x2vecerr(n * n, F[0], F[1]);
+}
+
+void tests() {
+    lWork = 2 * n * n;
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+    B[0] = xmalloc(n * n);
+    B[1] = xmalloc(n * n);
+    C[0] = xmalloc(n * n);
+    C[1] = xmalloc(n * n);
+    D[0] = xmalloc(n * n);
+    D[1] = xmalloc(n * n);
+    E[0] = xmalloc(n * n);
+    E[1] = xmalloc(n * n);
+    F[0] = xmalloc(n * n);
+    F[1] = xmalloc(n * n);
+    Work = xmalloc(lWork);
+    iWork = imalloc(n + n + 2);
+
+    #define ROUTINE XPREF(tgsyl)
+
+    TEST("N", iZERO, &n, &n, A[i], &n, B[i], &n, C[i], &n, D[i], &n, E[i], &n, F[i], &n, &scale[i], &dif[i], Work, &lWork, iWork, &info);
+    TEST("N", iZERO, &n2, &n, A[i], &n, B[i], &n, C[i], &n, D[i], &n, E[i], &n, F[i], &n, &scale[i], &dif[i], Work, &lWork, iWork, &info);
+    TEST("N", iZERO, &n, &n2, A[i], &n, B[i], &n, C[i], &n, D[i], &n, E[i], &n, F[i], &n, &scale[i], &dif[i], Work, &lWork, iWork, &info);
+    TEST("N", iONE, &n, &n, A[i], &n, B[i], &n, C[i], &n, D[i], &n, E[i], &n, F[i], &n, &scale[i], &dif[i], Work, &lWork, iWork, &info);
+    TEST("N", iTWO, &n, &n, A[i], &n, B[i], &n, C[i], &n, D[i], &n, E[i], &n, F[i], &n, &scale[i], &dif[i], Work, &lWork, iWork, &info);
+    TEST("N", iTHREE, &n, &n, A[i], &n, B[i], &n, C[i], &n, D[i], &n, E[i], &n, F[i], &n, &scale[i], &dif[i], Work, &lWork, iWork, &info);
+    TEST("N", iFOUR, &n, &n, A[i], &n, B[i], &n, C[i], &n, D[i], &n, E[i], &n, F[i], &n, &scale[i], &dif[i], Work, &lWork, iWork, &info);
+    TEST(xCTRANS, iZERO, &n, &n, A[i], &n, B[i], &n, C[i], &n, D[i], &n, E[i], &n, F[i], &n, &scale[i], &dif[i], Work, &lWork, iWork, &info);
+
+    free(A[0]);
+    free(A[1]);
+    free(B[0]);
+    free(B[1]);
+    free(C[0]);
+    free(C[1]);
+    free(D[0]);
+    free(D[1]);
+    free(E[0]);
+    free(E[1]);
+    free(F[0]);
+    free(F[1]);
+    free(Work);
+    free(iWork);
+}

relapack/test/xtrsyl.c

@@ -0,0 +1,65 @@
+#include "test.h"
+
+datatype *A[2], *B[2], *C[2], *Work, scale[2];
+int info;
+
+#define xlascl XPREF(LAPACK(lascl))
+void xlascl(const char *, const int *, const int *, const datatype *, const
+        datatype *, const int *, const int *, datatype *, const int *, int *);
+
+void pre() {
+    int i;
+
+    x2matgen(n, n, A[0], A[1]);
+    x2matgen(n, n, B[0], B[1]);
+    x2matgen(n, n, C[0], C[1]);
+
+    for (i = 0; i < n; i++) {
+        // set diagonal
+        A[0][x1 * (i + n * i)] =
+        A[1][x1 * (i + n * i)] = (datatype) rand() / RAND_MAX;
+        // clear first subdiagonal
+        A[0][x1 * (i + 1 + n * i)] =
+        A[1][x1 * (i + 1 + n * i)] =
+        B[0][x1 * (i + 1 + n * i)] =
+        B[1][x1 * (i + 1 + n * i)] =
+        A[0][x1 * (i + 1 + n * i) + x1 - 1] =
+        A[1][x1 * (i + 1 + n * i) + x1 - 1] =
+        B[0][x1 * (i + 1 + n * i) + x1 - 1] =
+        B[1][x1 * (i + 1 + n * i) + x1 - 1] = 0;
+    }
+}
+
+void post() {
+    if (scale[0] != 1 || scale[0] != 1)
+        printf("scale[RELAPACK] = %12g\tscale[LAPACK] = %12g\n", scale[0], scale[1]);
+    if (scale[0])
+        xlascl("G", iZERO, iZERO, &scale[0], &scale[1], &n, &n, C[0], &n, &info);
+    error = x2vecerr(n * n, C[0], C[1]);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+    B[0] = xmalloc(n * n);
+    B[1] = xmalloc(n * n);
+    C[0] = xmalloc(n * n);
+    C[1] = xmalloc(n * n);
+
+    #define ROUTINE XPREF(trsyl)
+
+    TEST("N", "N", iONE, &n, &n, A[i], &n, B[i], &n, C[i], &n, &scale[i], &info);
+    TEST("N", "N", iONE, &n2, &n, A[i], &n, B[i], &n, C[i], &n, &scale[i], &info);
+    TEST("N", "N", iONE, &n, &n2, A[i], &n, B[i], &n, C[i], &n, &scale[i], &info);
+    TEST("C", "N", iONE, &n, &n, A[i], &n, B[i], &n, C[i], &n, &scale[i], &info);
+    TEST("N", "C", iONE, &n, &n, A[i], &n, B[i], &n, C[i], &n, &scale[i], &info);
+    TEST("C", "C", iONE, &n, &n, A[i], &n, B[i], &n, C[i], &n, &scale[i], &info);
+    TEST("N", "N", iMONE, &n, &n, A[i], &n, B[i], &n, C[i], &n, &scale[i], &info);
+
+    free(A[0]);
+    free(A[1]);
+    free(B[0]);
+    free(B[1]);
+    free(C[0]);
+    free(C[1]);
+}

relapack/test/xtrtri.c

@@ -0,0 +1,25 @@
+#include "test.h"
+
+datatype *A[2];
+int info;
+
+void pre() {
+    x2matgen(n, n, A[0], A[1]);
+}
+
+void post() {
+    error = x2vecerr(n * n, A[0], A[1]);
+}
+
+void tests() {
+    A[0] = xmalloc(n * n);
+    A[1] = xmalloc(n * n);
+
+    #define ROUTINE XPREF(trtri)
+
+    TEST("L", "N", &n, A[i], &n, &info);
+    TEST("U", "N", &n, A[i], &n, &info);
+
+    free(A[0]);
+    free(A[1]);
+}