From b122413fb01aacff3112fdb387e74126539a8e5f Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Wed, 28 Jun 2017 18:13:14 +0200 Subject: [PATCH] Restore ReLAPACK test folder --- relapack/test/README.md | 48 +++++++++++++++++ relapack/test/config.h | 13 +++++ relapack/test/lapack.h | 64 +++++++++++++++++++++++ relapack/test/test.h | 136 ++++++++++++++++++++++++++++++++++++++++++++++++ relapack/test/util.c | 116 +++++++++++++++++++++++++++++++++++++++++ relapack/test/util.h | 15 ++++++ relapack/test/xgbtrf.c | 43 +++++++++++++++ relapack/test/xgemmt.c | 65 +++++++++++++++++++++++ relapack/test/xgetrf.c | 32 ++++++++++++ relapack/test/xhegst.c | 32 ++++++++++++ relapack/test/xhetrf.c | 40 ++++++++++++++ relapack/test/xlauum.c | 25 +++++++++ relapack/test/xpbtrf.c | 40 ++++++++++++++ relapack/test/xpotrf.c | 25 +++++++++ relapack/test/xsygst.c | 32 ++++++++++++ relapack/test/xsytrf.c | 40 ++++++++++++++ relapack/test/xtgsyl.c | 94 +++++++++++++++++++++++++++++++++ relapack/test/xtrsyl.c | 65 +++++++++++++++++++++++ relapack/test/xtrtri.c | 25 +++++++++ 19 files changed, 950 insertions(+) create mode 100644 relapack/test/README.md create mode 100644 relapack/test/config.h create mode 100644 relapack/test/lapack.h create mode 100644 relapack/test/test.h create mode 100644 relapack/test/util.c create mode 100644 relapack/test/util.h create mode 100644 relapack/test/xgbtrf.c create mode 100644 relapack/test/xgemmt.c create mode 100644 relapack/test/xgetrf.c create mode 100644 relapack/test/xhegst.c create mode 100644 relapack/test/xhetrf.c create mode 100644 relapack/test/xlauum.c create mode 100644 relapack/test/xpbtrf.c create mode 100644 relapack/test/xpotrf.c create mode 100644 relapack/test/xsygst.c create mode 100644 relapack/test/xsytrf.c create mode 100644 relapack/test/xtgsyl.c create mode 100644 relapack/test/xtrsyl.c create mode 100644 relapack/test/xtrtri.c diff --git a/relapack/test/README.md b/relapack/test/README.md new file mode 100644 index 0000000..48434b3 --- /dev/null +++ b/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. diff --git a/relapack/test/config.h b/relapack/test/config.h new file mode 100644 index 0000000..ab06a2f --- /dev/null +++ b/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 */ diff --git a/relapack/test/lapack.h b/relapack/test/lapack.h new file mode 100644 index 0000000..80f5c41 --- /dev/null +++ b/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 */ diff --git a/relapack/test/test.h b/relapack/test/test.h new file mode 100644 index 0000000..24089f3 --- /dev/null +++ b/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 +#include +#include + +// 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 */ diff --git a/relapack/test/util.c b/relapack/test/util.c new file mode 100644 index 0000000..e0fca3e --- /dev/null +++ b/relapack/test/util.c @@ -0,0 +1,116 @@ +#include "util.h" +#include +#include +#include + +#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; +} diff --git a/relapack/test/util.h b/relapack/test/util.h new file mode 100644 index 0000000..11d2999 --- /dev/null +++ b/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 */ diff --git a/relapack/test/xgbtrf.c b/relapack/test/xgbtrf.c new file mode 100644 index 0000000..f255006 --- /dev/null +++ b/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]); +} diff --git a/relapack/test/xgemmt.c b/relapack/test/xgemmt.c new file mode 100644 index 0000000..ffc3704 --- /dev/null +++ b/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); +} diff --git a/relapack/test/xgetrf.c b/relapack/test/xgetrf.c new file mode 100644 index 0000000..4484a24 --- /dev/null +++ b/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]); +} diff --git a/relapack/test/xhegst.c b/relapack/test/xhegst.c new file mode 100644 index 0000000..c318ef5 --- /dev/null +++ b/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]); +} diff --git a/relapack/test/xhetrf.c b/relapack/test/xhetrf.c new file mode 100644 index 0000000..b5d54bd --- /dev/null +++ b/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); +} diff --git a/relapack/test/xlauum.c b/relapack/test/xlauum.c new file mode 100644 index 0000000..d2c42fa --- /dev/null +++ b/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]); +} diff --git a/relapack/test/xpbtrf.c b/relapack/test/xpbtrf.c new file mode 100644 index 0000000..9a9babb --- /dev/null +++ b/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]); +} diff --git a/relapack/test/xpotrf.c b/relapack/test/xpotrf.c new file mode 100644 index 0000000..5e04d42 --- /dev/null +++ b/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]); +} diff --git a/relapack/test/xsygst.c b/relapack/test/xsygst.c new file mode 100644 index 0000000..b473a59 --- /dev/null +++ b/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]); +} diff --git a/relapack/test/xsytrf.c b/relapack/test/xsytrf.c new file mode 100644 index 0000000..82d626f --- /dev/null +++ b/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); +} diff --git a/relapack/test/xtgsyl.c b/relapack/test/xtgsyl.c new file mode 100644 index 0000000..74db500 --- /dev/null +++ b/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); +} diff --git a/relapack/test/xtrsyl.c b/relapack/test/xtrsyl.c new file mode 100644 index 0000000..358a892 --- /dev/null +++ b/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]); +} diff --git a/relapack/test/xtrtri.c b/relapack/test/xtrtri.c new file mode 100644 index 0000000..106391b --- /dev/null +++ b/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]); +} -- 2.7.4