diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0ef653d8cb..bcea2f625c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -24,7 +24,7 @@ set(CMAKE_CXX_EXTENSIONS OFF)
find_package(Boost)
find_package(MPI)
find_package(GDAL CONFIG)
-find_package(GSL)
+find_package(Eigen3)
find_package(BZip2)
# Build options - which programming models to enable
diff --git a/README.md b/README.md
index 2d5890abf7..85472bc473 100644
--- a/README.md
+++ b/README.md
@@ -10,11 +10,11 @@ Z. Jin and J. S. Vetter, "A Benchmark Suite for Improving Performance Portabilit
[NVIDIA HPC SDK](https://developer.nvidia.com/hpc-sdk)
# Dependencies
-Certain benchmarks require [Boost](https://www.boost.org/releases/latest/), [GSL](https://www.gnu.org/software/gsl), [GDAL](https://github.com/OSGeo/gdal), GPU-aware Message Passing Interface(MPI) or vendors' collective communication libraries (e.g. NCCL).
+Certain benchmarks require [Boost](https://www.boost.org/releases/latest/), [Eigen](https://eigen.tuxfamily.org), [GDAL](https://github.com/OSGeo/gdal), GPU-aware Message Passing Interface(MPI) or vendors' collective communication libraries (e.g. NCCL).
Boost: hbc, ge-spmm, mmcsf, warpsort, gerbil
MPI: miniDGS, miniWeather, pingpong, sparkler, allreduce, ccl, halo-finder
CCL: ccl
-GSL: sss, xlqc
+Eigen: xlqc
GDAL: stsg
BZip2: gerbil
diff --git a/cmake/modules/BenchmarkMacros.cmake b/cmake/modules/BenchmarkMacros.cmake
index 4f6d18bb2d..63528e0b6e 100644
--- a/cmake/modules/BenchmarkMacros.cmake
+++ b/cmake/modules/BenchmarkMacros.cmake
@@ -11,8 +11,8 @@ set(DEPEND_ON_BOOST "hbc" "ge-spmm" "mmcsf" "warpsort" "gerbil")
# Global list for benchmarks that require MPI
set(DEPEND_ON_MPI "miniDGS" "miniWeather" "pingpong" "sparkler" "allreduce" "ccl" "halo-finder")
-# Global list for benchmarks that require GSL
-set(DEPEND_ON_GSL "sss" "xlqc")
+# Global list for benchmarks that require Eigen
+set(DEPEND_ON_EIGEN "xlqc")
# Global list for benchmarks that require GDAL
set(DEPEND_ON_GDAL "stsg")
@@ -116,9 +116,10 @@ function(add_hecbench_benchmark)
endif()
endif()
- if(${BENCH_NAME} IN_LIST DEPEND_ON_GSL)
- if(NOT GSL_FOUND)
- message(STATUS "Skipping ${BENCH_NAME}-${BENCH_MODEL_LOWER} (GSL not found)")
+ if(${BENCH_NAME} IN_LIST DEPEND_ON_EIGEN)
+ # Eigen is header-only; the imported target carries the include path
+ if(NOT TARGET Eigen3::Eigen)
+ message(STATUS "Skipping ${BENCH_NAME}-${BENCH_MODEL_LOWER} (Eigen not found)")
return()
endif()
endif()
diff --git a/src/sss-cuda/CMakeLists.txt b/src/sss-cuda/CMakeLists.txt
index bbb0a0a282..1557a21b3e 100644
--- a/src/sss-cuda/CMakeLists.txt
+++ b/src/sss-cuda/CMakeLists.txt
@@ -4,7 +4,5 @@ add_hecbench_benchmark(
NAME sss
MODEL cuda
SOURCES main.cu
- INCLUDE_DIRS ${GSL_INCLUDE_DIRS}
CATEGORIES algorithms
- LINK_LIBRARIES gsl
)
diff --git a/src/sss-cuda/DPmixGGM.cpp b/src/sss-cuda/DPmixGGM.cpp
index 030ef06036..8dccb876e7 100644
--- a/src/sss-cuda/DPmixGGM.cpp
+++ b/src/sss-cuda/DPmixGGM.cpp
@@ -1,3 +1,4 @@
+#include
#define DPMIXGGM_CPP
#ifndef GRAPH_CPP
#include "graph.cpp"
diff --git a/src/sss-cuda/Makefile b/src/sss-cuda/Makefile
index 405688c7c9..4e298695ce 100755
--- a/src/sss-cuda/Makefile
+++ b/src/sss-cuda/Makefile
@@ -7,7 +7,6 @@ CC = nvcc
OPTIMIZE = yes
DEBUG = no
ARCH = sm_60
-GSL = /path/to/GSL
LAUNCHER =
#===============================================================================
@@ -25,11 +24,10 @@ obj = $(source:.cu=.o)
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Xcompiler -Wall -arch=$(ARCH) \
- -I$(GSL)/include
+CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Xcompiler -Wall -arch=$(ARCH)
# Linker Flags
-LDFLAGS = -L$(GSL)/lib -lgsl
+LDFLAGS =
# Debug Flags
ifeq ($(DEBUG),yes)
@@ -49,11 +47,11 @@ $(program): $(obj) Makefile
$(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS)
%.o: %.cu DPmixGGM.cpp DPmixGGM_Lists.cpp DPmixGGM_SSSmoves.cpp \
- kernels.cu graph.cpp gwish.cpp utilities.cpp Makefile
+ kernels.cu graph.cpp gwish.cpp utilities.cpp gsl_compat.h Makefile
$(CC) $(CFLAGS) -c $< -o $@
clean:
rm -rf $(program) $(obj) RES/f9_n150_p50_modes_GPU.txt
run: $(program)
- $(LAUNCHER) GSL_RNG_SEED=123 ./$(program) f9_n150_p50
+ $(LAUNCHER) ./$(program) f9_n150_p50
diff --git a/src/sss-cuda/README.md b/src/sss-cuda/README.md
index d020e7e2f9..d313bd8d90 100644
--- a/src/sss-cuda/README.md
+++ b/src/sss-cuda/README.md
@@ -1,22 +1,28 @@
DPmixGGM
========
-This folder contains source codes for the "GPU-powered Stochastic Shotgun Search for Dirichlet proces mixtures of Gaussian Graphical Models"
+This folder contains source codes for the "GPU-powered Stochastic Shotgun Search for Dirichlet proces mixtures of Gaussian Graphical Models"
by Chiranjit Mukherjee and Abel Rodriguez.
-
-The "DPmixGGM_SSS_main.cpp" file contains tuning parameters for the algorithm, as elaborated below:
+
+The "DPmixGGM_SSS_main.cpp" file contains tuning parameters for the algorithm, as elaborated below:
1. Run the SSS
-2. Run GPU/CPU versions of the SSS by enabling / disabling the macro CUDA.
-3. Specify maximum number of mixture components that the model should accommodate (for pre-allocation of memory).
-4. Set SSS runtime parameters C, D, R, S, M, g, h, f, t.
-5. Set SSS number of chain parameters. User needs to provide at least one initial point.
-7. Set hyperparameters of for the prior on (mu, K | G) with N0, DELTA0.
-
+2. Run GPU/CPU versions of the SSS by enabling / disabling the macro CUDA.
+3. Specify maximum number of mixture components that the model should accommodate (for pre-allocation of memory).
+4. Set SSS runtime parameters C, D, R, S, M, g, h, f, t.
+5. Set SSS number of chain parameters. User needs to provide at least one initial point.
+7. Set hyperparameters of for the prior on (mu, K | G) with N0, DELTA0.
+
Complie source codes using the "make" command and run with "./main f9_n150_p50" command.
-
-The program expects an input-data file (e.g. f9_n150_p50) in the DATA/ folder and at least one initialization point (e.g. f9_n150_p50_init1).
-The input-data file should specify n and p in the first row and then provide n rows of length p. The initial point data-file should specify n, p
-and L of the initial model configuration in the first row and xi-indices of the initial point in the second row. Subsequent L rows specify
-G_l (l=1:L).
-
+
+The program expects an input-data file (e.g. f9_n150_p50) in the DATA/ folder and at least one initialization point (e.g. f9_n150_p50_init1).
+The input-data file should specify n and p in the first row and then provide n rows of length p. The initial point data-file should specify n, p
+and L of the initial model configuration in the first row and xi-indices of the initial point in the second row. Subsequent L rows specify
+G_l (l=1:L).
+
A list of highest-score models is stored in folder RES/.
+
+Dependencies
+------------
+None beyond a C++17 compiler. The uniform random number generator (MT19937) and
+the adaptive Gauss-Kronrod quadrature the sampler needs are implemented in
+`gsl_compat.h` in this folder.
diff --git a/src/sss-cuda/gsl_compat.h b/src/sss-cuda/gsl_compat.h
new file mode 100644
index 0000000000..4369bd91ae
--- /dev/null
+++ b/src/sss-cuda/gsl_compat.h
@@ -0,0 +1,872 @@
+/*****************************************************************************
+ Self-contained replacement for the subset of the GNU Scientific Library used
+ by this benchmark, added so that the benchmark no longer links against the
+ GPL-licensed GSL (see HeCBench issue #319).
+
+ Two facilities are provided, both under the BSD 3-Clause license of this
+ benchmark:
+
+ 1. The MT19937 generator of Matsumoto and Nishimura, seeded and tempered
+ exactly as GSL's gsl_rng_mt19937 does, so that the random stream (and
+ therefore the stochastic search performed by this benchmark) is
+ unchanged.
+
+ 2. Adaptive Gauss-Kronrod quadrature over a semi-infinite interval,
+ following the QAGS/QAGI algorithm of QUADPACK (Piessens, de
+ Doncker-Kapenga, Ueberhuber and Kahaner; public domain), which is the
+ algorithm GSL's gsl_integration_qagiu implements.
+ *****************************************************************************/
+
+#ifndef GSL_COMPAT_H
+#define GSL_COMPAT_H
+
+#include
+#include
+#include
+
+#define GSL_DBL_EPSILON 2.2204460492503131e-16
+#define GSL_DBL_MIN 2.2250738585072014e-308
+#define GSL_DBL_MAX 1.7976931348623157e+308
+
+#define GSL_MAX_DBL(a, b) ((a) > (b) ? (a) : (b))
+
+#define GSL_SUCCESS 0
+#define GSL_EINVAL 4
+#define GSL_EMAXITER 11
+#define GSL_EBADTOL 13
+#define GSL_EROUND 18
+#define GSL_ESING 21
+#define GSL_EDIVERGE 22
+
+//--------------------------------------------------------------------
+// Random number generation: MT19937
+//--------------------------------------------------------------------
+
+typedef struct {
+ const char *name;
+ unsigned long int max;
+ unsigned long int min;
+} gsl_rng_type;
+
+#define GSL_COMPAT_MT_N 624
+#define GSL_COMPAT_MT_M 397
+
+typedef struct {
+ const gsl_rng_type *type;
+ unsigned long int mt[GSL_COMPAT_MT_N];
+ int mti;
+} gsl_rng;
+
+inline const gsl_rng_type gsl_rng_mt19937_type = {"mt19937", 0xffffffffUL, 0UL};
+inline const gsl_rng_type *gsl_rng_mt19937 = &gsl_rng_mt19937_type;
+inline const gsl_rng_type *gsl_rng_default = &gsl_rng_mt19937_type;
+inline unsigned long int gsl_rng_default_seed = 0;
+
+inline void gsl_rng_set(gsl_rng *r, unsigned long int s) {
+ // 4357 is the seed GSL substitutes for zero
+ if (s == 0) s = 4357;
+
+ r->mt[0] = s & 0xffffffffUL;
+
+ int i;
+ for (i = 1; i < GSL_COMPAT_MT_N; i++) {
+ r->mt[i] = (1812433253UL * (r->mt[i - 1] ^ (r->mt[i - 1] >> 30)) +
+ (unsigned long int)i) & 0xffffffffUL;
+ }
+ r->mti = i;
+}
+
+inline unsigned long int gsl_rng_get(gsl_rng *r) {
+ const unsigned long int UPPER_MASK = 0x80000000UL;
+ const unsigned long int LOWER_MASK = 0x7fffffffUL;
+ unsigned long int *const mt = r->mt;
+
+ if (r->mti >= GSL_COMPAT_MT_N) {
+ int kk;
+ unsigned long int y;
+
+ for (kk = 0; kk < GSL_COMPAT_MT_N - GSL_COMPAT_MT_M; kk++) {
+ y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
+ mt[kk] = mt[kk + GSL_COMPAT_MT_M] ^ (y >> 1) ^
+ ((y & 0x1UL) ? 0x9908b0dfUL : 0UL);
+ }
+ for (; kk < GSL_COMPAT_MT_N - 1; kk++) {
+ y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
+ mt[kk] = mt[kk + (GSL_COMPAT_MT_M - GSL_COMPAT_MT_N)] ^ (y >> 1) ^
+ ((y & 0x1UL) ? 0x9908b0dfUL : 0UL);
+ }
+ y = (mt[GSL_COMPAT_MT_N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
+ mt[GSL_COMPAT_MT_N - 1] = mt[GSL_COMPAT_MT_M - 1] ^ (y >> 1) ^
+ ((y & 0x1UL) ? 0x9908b0dfUL : 0UL);
+
+ r->mti = 0;
+ }
+
+ // Tempering. The state is held masked to 32 bits, and both shift-and-mask
+ // steps use 32-bit masks, so k stays within 32 bits throughout.
+ unsigned long int k = mt[r->mti];
+ k ^= (k >> 11);
+ k ^= (k << 7) & 0x9d2c5680UL;
+ k ^= (k << 15) & 0xefc60000UL;
+ k ^= (k >> 18);
+
+ r->mti++;
+
+ return k;
+}
+
+inline double gsl_rng_uniform(gsl_rng *r) {
+ return gsl_rng_get(r) / 4294967296.0;
+}
+
+inline gsl_rng *gsl_rng_alloc(const gsl_rng_type *T) {
+ gsl_rng *r = (gsl_rng *)malloc(sizeof(gsl_rng));
+ if (r == NULL) return NULL;
+ r->type = T;
+ gsl_rng_set(r, gsl_rng_default_seed);
+ return r;
+}
+
+inline void gsl_rng_free(gsl_rng *r) { free(r); }
+
+// Honours GSL_RNG_SEED so that the command lines shipped with the benchmark
+// keep working. GSL_RNG_TYPE is ignored: only the default generator exists.
+inline const gsl_rng_type *gsl_rng_env_setup(void) {
+ const char *seed = getenv("GSL_RNG_SEED");
+ if (seed != NULL) {
+ gsl_rng_default_seed = strtoul(seed, NULL, 0);
+ }
+ return gsl_rng_default;
+}
+
+inline double gsl_ran_flat(gsl_rng *r, const double a, const double b) {
+ double u = gsl_rng_uniform(r);
+ return a * (1 - u) + b * u;
+}
+
+//--------------------------------------------------------------------
+// Numerical integration
+//--------------------------------------------------------------------
+
+typedef struct {
+ double (*function)(double x, void *params);
+ void *params;
+} gsl_function;
+
+#define GSL_FN_EVAL(F, x) (*((F)->function))((x), (F)->params)
+
+typedef struct {
+ size_t limit;
+ size_t size;
+ size_t nrmax;
+ size_t i;
+ size_t maximum_level;
+ double *alist;
+ double *blist;
+ double *rlist;
+ double *elist;
+ size_t *order;
+ size_t *level;
+} gsl_integration_workspace;
+
+inline gsl_integration_workspace *
+gsl_integration_workspace_alloc(const size_t n) {
+ if (n == 0) return NULL;
+
+ gsl_integration_workspace *w =
+ (gsl_integration_workspace *)malloc(sizeof(gsl_integration_workspace));
+ if (w == NULL) return NULL;
+
+ w->alist = (double *)malloc(n * sizeof(double));
+ w->blist = (double *)malloc(n * sizeof(double));
+ w->rlist = (double *)malloc(n * sizeof(double));
+ w->elist = (double *)malloc(n * sizeof(double));
+ w->order = (size_t *)malloc(n * sizeof(size_t));
+ w->level = (size_t *)malloc(n * sizeof(size_t));
+
+ w->limit = n;
+ w->size = 0;
+ w->nrmax = 0;
+ w->i = 0;
+ w->maximum_level = 0;
+
+ return w;
+}
+
+inline void gsl_integration_workspace_free(gsl_integration_workspace *w) {
+ if (w == NULL) return;
+ free(w->level);
+ free(w->order);
+ free(w->elist);
+ free(w->rlist);
+ free(w->blist);
+ free(w->alist);
+ free(w);
+}
+
+//---------------- 15-point Gauss-Kronrod rule -----------------------
+
+inline double gsl_compat_rescale_error(double err, const double result_abs,
+ const double result_asc) {
+ err = fabs(err);
+
+ if (result_asc != 0 && err != 0) {
+ double scale = pow((200 * err / result_asc), 1.5);
+ if (scale < 1)
+ err = result_asc * scale;
+ else
+ err = result_asc;
+ }
+
+ if (result_abs > GSL_DBL_MIN / (50 * GSL_DBL_EPSILON)) {
+ double min_err = 50 * GSL_DBL_EPSILON * result_abs;
+ if (min_err > err) err = min_err;
+ }
+
+ return err;
+}
+
+inline void gsl_integration_qk15(const gsl_function *f, double a, double b,
+ double *result, double *abserr, double *resabs,
+ double *resasc) {
+ // Abscissae of the 15-point Kronrod rule
+ static const double xgk[8] = {
+ 0.991455371120812639206854697526329,
+ 0.949107912342758524526189684047851,
+ 0.864864423359769072789712788640926,
+ 0.741531185599394439863864773280788,
+ 0.586087235467691130294144838258730,
+ 0.405845151377397166906606412076961,
+ 0.207784955007898467600689403773245,
+ 0.000000000000000000000000000000000};
+
+ // Weights of the 15-point Kronrod rule
+ static const double wgk[8] = {
+ 0.022935322010529224963732008058970,
+ 0.063092092629978553290700663189204,
+ 0.104790010322250183839876322541518,
+ 0.140653259715525918745189590510238,
+ 0.169004726639267902826583426598550,
+ 0.190350578064785409913256402421014,
+ 0.204432940075298892414161999234649,
+ 0.209482141084727828012999174891714};
+
+ // Weights of the embedded 7-point Gauss rule
+ static const double wg[4] = {0.129484966168869693270611432679082,
+ 0.279705391489276667901467771423780,
+ 0.381830050505118944950369775488975,
+ 0.417959183673469387755102040816327};
+
+ const int n = 8;
+ double fv1[8], fv2[8];
+
+ const double center = 0.5 * (a + b);
+ const double half_length = 0.5 * (b - a);
+ const double abs_half_length = fabs(half_length);
+ const double f_center = GSL_FN_EVAL(f, center);
+
+ double result_gauss = f_center * wg[n / 2 - 1];
+ double result_kronrod = f_center * wgk[n - 1];
+ double result_abs = fabs(result_kronrod);
+ double result_asc = 0;
+
+ int j;
+
+ for (j = 0; j < (n - 1) / 2; j++) {
+ const int jtw = j * 2 + 1;
+ const double abscissa = half_length * xgk[jtw];
+ const double fval1 = GSL_FN_EVAL(f, center - abscissa);
+ const double fval2 = GSL_FN_EVAL(f, center + abscissa);
+ const double fsum = fval1 + fval2;
+ fv1[jtw] = fval1;
+ fv2[jtw] = fval2;
+ result_gauss += wg[j] * fsum;
+ result_kronrod += wgk[jtw] * fsum;
+ result_abs += wgk[jtw] * (fabs(fval1) + fabs(fval2));
+ }
+
+ for (j = 0; j < n / 2; j++) {
+ const int jtwm1 = j * 2;
+ const double abscissa = half_length * xgk[jtwm1];
+ const double fval1 = GSL_FN_EVAL(f, center - abscissa);
+ const double fval2 = GSL_FN_EVAL(f, center + abscissa);
+ fv1[jtwm1] = fval1;
+ fv2[jtwm1] = fval2;
+ result_kronrod += wgk[jtwm1] * (fval1 + fval2);
+ result_abs += wgk[jtwm1] * (fabs(fval1) + fabs(fval2));
+ }
+
+ const double mean = result_kronrod * 0.5;
+
+ result_asc = wgk[n - 1] * fabs(f_center - mean);
+ for (j = 0; j < n - 1; j++) {
+ result_asc += wgk[j] * (fabs(fv1[j] - mean) + fabs(fv2[j] - mean));
+ }
+
+ const double err = (result_kronrod - result_gauss) * half_length;
+
+ result_kronrod *= half_length;
+ result_abs *= abs_half_length;
+ result_asc *= abs_half_length;
+
+ *result = result_kronrod;
+ *resabs = result_abs;
+ *resasc = result_asc;
+ *abserr = gsl_compat_rescale_error(err, result_abs, result_asc);
+}
+
+//---------------- Subinterval bookkeeping ---------------------------
+
+inline void gsl_compat_initialise(gsl_integration_workspace *w, double a,
+ double b) {
+ w->size = 0;
+ w->nrmax = 0;
+ w->i = 0;
+ w->alist[0] = a;
+ w->blist[0] = b;
+ w->rlist[0] = 0.0;
+ w->elist[0] = 0.0;
+ w->order[0] = 0;
+ w->level[0] = 0;
+ w->maximum_level = 0;
+}
+
+inline void gsl_compat_set_initial_result(gsl_integration_workspace *w,
+ double result, double error) {
+ w->size = 1;
+ w->rlist[0] = result;
+ w->elist[0] = error;
+}
+
+// Maintains the list of subintervals sorted by descending error estimate.
+inline void gsl_compat_qpsrt(gsl_integration_workspace *w) {
+ const size_t last = w->size - 1;
+ const size_t limit = w->limit;
+ double *elist = w->elist;
+ size_t *order = w->order;
+
+ size_t i_nrmax = w->nrmax;
+ size_t i_maxerr = order[i_nrmax];
+
+ if (last < 2) {
+ order[0] = 0;
+ order[1] = 1;
+ w->i = i_maxerr;
+ return;
+ }
+
+ const double errmax = elist[i_maxerr];
+
+ while (i_nrmax > 0 && errmax > elist[order[i_nrmax - 1]]) {
+ order[i_nrmax] = order[i_nrmax - 1];
+ i_nrmax--;
+ }
+
+ int top;
+ if (last < (limit / 2 + 2))
+ top = (int)last;
+ else
+ top = (int)(limit - last + 1);
+
+ int i = (int)i_nrmax + 1;
+
+ while (i < top && errmax < elist[order[i]]) {
+ order[i - 1] = order[i];
+ i++;
+ }
+
+ order[i - 1] = i_maxerr;
+
+ const double errmin = elist[last];
+
+ int k = top - 1;
+
+ while (k > i - 2 && errmin >= elist[order[k]]) {
+ order[k + 1] = order[k];
+ k--;
+ }
+
+ order[k + 1] = last;
+
+ i_maxerr = order[i_nrmax];
+
+ w->i = i_maxerr;
+ w->nrmax = i_nrmax;
+}
+
+inline void gsl_compat_retrieve(const gsl_integration_workspace *w, double *a,
+ double *b, double *r, double *e) {
+ const size_t i = w->i;
+ *a = w->alist[i];
+ *b = w->blist[i];
+ *r = w->rlist[i];
+ *e = w->elist[i];
+}
+
+inline void gsl_compat_update(gsl_integration_workspace *w, double a1, double b1,
+ double area1, double error1, double a2, double b2,
+ double area2, double error2) {
+ const size_t i_max = w->i;
+ const size_t i_new = w->size;
+ const size_t new_level = w->level[i_max] + 1;
+
+ if (error2 > error1) {
+ w->alist[i_max] = a2; // blist[i_max] already holds b2
+ w->rlist[i_max] = area2;
+ w->elist[i_max] = error2;
+ w->level[i_max] = new_level;
+
+ w->alist[i_new] = a1;
+ w->blist[i_new] = b1;
+ w->rlist[i_new] = area1;
+ w->elist[i_new] = error1;
+ w->level[i_new] = new_level;
+ } else {
+ w->blist[i_max] = b1; // alist[i_max] already holds a1
+ w->rlist[i_max] = area1;
+ w->elist[i_max] = error1;
+ w->level[i_max] = new_level;
+
+ w->alist[i_new] = a2;
+ w->blist[i_new] = b2;
+ w->rlist[i_new] = area2;
+ w->elist[i_new] = error2;
+ w->level[i_new] = new_level;
+ }
+
+ w->size++;
+
+ if (new_level > w->maximum_level) w->maximum_level = new_level;
+
+ gsl_compat_qpsrt(w);
+}
+
+inline double gsl_compat_sum_results(const gsl_integration_workspace *w) {
+ double result_sum = 0;
+ for (size_t k = 0; k < w->size; k++) result_sum += w->rlist[k];
+ return result_sum;
+}
+
+inline int gsl_compat_subinterval_too_small(double a1, double a2, double b2) {
+ const double tmp =
+ (1 + 100 * GSL_DBL_EPSILON) * (fabs(a2) + 1000 * GSL_DBL_MIN);
+ return fabs(a1) <= tmp && fabs(b2) <= tmp;
+}
+
+inline void gsl_compat_reset_nrmax(gsl_integration_workspace *w) {
+ w->nrmax = 0;
+ w->i = w->order[0];
+}
+
+inline int gsl_compat_increase_nrmax(gsl_integration_workspace *w) {
+ const size_t id = w->nrmax;
+ const size_t last = w->size - 1;
+
+ size_t jupbnd;
+ if (last > (1 + w->limit / 2))
+ jupbnd = w->limit + 1 - last;
+ else
+ jupbnd = last;
+
+ for (size_t k = id; k <= jupbnd; k++) {
+ const size_t i_max = w->order[w->nrmax];
+ w->i = i_max;
+ if (w->level[i_max] < w->maximum_level) return 1;
+ w->nrmax++;
+ }
+ return 0;
+}
+
+inline int gsl_compat_large_interval(gsl_integration_workspace *w) {
+ return w->level[w->i] < w->maximum_level;
+}
+
+//---------------- Wynn epsilon extrapolation ------------------------
+
+typedef struct {
+ size_t n;
+ double rlist2[52];
+ size_t nres;
+ double res3la[3];
+} gsl_compat_extrap_table;
+
+inline void gsl_compat_initialise_table(gsl_compat_extrap_table *table) {
+ table->n = 0;
+ table->nres = 0;
+}
+
+inline void gsl_compat_append_table(gsl_compat_extrap_table *table, double y) {
+ table->rlist2[table->n] = y;
+ table->n++;
+}
+
+inline void gsl_compat_qelg(gsl_compat_extrap_table *table, double *result,
+ double *abserr) {
+ double *epstab = table->rlist2;
+ double *res3la = table->res3la;
+ const size_t n = table->n - 1;
+
+ const double current = epstab[n];
+
+ double absolute = GSL_DBL_MAX;
+ double relative = 5 * GSL_DBL_EPSILON * fabs(current);
+
+ const size_t newelm = n / 2;
+ const size_t n_orig = n;
+ size_t n_final = n;
+ const size_t nres_orig = table->nres;
+
+ *result = current;
+ *abserr = GSL_DBL_MAX;
+
+ if (n < 2) {
+ *abserr = GSL_MAX_DBL(absolute, relative);
+ return;
+ }
+
+ epstab[n + 2] = epstab[n];
+ epstab[n] = GSL_DBL_MAX;
+
+ for (size_t i = 0; i < newelm; i++) {
+ double res = epstab[n - 2 * i + 2];
+ const double e0 = epstab[n - 2 * i - 2];
+ const double e1 = epstab[n - 2 * i - 1];
+ const double e2 = res;
+
+ const double e1abs = fabs(e1);
+ const double delta2 = e2 - e1;
+ const double err2 = fabs(delta2);
+ const double tol2 = GSL_MAX_DBL(fabs(e2), e1abs) * GSL_DBL_EPSILON;
+ const double delta3 = e1 - e0;
+ const double err3 = fabs(delta3);
+ const double tol3 = GSL_MAX_DBL(e1abs, fabs(e0)) * GSL_DBL_EPSILON;
+
+ if (err2 <= tol2 && err3 <= tol3) {
+ // e0, e1 and e2 agree to machine accuracy: assume convergence
+ *result = res;
+ absolute = err2 + err3;
+ relative = 5 * GSL_DBL_EPSILON * fabs(res);
+ *abserr = GSL_MAX_DBL(absolute, relative);
+ return;
+ }
+
+ const double e3 = epstab[n - 2 * i];
+ epstab[n - 2 * i] = e1;
+ const double delta1 = e1 - e3;
+ const double err1 = fabs(delta1);
+ const double tol1 = GSL_MAX_DBL(e1abs, fabs(e3)) * GSL_DBL_EPSILON;
+
+ // Drop part of the table when two elements are nearly equal, or when the
+ // table behaves irregularly
+ if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3) {
+ n_final = 2 * i;
+ break;
+ }
+
+ const double ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
+
+ if (fabs(ss * e1) <= 0.0001) {
+ n_final = 2 * i;
+ break;
+ }
+
+ res = e1 + 1 / ss;
+ epstab[n - 2 * i] = res;
+
+ const double error = err2 + fabs(res - e2) + err3;
+ if (error <= *abserr) {
+ *abserr = error;
+ *result = res;
+ }
+ }
+
+ const size_t limexp = 50 - 1;
+ if (n_final == limexp) n_final = 2 * (limexp / 2);
+
+ if (n_orig % 2 == 1) {
+ for (size_t i = 0; i <= newelm; i++) epstab[1 + i * 2] = epstab[i * 2 + 3];
+ } else {
+ for (size_t i = 0; i <= newelm; i++) epstab[i * 2] = epstab[i * 2 + 2];
+ }
+
+ if (n_orig != n_final) {
+ for (size_t i = 0; i <= n_final; i++)
+ epstab[i] = epstab[n_orig - n_final + i];
+ }
+
+ table->n = n_final + 1;
+
+ if (nres_orig < 3) {
+ res3la[nres_orig] = *result;
+ *abserr = GSL_DBL_MAX;
+ } else {
+ *abserr = (fabs(*result - res3la[2]) + fabs(*result - res3la[1]) +
+ fabs(*result - res3la[0]));
+ res3la[0] = res3la[1];
+ res3la[1] = res3la[2];
+ res3la[2] = *result;
+ }
+
+ table->nres = nres_orig + 1;
+
+ *abserr = GSL_MAX_DBL(*abserr, 5 * GSL_DBL_EPSILON * fabs(*result));
+}
+
+//---------------- Adaptive integration with extrapolation -----------
+
+inline int gsl_compat_positive(double result, double resabs) {
+ return fabs(result) >= (1 - 50 * GSL_DBL_EPSILON) * resabs;
+}
+
+inline int gsl_compat_qags(const gsl_function *f, const double a, const double b,
+ const double epsabs, const double epsrel,
+ const size_t limit, gsl_integration_workspace *w,
+ double *result, double *abserr) {
+ double result0, abserr0, resabs0, resasc0;
+ double tolerance;
+ double ertest = 0;
+ double error_over_large_intervals = 0;
+ double reseps = 0, abseps = 0, correc = 0;
+ size_t ktmin = 0;
+ int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
+ int error_type = 0, error_type2 = 0;
+ int extrapolate = 0, disallow_extrapolation = 0;
+
+ gsl_compat_extrap_table table;
+
+ gsl_compat_initialise(w, a, b);
+
+ *result = 0;
+ *abserr = 0;
+
+ if (limit > w->limit) return GSL_EINVAL;
+
+ if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
+ return GSL_EBADTOL;
+
+ gsl_integration_qk15(f, a, b, &result0, &abserr0, &resabs0, &resasc0);
+
+ gsl_compat_set_initial_result(w, result0, abserr0);
+
+ tolerance = GSL_MAX_DBL(epsabs, epsrel * fabs(result0));
+
+ if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance) {
+ *result = result0;
+ *abserr = abserr0;
+ return GSL_EROUND;
+ } else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0) {
+ *result = result0;
+ *abserr = abserr0;
+ return GSL_SUCCESS;
+ } else if (limit == 1) {
+ *result = result0;
+ *abserr = abserr0;
+ return GSL_EMAXITER;
+ }
+
+ gsl_compat_initialise_table(&table);
+ gsl_compat_append_table(&table, result0);
+
+ double area = result0;
+ double errsum = abserr0;
+ double res_ext = result0;
+ double err_ext = GSL_DBL_MAX;
+
+ const int positive_integrand = gsl_compat_positive(result0, resabs0);
+
+ size_t iteration = 1;
+ int compute_result = 0;
+
+ do {
+ double a_i, b_i, r_i, e_i;
+ double area1 = 0, area2 = 0;
+ double error1 = 0, error2 = 0;
+ double resabs1, resabs2, resasc1, resasc2;
+
+ // Bisect the subinterval carrying the largest error estimate
+ gsl_compat_retrieve(w, &a_i, &b_i, &r_i, &e_i);
+
+ const size_t current_level = w->level[w->i] + 1;
+
+ const double a1 = a_i;
+ const double b1 = 0.5 * (a_i + b_i);
+ const double a2 = b1;
+ const double b2 = b_i;
+
+ iteration++;
+
+ gsl_integration_qk15(f, a1, b1, &area1, &error1, &resabs1, &resasc1);
+ gsl_integration_qk15(f, a2, b2, &area2, &error2, &resabs2, &resasc2);
+
+ const double area12 = area1 + area2;
+ const double error12 = error1 + error2;
+ const double last_e_i = e_i;
+
+ errsum = errsum + error12 - e_i;
+ area = area + area12 - r_i;
+
+ tolerance = GSL_MAX_DBL(epsabs, epsrel * fabs(area));
+
+ if (resasc1 != error1 && resasc2 != error2) {
+ const double delta = r_i - area12;
+
+ if (fabs(delta) <= 1.0e-5 * fabs(area12) && error12 >= 0.99 * e_i) {
+ if (!extrapolate)
+ roundoff_type1++;
+ else
+ roundoff_type2++;
+ }
+ if (iteration > 10 && error12 > e_i) roundoff_type3++;
+ }
+
+ if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
+ error_type = 2;
+
+ if (roundoff_type2 >= 5) error_type2 = 1;
+
+ // Bad integrand behaviour at a point of the integration range
+ if (gsl_compat_subinterval_too_small(a1, a2, b2)) error_type = 4;
+
+ gsl_compat_update(w, a1, b1, area1, error1, a2, b2, area2, error2);
+
+ if (errsum <= tolerance) {
+ compute_result = 1;
+ break;
+ }
+
+ if (error_type) break;
+
+ if (iteration >= limit - 1) {
+ error_type = 1;
+ break;
+ }
+
+ if (iteration == 2) {
+ error_over_large_intervals = errsum;
+ ertest = tolerance;
+ gsl_compat_append_table(&table, area);
+ continue;
+ }
+
+ if (disallow_extrapolation) continue;
+
+ error_over_large_intervals += -last_e_i;
+
+ if (current_level < w->maximum_level)
+ error_over_large_intervals += error12;
+
+ if (!extrapolate) {
+ // Extrapolate only once the interval to be bisected next is the
+ // smallest one
+ if (gsl_compat_large_interval(w)) continue;
+
+ extrapolate = 1;
+ w->nrmax = 1;
+ }
+
+ if (!error_type2 && error_over_large_intervals > ertest) {
+ if (gsl_compat_increase_nrmax(w)) continue;
+ }
+
+ gsl_compat_append_table(&table, area);
+
+ gsl_compat_qelg(&table, &reseps, &abseps);
+
+ ktmin++;
+
+ if (ktmin > 5 && err_ext < 0.001 * errsum) error_type = 5;
+
+ if (abseps < err_ext) {
+ ktmin = 0;
+ err_ext = abseps;
+ res_ext = reseps;
+ correc = error_over_large_intervals;
+ ertest = GSL_MAX_DBL(epsabs, epsrel * fabs(reseps));
+ if (err_ext <= ertest) break;
+ }
+
+ if (table.n == 1) disallow_extrapolation = 1;
+
+ if (error_type == 5) break;
+
+ gsl_compat_reset_nrmax(w);
+ extrapolate = 0;
+ error_over_large_intervals = errsum;
+
+ } while (iteration < limit);
+
+ if (!compute_result) {
+ *result = res_ext;
+ *abserr = err_ext;
+
+ if (err_ext == GSL_DBL_MAX) {
+ compute_result = 1;
+ } else if (error_type || error_type2) {
+ if (error_type2) err_ext += correc;
+
+ if (error_type == 0) error_type = 3;
+
+ if (res_ext != 0.0 && area != 0.0) {
+ if (err_ext / fabs(res_ext) > errsum / fabs(area)) compute_result = 1;
+ } else if (err_ext > errsum) {
+ compute_result = 1;
+ }
+ }
+
+ if (!compute_result) {
+ // Test on divergence
+ const double max_area = GSL_MAX_DBL(fabs(res_ext), fabs(area));
+ if (!(!positive_integrand && max_area < 0.01 * resabs0)) {
+ const double ratio = res_ext / area;
+ if (ratio < 0.01 || ratio > 100.0 || errsum > fabs(area))
+ error_type = 6;
+ }
+ }
+ }
+
+ if (compute_result) {
+ *result = gsl_compat_sum_results(w);
+ *abserr = errsum;
+ }
+
+ if (error_type > 2) error_type--;
+
+ switch (error_type) {
+ case 0: return GSL_SUCCESS;
+ case 1: return GSL_EMAXITER;
+ case 2: return GSL_EROUND;
+ case 3: return GSL_ESING;
+ default: return GSL_EDIVERGE;
+ }
+}
+
+typedef struct {
+ double a;
+ const gsl_function *f;
+} gsl_compat_iu_params;
+
+// Maps (a, +infinity) onto (0, 1] through x = a + (1 - t) / t
+inline double gsl_compat_iu_transform(double t, void *params) {
+ gsl_compat_iu_params *p = (gsl_compat_iu_params *)params;
+ const double x = p->a + (1 - t) / t;
+ const double y = GSL_FN_EVAL(p->f, x);
+ return (y / t) / t;
+}
+
+inline int gsl_integration_qagiu(gsl_function *f, double a, double epsabs,
+ double epsrel, size_t limit,
+ gsl_integration_workspace *w, double *result,
+ double *abserr) {
+ gsl_compat_iu_params transform_params;
+ transform_params.a = a;
+ transform_params.f = f;
+
+ gsl_function f_transform;
+ f_transform.function = &gsl_compat_iu_transform;
+ f_transform.params = &transform_params;
+
+ return gsl_compat_qags(&f_transform, 0.0, 1.0, epsabs, epsrel, limit, w,
+ result, abserr);
+}
+
+#endif // GSL_COMPAT_H
diff --git a/src/sss-cuda/gwish.cpp b/src/sss-cuda/gwish.cpp
index 343adf69ef..03c44cbafe 100644
--- a/src/sss-cuda/gwish.cpp
+++ b/src/sss-cuda/gwish.cpp
@@ -1,3 +1,4 @@
+#include
// gwish.cpp: This is a collection of device functions that manipulate graph
// objects according to the G-Wishart distribution. Note: this depends on the
// graph.cpp library Chiranjit Mukherjee : chiranjit@soe.ucsc.edu -- based on
diff --git a/src/sss-cuda/kernels.cu b/src/sss-cuda/kernels.cu
index e1218b4da0..d5b33e29b7 100644
--- a/src/sss-cuda/kernels.cu
+++ b/src/sss-cuda/kernels.cu
@@ -23,12 +23,15 @@ __global__ void CanDeleteEdge(myInt *d_in_delete, myInt *isDecomposable) {
}
__shared__ myInt contain_a, contain_b, which_ab;
+ SYNC;
+
for (i = 0; i < nCliques; i++) {
ii = i * n;
if (tid == 0) {
contain_a = 0;
contain_b = 0;
}
+ SYNC;
for (j = tid; j < CliquesDimens[i]; j += bdim) {
k = Cliques[ii + j];
if (k == a) {
@@ -38,12 +41,16 @@ __global__ void CanDeleteEdge(myInt *d_in_delete, myInt *isDecomposable) {
contain_b = 1;
}
}
+ SYNC;
if (tid == 0) {
if (contain_a && contain_b) {
count++;
which_ab = i;
}
}
+ // count is updated by thread 0 only, so it must be visible to every thread
+ // before the loop-exit test below, or threads leave the loop divergently.
+ SYNC;
if (count > 1) {
break;
}
@@ -145,6 +152,7 @@ __global__ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add,
contain_b = 0;
};
contain_S[tid] = 0;
+ SYNC;
t = i * n;
for (j = tid; j < CliquesDimens[i]; j += bdim) {
c = Cliques[t + j];
@@ -160,6 +168,9 @@ __global__ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add,
}
}
}
+ // every thread contributes to contain_S[] and contain_a/contain_b above,
+ // so those writes must be visible before thread 0 reduces them.
+ SYNC;
if (tid == 0) {
k = 0;
for (j = 0; j < BLOCK_SIZE; j++) {
@@ -172,6 +183,7 @@ __global__ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add,
bSi = i;
}
}
+ SYNC;
}
if (tid == 0) { // find the path from aSi to root
@@ -208,6 +220,8 @@ __global__ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add,
}
}
}
+ // R/pR are produced by thread 0 and T/pT by thread 1; both are read below.
+ SYNC;
if (tid == 0) {
// find the branching point
diff --git a/src/sss-cuda/main.cu b/src/sss-cuda/main.cu
index c06d9325c2..00898130ff 100644
--- a/src/sss-cuda/main.cu
+++ b/src/sss-cuda/main.cu
@@ -66,13 +66,11 @@
#define ISFLOAT 10
using namespace std;
-#include
-#include
+#include "gsl_compat.h"
#define GSL_INTEGRATION_GRIDSIZE 1000
gsl_integration_workspace *w;
gsl_function F;
-#include
#define RANDOMSEED 314159265
// Define hyperparameters for the prior distribution of (mu, K | G)
diff --git a/src/sss-hip/CMakeLists.txt b/src/sss-hip/CMakeLists.txt
index 7bc9f0978a..e8884df559 100644
--- a/src/sss-hip/CMakeLists.txt
+++ b/src/sss-hip/CMakeLists.txt
@@ -4,7 +4,6 @@ add_hecbench_benchmark(
NAME sss
MODEL hip
SOURCES main.cu
- INCLUDE_DIRS ${GSL_INCLUDE_DIRS}
+ INCLUDE_DIRS ${CMAKE_CURRENT_LIST_DIR}/../sss-cuda
CATEGORIES algorithms
- LINK_LIBRARIES gsl
)
diff --git a/src/sss-hip/DPmixGGM.cpp b/src/sss-hip/DPmixGGM.cpp
deleted file mode 100644
index 030ef06036..0000000000
--- a/src/sss-hip/DPmixGGM.cpp
+++ /dev/null
@@ -1,515 +0,0 @@
-#define DPMIXGGM_CPP
-#ifndef GRAPH_CPP
-#include "graph.cpp"
-#endif
-#ifndef GWISH_CPP
-#include "gwish.cpp"
-#endif
-
-typedef class DPmixGGM *State;
-
-class DPmixGGM {
- //-------- Variables ---------------
-public:
- int n;
- int p;
- Real *X;
- myInt *xi;
- myInt L;
- LPGraph *graphlist;
- Real plp; // log-prior of the partition
- Real *pll; // array of log-likelihood for the cluster
- Real alpha;
-
- //--------- Functions --------------
-public:
- DPmixGGM(Real *data, myInt L_start, myInt n_obs, myInt p_model,
- Real edgeInclusionProb, ifstream &initfile);
- DPmixGGM(State a);
- ~DPmixGGM();
-
- void RandomStartAllXi(myInt L);
- void RandomStartAllG(myInt L, Real edgeInclusionProb);
-
- void RandomStart(myInt L, Real edgeInclusionProb);
- void InformedStart(ifstream &initfile);
-
- void ReadState(ifstream &initfile);
- void WriteState(ofstream &out, int itr);
- void CopyState(State a);
-
- Real partitionlogPrior(myInt theL, myInt *thexi, Real thealpha);
- Real lA(myInt the_n, myInt m, myInt k);
-
- Real cluster_k_loglikelihood(myInt k, myInt *thexi, LPGraph thegraph);
- Real predictiveDistribution(myInt k, myInt l, myInt *thexi, LPGraph thegraph);
-};
-
-//--------- This is the initializer ----------------------
-DPmixGGM::DPmixGGM(Real *data, myInt L_start, myInt n_obs, myInt p_model,
- Real edgeInclusionProb, ifstream &initfile) {
- myInt i;
- X = data;
- n = n_obs;
- p = p_model;
- L = L_start;
- xi = new myInt[n];
- graphlist = new LPGraph[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- };
- pll = new Real[L];
-
- // alpha = 20;
- alpha = 1;
-
-#ifdef RANDOMSTART
- RandomStartAllXi(L);
- RandomStartAllG(L, edgeInclusionProb);
-#else
- ReadState(initfile);
-#endif
-
- plp = partitionlogPrior(L, xi, alpha);
- for (i = 0; i < L; i++) {
- pll[i] = cluster_k_loglikelihood(i, xi, graphlist[i]);
- }
-}
-
-DPmixGGM::DPmixGGM(State a) {
- myInt i;
-
- n = a->n;
- p = a->p;
- L = a->L;
- X = a->X;
- alpha = a->alpha;
- plp = a->plp;
- xi = new myInt[n];
- for (i = 0; i < n; i++) {
- xi[i] = a->xi[i];
- }
- graphlist = new LPGraph[L];
- pll = new Real[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- graphlist[i]->CopyGraph(a->graphlist[i]);
- pll[i] = a->pll[i];
- }
-}
-
-//-------- Destructor ------------------------------
-DPmixGGM::~DPmixGGM() {
- myInt i;
- delete[] xi;
- for (i = 0; i < L; i++) {
- delete graphlist[i];
- };
- delete[] graphlist;
- delete[] pll;
-}
-
-void DPmixGGM::RandomStartAllXi(myInt L_start) {
- for (myInt i = 0; i < n; i++) {
- xi[i] = (myInt)(L_start * gsl_ran_flat(rnd, 0.0, 1.0)) / 1;
- }
-}
-
-void DPmixGGM::RandomStartAllG(myInt L_start, Real edgeInclusionProb) {
- myInt i, k, l;
- bool temp;
-
- for (i = 0; i < L_start; i++) {
- for (k = 0; k < p - 1; k++) {
- graphlist[i]->Edge[k][k] = 0;
- for (l = k + 1; l < p; l++) {
- temp = (gsl_ran_flat(rnd, 0.0, 1.0) < edgeInclusionProb);
- graphlist[i]->Edge[k][l] = temp;
- graphlist[i]->Edge[l][k] = temp;
- }
- }
- graphlist[i]->Edge[p - 1][p - 1] = 0;
- TurnFillInGraph(graphlist[i]);
- graphlist[i]->GenerateAllCliques();
- }
-}
-
-void DPmixGGM::RandomStart(myInt L_start, Real edgeInclusionProb) {
- myInt i;
- for (i = 0; i < L; i++) {
- delete graphlist[i];
- };
- delete graphlist;
- delete pll;
- L = L_start;
- graphlist = new LPGraph[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- };
- pll = new Real[L];
-
- RandomStartAllXi(L);
- RandomStartAllG(L, edgeInclusionProb);
-
- plp = partitionlogPrior(L, xi, alpha);
- for (i = 0; i < L; i++) {
- pll[i] = cluster_k_loglikelihood(i, xi, graphlist[i]);
- }
-}
-
-void DPmixGGM::InformedStart(ifstream &initfile) {
- ReadState(initfile);
- plp = partitionlogPrior(L, xi, alpha);
- for (myInt i = 0; i < L; i++) {
- pll[i] = cluster_k_loglikelihood(i, xi, graphlist[i]);
- }
-}
-
-void DPmixGGM::ReadState(ifstream &initfile) {
- myInt i, j, k, l;
-
- for (l = 0; l < L; l++) {
- delete graphlist[l];
- }
- delete[] graphlist;
- delete[] pll;
-
- int supern, superp;
- initfile >> supern;
- initfile >> superp;
- initfile >> L; // cout << "supern = " << supern << " superp = " << superp << "
- // L = " << L << endl; fflush(stdout);
- for (i = 0; i < n; i++) {
- initfile >> xi[i];
- xi[i]--;
- };
- for (i = n; i < supern; i++) {
- initfile >> j;
- }
- graphlist = new LPGraph[L];
- pll = new Real[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- for (k = 0; k < p; k++) {
- for (l = 0; l < p; l++) {
- initfile >> graphlist[i]->Edge[k][l];
- };
- for (l = p; l < superp; l++) {
- initfile >> j;
- }
- }
- for (k = p; k < superp; k++) {
- for (l = 0; l < superp; l++) {
- initfile >> j;
- }
- }
- TurnFillInGraph(graphlist[i]);
- graphlist[i]->GenerateAllCliques();
- }
-
- // cout << "end of ReadState" << endl; fflush(stdout);
-}
-
-void DPmixGGM::WriteState(ofstream &out, int itr) {
- myInt i, j, l, q, r;
- Real score = plp;
- for (l = 0; l < L; l++) {
- score += pll[l];
- }
-
- out << L << " " << score << " " << itr << " ";
- for (i = 0; i < n; i++) {
- out << xi[i] << " ";
- }
- for (l = 0; l < L; l++) {
- for (q = 0; q < p - 1; q++) {
- for (r = q + 1; r < p; r++) {
- out << graphlist[l]->Edge[q][r] << " ";
- }
- }
- }
- for (l = 0; l < L; l++) {
- out << graphlist[l]->nCliques << " ";
- for (i = 0; i < (graphlist[l]->nCliques); i++) {
- out << graphlist[l]->CliquesDimens[i] << " ";
- for (j = 0; j < (graphlist[l]->CliquesDimens[i]); j++) {
- out << graphlist[l]->Cliques[i][j] << " ";
- }
- }
-
- out << graphlist[l]->nTreeEdges << " ";
- for (i = 0; i < (graphlist[l]->nTreeEdges); i++) {
- out << graphlist[l]->TreeEdgeA[i] << " " << graphlist[l]->TreeEdgeB[i]
- << " ";
- }
-
- out << (graphlist[l]->nSeparators) << " ";
- for (i = 0; i < (graphlist[l]->nSeparators); i++) {
- out << graphlist[l]->SeparatorsDimens[i] << " ";
- for (j = 0; j < (graphlist[l]->SeparatorsDimens[i]); j++) {
- out << graphlist[l]->Separators[i][j] << " ";
- }
- }
- }
- out << endl;
-}
-
-void DPmixGGM::CopyState(State a) {
- myInt i;
- myInt oldL = L;
-
- n = a->n;
- p = a->p;
- X = a->X;
- alpha = a->alpha;
- plp = a->plp;
- for (i = 0; i < n; i++) {
- xi[i] = a->xi[i];
- };
- if (L != a->L) {
- L = a->L;
- pll = new Real[L];
- for (i = 0; i < oldL; i++) {
- delete graphlist[i];
- };
- delete[] graphlist;
- graphlist = new LPGraph[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- }
- }
- for (i = 0; i < L; i++) {
- graphlist[i]->CopyGraph(a->graphlist[i]);
- pll[i] = a->pll[i];
- }
-}
-
-#ifdef JEFFREYS_PRIOR
-struct f_params {
- int n;
- int m;
- int k;
-};
-double f(double beta, void *params) {
- struct f_params *iparams = (f_params *)params;
- int k = iparams->k;
- int n = iparams->n;
- int m = iparams->m;
- double s;
- double sum = 0;
- for (int j = 1; j < m; j++) {
- s = beta + j;
- sum += ((Real)j) / (s * s);
- }
- return exp(lgamma(beta) - lgamma(beta + n) + lgamma(n + 1.0) +
- (k - 0.5) * log(beta) + log(sqrt(sum)));
-}
-
-Real DPmixGGM::lA(myInt the_n, myInt m, myInt k) {
- struct f_params params = {the_n, m, k};
- F.function = &f;
- F.params = ¶ms;
- double result, error;
- gsl_integration_qagiu(&F, 0.0, 1e-7, 1e-7, GSL_INTEGRATION_GRIDSIZE, w,
- &result, &error);
- return log(result);
-}
-#endif
-
-//--------------------------------------------------------------------
-// This returns the partition prior probability -- normalising constant ignored
-// assuming alpha is fixed maxL is maximum number of clusters, effective number
-// of clusters can be smaller
-Real DPmixGGM::partitionlogPrior(myInt maxL, myInt *thexi, Real thealpha) {
- myInt i, j;
- Real siz;
- Real pri = 0;
- myInt effectiveL = 0;
- for (i = 0; i < maxL; i++) {
- siz = 0.0;
- for (j = 0; j < n; j++) {
- if (thexi[j] == i)
- siz = siz + 1;
- };
- if (siz > 0) {
- pri += lgamma(siz);
- effectiveL++;
- }
- }
-
-#ifdef JEFFREYS_PRIOR
- struct f_params params = {n, n, effectiveL};
- F.function = &f;
- F.params = ¶ms;
- double result, error;
- gsl_integration_qagiu(&F, 0.0, 1e-7, 1e-7, GSL_INTEGRATION_GRIDSIZE, w,
- &result, &error);
- return pri + log(result);
-#else
- return pri + effectiveL * log(thealpha);
-#endif
-}
-
-// This returns the likelihood of the observations in cluster k
-Real DPmixGGM::cluster_k_loglikelihood(myInt k, myInt *thexi,
- LPGraph thegraph) {
- //-------- Variables, follows paper ---------------
- int i, j, ii;
- Real n0 = N0;
- Real *mu0 = new Real[p];
- for (i = 0; i < p; i++)
- mu0[i] = 0;
- Real *xbar = new Real[p];
- Real *mu_bar = new Real[p];
- Real *D_prior = new Real[p * (p + 1) / 2];
- Real *D_post = new Real[p * (p + 1) / 2];
- myInt n_sub = 0;
- Real J_G;
- Real Norm_terms;
-
- for (i = 0; i < n; i++) {
- if (thexi[i] == k)
- n_sub++;
- } // count the number in cluster k
- if (n_sub == 0) {
- delete[] D_post;
- delete[] D_prior;
- delete[] xbar;
- delete[] mu_bar;
- delete[] mu0;
- return (0);
- }
- // cout << "l = " << k << ", n_sub = " << n_sub << endl; fflush(stdout);
-
- //----------- Form the sufficient statistics ----------
- make_sub_means_and_cov(X, thexi, k, p, n, n_sub, xbar, D_post);
-
- //------------------ Update these ---------------------
- for (i = 0; i < p; i++)
- mu_bar[i] = (n_sub * xbar[i] + n0 * mu0[i]) / (n_sub + n0);
- for (i = 0; i < p * (p + 1) / 2; i++)
- D_prior[i] = 0;
- for (i = 0; i < p; i++)
- D_prior[i * (i + 1) / 2 + i] = 1;
- for (i = 0; i < p; i++)
- D_post[i * (i + 1) / 2 + i] += 1;
-
- //----------------- Factor in the mean for the matrix D_posterior
- //-------------
- for (i = 0; i < p; i++) {
- ii = i * (i + 1) / 2;
- for (j = 0; j <= i; j++) {
- D_post[ii + j] += (-(n_sub + n0) * mu_bar[i] * mu_bar[j] +
- n_sub * xbar[i] * xbar[j] + n0 * mu0[i] * mu0[j]);
- }
- }
-
- //------------ Calculate the score
- //---------------------------------------------
- J_G = j_g_decomposable(thegraph, D_prior, D_post, DELTA0, n_sub, 0);
- Norm_terms =
- -(Real(n_sub * p) / 2) * log_2_pi + Real(p) / 2 * log(n0 / (n_sub + n0));
-
- delete[] D_post;
- delete[] D_prior;
- delete[] xbar;
- delete[] mu_bar;
- delete[] mu0;
- return (Norm_terms + J_G);
-}
-
-//--------------------------------------------------------------------
-// This returns the predictive distribution score for observation k, if it were
-// to belong to cluster l thegraph is the graph associated with cluster l
-Real DPmixGGM::predictiveDistribution(myInt k, myInt l, myInt *thexi,
- LPGraph thegraph) {
- int i, j, ii;
- Real a;
- Real n0 = N0;
- Real *mu0 = new Real[p];
- Real *xbar = new Real[p];
- Real *mu_bar = new Real[p];
- Real *mu_tilde = new Real[p];
- for (i = 0; i < p; i++) {
- mu0[i] = 0;
- }
- Real *D_prior = new Real[p * (p + 1) / 2];
- Real *D_post = new Real[p * (p + 1) / 2];
- myInt n_sub = 0;
- Real coef;
- Real J_G;
- Real Norm_terms;
-
- for (i = 0; i < n; i++) {
- if (thexi[i] == l) {
- n_sub++;
- }
- } // count the number in this cluster
-#ifdef JEFFREYS_PRIOR
- if (l == L) {
- coef = lA(n, n, L + 1) - lA(n - 1, n, L);
- } else if (n_sub > 0) {
- coef = log(n_sub) + lA(n, n, L) - lA(n - 1, n, L);
- } else {
- coef = 0;
- }
-#else
- if (l == L) {
- coef = log(alpha);
- } else if (n_sub > 0) {
- coef = log(n_sub);
- } else {
- coef = 0;
- }
-#endif
-
- //----------- Form the sufficient statistics ----------
- make_sub_means_and_cov(X, thexi, l, p, n, n_sub, xbar, D_prior);
-
- //------------------ Update these ---------------------
- for (i = 0; i < p; i++)
- mu_bar[i] = (n_sub * xbar[i] + n0 * mu0[i]) / (n_sub + n0);
- for (i = 0; i < p; i++)
- D_prior[i * (i + 1) / 2 + i] += 1;
- //----------------- Factor in the mean for the matrix D_prior -------------
- for (i = 0; i < p; i++) {
- ii = i * (i + 1) / 2;
- for (j = 0; j <= i; j++) {
- D_prior[ii + j] += -(n_sub + n0) * mu_bar[i] * mu_bar[j] +
- n_sub * xbar[i] * xbar[j] + n0 * mu0[i] * mu0[j];
- }
- }
-
- //-------------- Get the posterior information ----------------------------
- for (i = 0; i < p * (p + 1) / 2; i++)
- D_post[i] = D_prior[i];
- for (i = 0; i < p; i++)
- mu_tilde[i] = (X[k * p + i] + (n_sub + n0) * mu_bar[i]) / (n_sub + n0 + 1);
- for (i = 0; i < p; i++) {
- ii = i * (i + 1) / 2;
- for (j = 0; j <= i; j++) {
- a = -(n_sub + 1 + n0) * mu_tilde[i] * mu_tilde[j] +
- X[k * p + i] * X[k * p + j] + (n_sub + n0) * mu_bar[i] * mu_bar[j];
- D_post[ii + j] += a;
- }
- }
-
- //------------ Calculate the score
- //---------------------------------------------
- J_G = j_g_decomposable(thegraph, D_prior, D_post, DELTA0 + n_sub, 1, 1);
- Norm_terms =
- -(p / 2) * log_2_pi + (p / 2) * log((n_sub + n0) / (n_sub + 1 + n0));
-
- delete[] D_post;
- delete[] D_prior;
- delete[] xbar;
- delete[] mu_bar;
- delete[] mu_tilde;
- delete[] mu0;
- return (coef + Norm_terms + J_G);
-}
\ No newline at end of file
diff --git a/src/sss-hip/DPmixGGM_Lists.cpp b/src/sss-hip/DPmixGGM_Lists.cpp
deleted file mode 100644
index 191ced7e5e..0000000000
--- a/src/sss-hip/DPmixGGM_Lists.cpp
+++ /dev/null
@@ -1,413 +0,0 @@
-#define LISTS_CPP
-#ifndef GRAPH_CPP
-#include "graph.cpp"
-#endif
-#ifndef GWISH_CPP
-#include "gwish.cpp"
-#endif
-#ifndef DPMIXGGM_CPP
-#include "DPmixGGM.cpp"
-#endif
-
-typedef class DPmixGGMlist *List;
-
-class DPmixGGMlist {
- // variables
-public:
- int n;
- int p;
- myInt M;
-
- myInt *L_list;
- myInt *xi_list;
- myInt *edge_list;
- Real *score_list;
-
- myInt ind_minOfListScores;
- Real minOfListScores;
-
- // functions
-public:
- DPmixGGMlist(myInt size, int n, int p);
- // DPmixGGMlist (myInt size, State a);
- ~DPmixGGMlist();
-
- void FlushList(State a);
- void WriteList(ofstream &out);
- void ReadFromList(State a, myInt m);
-
- void UpdateList(myInt L, myInt *xi, LPGraph *graphlist, Real score);
- void UpdateList(State a);
-
- LPGraph ProposeGraph(myInt NN_xi, myInt *which_xi, Real sFactor);
-};
-
-DPmixGGMlist::DPmixGGMlist(myInt size, int in_n, int in_p) {
- int i;
- n = in_n;
- p = in_p;
- M = size;
- int ee = p * (p - 1) / 2;
-
- L_list = new myInt[M];
- for (i = 0; i < M; i++) {
- L_list[i] = -1;
- }
- xi_list = new myInt[M * n];
- for (i = 0; i < M * n; i++) {
- xi_list[i] = -1;
- }
- edge_list = new myInt[M * n * ee];
- for (i = 0; i < (M * n * ee); i++) {
- edge_list[i] = -1;
- }
- score_list = new Real[M];
- for (i = 0; i < M; i++) {
- score_list[i] = NEG_INF;
- }
-
- ind_minOfListScores = 0;
- minOfListScores = NEG_INF;
-}
-
-DPmixGGMlist::~DPmixGGMlist() {
- delete[] L_list;
- delete[] xi_list;
- delete[] edge_list;
- delete[] score_list;
-}
-
-void DPmixGGMlist::FlushList(State a) {
- int i;
- int ee = p * (p - 1) / 2;
-
- for (i = 0; i < M; i++) {
- L_list[i] = -1;
- }
- for (i = 0; i < M * n; i++) {
- xi_list[i] = -1;
- }
- for (i = 0; i < (M * n * ee); i++) {
- edge_list[i] = -1;
- }
- for (i = 0; i < M; i++) {
- score_list[i] = NEG_INF;
- }
-
- int q, r, l;
- L_list[0] = a->L;
- for (i = 0; i < n; i++) {
- xi_list[i] = a->xi[i];
- }
- for (i = 0; i < n; i++) {
- l = 0;
- for (q = 0; q < p - 1; q++) {
- for (r = q + 1; r < p; r++) {
- edge_list[i * ee + l] = a->graphlist[a->xi[i]]->Edge[q][r];
- l++;
- }
- }
- }
- score_list[0] = a->plp;
- for (i = 0; i < (a->L); i++) {
- score_list[0] += a->pll[i];
- }
-
- ind_minOfListScores = 1;
- minOfListScores = NEG_INF;
-}
-
-void DPmixGGMlist::WriteList(ofstream &out) {
- int i, j, l, L, q, r, t;
- int ee = p * (p - 1) / 2;
-
- for (t = 0; t < M; t++) {
- if (L_list[t] == -1) {
- continue;
- }
-
- L = L_list[t];
- out << L << " " << score_list[t] << " ";
- for (i = 0; i < n; i++) {
- out << xi_list[t * n + i] << " ";
- }
-
- LPGraph *graphlist = new LPGraph[L];
- for (l = 0; l < L; l++) {
- graphlist[l] = new Graph();
- graphlist[l]->InitGraph(p);
- }
- for (l = 0; l < L; l++) {
- for (i = 0; i < n; i++) {
- if (xi_list[t * n + i] == l)
- break;
- }
- j = 0;
- for (q = 0; q < p - 1; q++) {
- graphlist[l]->Edge[q][q] = 0;
- for (r = q + 1; r < p; r++) {
- graphlist[l]->Edge[q][r] = edge_list[t * n * ee + i * ee + j];
- graphlist[l]->Edge[r][q] = graphlist[l]->Edge[q][r];
- out << graphlist[l]->Edge[q][r] << " ";
- j++;
- }
- }
- graphlist[l]->Edge[p - 1][p - 1] = 0;
- graphlist[l]->GenerateAllCliques();
- }
-
- for (l = 0; l < L; l++) {
- out << graphlist[l]->nCliques << " ";
- for (i = 0; i < (graphlist[l]->nCliques); i++) {
- out << graphlist[l]->CliquesDimens[i] << " ";
- for (j = 0; j < (graphlist[l]->CliquesDimens[i]); j++) {
- out << graphlist[l]->Cliques[i][j] << " ";
- }
- }
-
- out << graphlist[l]->nTreeEdges << " ";
- for (i = 0; i < (graphlist[l]->nTreeEdges); i++) {
- out << graphlist[l]->TreeEdgeA[i] << " " << graphlist[l]->TreeEdgeB[i]
- << " ";
- }
-
- out << graphlist[l]->nSeparators << " ";
- for (i = 0; i < (graphlist[l]->nSeparators); i++) {
- out << graphlist[l]->SeparatorsDimens[i] << " ";
- for (j = 0; j < (graphlist[l]->SeparatorsDimens[i]); j++) {
- out << graphlist[l]->Separators[i][j] << " ";
- }
- }
- }
-
- out << endl;
- for (l = 0; l < L; l++) {
- delete graphlist[l];
- };
- delete[] graphlist;
- }
-}
-
-void DPmixGGMlist::ReadFromList(State a, myInt m) {
- int i, j, k, l, L, t;
- int ee = p * (p - 1) / 2;
-
- for (l = 0; l < (a->L); l++) {
- delete a->graphlist[l];
- };
- delete[] a->graphlist;
- delete[] a->pll; // delete previous memory spaces
- L = L_list[m];
- a->L = L;
- a->pll = new Real[L];
- a->graphlist = new LPGraph[L];
- for (l = 0; l < L; l++) {
- a->graphlist[l] = new Graph();
- a->graphlist[l]->InitGraph(p);
- }
- for (i = 0; i < n; i++) {
- a->xi[i] = xi_list[m * n + i];
- }
- for (l = 0; l < L; l++) {
- for (k = 0; k < n; k++) {
- if (a->xi[k] == l) {
- break;
- }
- };
- t = 0;
- for (i = 0; i < (p - 1); i++) {
- for (j = (i + 1); j < p; j++) {
- a->graphlist[l]->Edge[i][j] = edge_list[m * n * ee + k * ee + t];
- a->graphlist[l]->Edge[j][i] = a->graphlist[l]->Edge[i][j];
- t++;
- }
- }
- a->graphlist[l]->GenerateAllCliques();
- }
-}
-
-// OVERWRITES A GIVEN LIST OF MODELS WITH A GIVEN MODEL IF IT HAS HIGHER SCORE
-void DPmixGGMlist::UpdateList(myInt L, myInt *xi, LPGraph *graphlist,
- Real score) {
- int i, j, k, l, q, r, t;
- int ee = p * (p - 1) / 2;
-
- // check if the candidate model deserve to be in the list of models -- if so,
- // do the necessary
- bool flag;
- int count;
- myInt *xi_new = new myInt[n];
- myInt *index_set = new myInt[L];
- for (i = 0; i < n; i++) {
- xi_new[i] = -1;
- }
-
- if (score > minOfListScores) {
- // reindexing (lexicographically) the incoming model -- if there is a
- // redundant index, getting rid of it
- myInt xi_temp = xi[0];
- index_set[0] = xi_temp;
- for (i = 0; i < n; i++) {
- if (xi[i] == xi_temp) {
- xi_new[i] = 0;
- }
- }
-
- j = 1;
- l = 0;
- while (j < n) {
- k = j;
- l++;
- while (k < n) // finding a new cluster index
- {
- flag = 1;
- for (t = 0; t < l; t++) {
- if (xi[k] == index_set[t]) {
- flag = 0;
- break;
- }
- }
- if (flag) {
- xi_temp = xi[k];
- index_set[l] = xi_temp;
- for (i = 0; i < n; i++) {
- if (xi[i] == xi_temp) {
- xi_new[i] = l;
- }
- }
- break;
- } else {
- k++;
- }
- }
- j = k + 1;
- }
- k = l;
-
- // check if this model is already in the list of models
- count = 0;
- for (t = 0; t < M; t++) {
- flag = 0;
- if (L_list[t] != k) {
- flag = 1;
- }
-
- if (flag == 0) {
- for (i = 0; i < n; i++) {
- if (xi_list[t * n + i] != xi_new[i]) {
- flag = 1;
- break;
- }
- }
- };
-
- if (flag == 0) {
- for (l = 0; l < k; l++) {
- for (i = 0; i < n; i++) {
- if (xi_list[t * n + i] == l)
- break;
- } // finiding an observation with cluster index l
-
- j = 0;
- for (q = 0; q < p - 1; q++) {
- for (r = q + 1; r < p; r++) {
- if (edge_list[t * n * ee + i * ee + j] !=
- graphlist[xi[i]]->Edge[q][r]) {
- flag = 1;
- break;
- };
- j++;
- }
- if (flag)
- break;
- }
- if (flag)
- break;
- }
- }
-
- if (flag)
- count++;
- }
-
- // if the model is not in list, substitute the lowest scoring model with the
- // candidate
- if (count == M) {
- L_list[ind_minOfListScores] = k; // cout << "effective L = " << k << endl;
- for (i = 0; i < n; i++) {
- xi_list[ind_minOfListScores * n + i] = xi_new[i];
- t = 0;
- for (q = 0; q < p - 1; q++) {
- for (r = q + 1; r < p; r++) {
- edge_list[ind_minOfListScores * n * ee + i * ee + t] =
- graphlist[xi[i]]->Edge[q][r];
- t++;
- }
- }
- }
- score_list[ind_minOfListScores] = score;
-
- // searching for minimum scoring model in the list
- ind_minOfListScores = 0;
- minOfListScores = score_list[0];
- for (t = 1; t < M; t++) {
- if (score_list[t] < minOfListScores) {
- minOfListScores = score_list[t];
- ind_minOfListScores = t;
- }
- }
- }
- }
-
- delete[] xi_new;
- delete[] index_set;
-}
-
-void DPmixGGMlist::UpdateList(State a) {
- Real score = a->plp;
- for (myInt i = 0; i < a->L; i++) {
- score += a->pll[i];
- }
- UpdateList(a->L, a->xi, a->graphlist, score);
-}
-
-// PROPOSES A DECOMPOSABLE GRAPH
-LPGraph DPmixGGMlist::ProposeGraph(myInt NN_xi, myInt *which_xi, Real sFactor) {
- int i, j, k, l, m, t;
- Real a, s;
- bool temp;
- int ee = p * (p - 1) / 2;
-
- LPGraph newgraph = new Graph();
- newgraph->InitGraph(p);
-
- // feature extraction
- l = 0;
- for (i = 0; i < p - 1; i++) {
- newgraph->Edge[i][i] = 0;
- for (j = i + 1; j < p; j++) {
- s = NN_xi * M * sFactor;
- a = s;
-
- for (m = 0; m < M; m++) {
- for (k = 0; k < NN_xi; k++) {
- t = edge_list[m * n * ee + which_xi[k] * ee + l];
- if (t != -1) {
- a += Real((bool)t);
- s += 1.0;
- }
- }
- }
- l++;
-
- temp = (gsl_ran_flat(rnd, 0.0, 1.0) < a / s);
- newgraph->Edge[i][j] = temp;
- newgraph->Edge[j][i] = temp;
- }
- }
- newgraph->Edge[p - 1][p - 1] = 0;
- TurnFillInGraph(newgraph);
- newgraph->GenerateAllCliques();
-
- return newgraph;
-}
\ No newline at end of file
diff --git a/src/sss-hip/Makefile b/src/sss-hip/Makefile
index b55419e4c2..c49d14bbf0 100755
--- a/src/sss-hip/Makefile
+++ b/src/sss-hip/Makefile
@@ -6,7 +6,6 @@
CC = hipcc
OPTIMIZE = yes
DEBUG = no
-GSL = /path/to/GSL
LAUNCHER =
#===============================================================================
@@ -24,10 +23,10 @@ obj = $(source:.cu=.o)
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I$(GSL)/include
+CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../sss-cuda
# Linker Flags
-LDFLAGS = -L$(GSL)/lib -lgsl
+LDFLAGS =
# Debug Flags
ifeq ($(DEBUG),yes)
@@ -46,12 +45,16 @@ endif
$(program): $(obj) Makefile
$(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS)
-%.o: %.cu DPmixGGM.cpp DPmixGGM_Lists.cpp DPmixGGM_SSSmoves.cpp \
- kernels.cu graph.cpp gwish.cpp utilities.cpp Makefile
+SHARED_SRC = ../sss-cuda/DPmixGGM.cpp ../sss-cuda/DPmixGGM_Lists.cpp \
+ ../sss-cuda/graph.cpp ../sss-cuda/gwish.cpp \
+ ../sss-cuda/utilities.cpp ../sss-cuda/graph.h \
+ ../sss-cuda/gsl_compat.h ../sss-cuda/kernels.cu
+
+%.o: %.cu DPmixGGM_SSSmoves.cpp $(SHARED_SRC) Makefile
$(CC) $(CFLAGS) -c $< -o $@
clean:
rm -rf $(program) $(obj) RES/f9_n150_p50_modes_GPU.txt
run: $(program)
- $(LAUNCHER) GSL_RNG_SEED=123 ./$(program) f9_n150_p50
+ $(LAUNCHER) ./$(program) f9_n150_p50
diff --git a/src/sss-hip/README.md b/src/sss-hip/README.md
index d020e7e2f9..d313bd8d90 100644
--- a/src/sss-hip/README.md
+++ b/src/sss-hip/README.md
@@ -1,22 +1,28 @@
DPmixGGM
========
-This folder contains source codes for the "GPU-powered Stochastic Shotgun Search for Dirichlet proces mixtures of Gaussian Graphical Models"
+This folder contains source codes for the "GPU-powered Stochastic Shotgun Search for Dirichlet proces mixtures of Gaussian Graphical Models"
by Chiranjit Mukherjee and Abel Rodriguez.
-
-The "DPmixGGM_SSS_main.cpp" file contains tuning parameters for the algorithm, as elaborated below:
+
+The "DPmixGGM_SSS_main.cpp" file contains tuning parameters for the algorithm, as elaborated below:
1. Run the SSS
-2. Run GPU/CPU versions of the SSS by enabling / disabling the macro CUDA.
-3. Specify maximum number of mixture components that the model should accommodate (for pre-allocation of memory).
-4. Set SSS runtime parameters C, D, R, S, M, g, h, f, t.
-5. Set SSS number of chain parameters. User needs to provide at least one initial point.
-7. Set hyperparameters of for the prior on (mu, K | G) with N0, DELTA0.
-
+2. Run GPU/CPU versions of the SSS by enabling / disabling the macro CUDA.
+3. Specify maximum number of mixture components that the model should accommodate (for pre-allocation of memory).
+4. Set SSS runtime parameters C, D, R, S, M, g, h, f, t.
+5. Set SSS number of chain parameters. User needs to provide at least one initial point.
+7. Set hyperparameters of for the prior on (mu, K | G) with N0, DELTA0.
+
Complie source codes using the "make" command and run with "./main f9_n150_p50" command.
-
-The program expects an input-data file (e.g. f9_n150_p50) in the DATA/ folder and at least one initialization point (e.g. f9_n150_p50_init1).
-The input-data file should specify n and p in the first row and then provide n rows of length p. The initial point data-file should specify n, p
-and L of the initial model configuration in the first row and xi-indices of the initial point in the second row. Subsequent L rows specify
-G_l (l=1:L).
-
+
+The program expects an input-data file (e.g. f9_n150_p50) in the DATA/ folder and at least one initialization point (e.g. f9_n150_p50_init1).
+The input-data file should specify n and p in the first row and then provide n rows of length p. The initial point data-file should specify n, p
+and L of the initial model configuration in the first row and xi-indices of the initial point in the second row. Subsequent L rows specify
+G_l (l=1:L).
+
A list of highest-score models is stored in folder RES/.
+
+Dependencies
+------------
+None beyond a C++17 compiler. The uniform random number generator (MT19937) and
+the adaptive Gauss-Kronrod quadrature the sampler needs are implemented in
+`gsl_compat.h` in this folder.
diff --git a/src/sss-hip/graph.cpp b/src/sss-hip/graph.cpp
deleted file mode 100644
index 55e75deb02..0000000000
--- a/src/sss-hip/graph.cpp
+++ /dev/null
@@ -1,973 +0,0 @@
-#include
-#include
-#include
-#include
-
-#define GRAPH_CPP
-#ifndef GRAPH_H
-#include "graph.h"
-#endif
-#ifndef GWISH_CPP
-#include "gwish.cpp"
-#endif
-
-// class Graph::Begins
-Graph::Graph() {
- nVertices = 0;
- Edge = NULL;
- Labels = NULL;
- nLabels = 0;
- Cliques = NULL;
- CliquesDimens = NULL;
- nCliques = 0;
- TreeEdgeA = NULL;
- TreeEdgeB = NULL;
- nTreeEdges = 0;
- Separators = NULL;
- SeparatorsDimens = NULL;
- nSeparators = 0;
- localord = NULL;
- return;
-}
-
-Graph::Graph(LPGraph InitialGraph) {
- nVertices = 0;
- Edge = NULL;
- Labels = NULL;
- nLabels = 0;
- Cliques = NULL;
- CliquesDimens = NULL;
- nCliques = 0;
- TreeEdgeA = NULL;
- TreeEdgeB = NULL;
- nTreeEdges = 0;
- Separators = NULL;
- SeparatorsDimens = NULL;
- nSeparators = 0;
- localord = NULL;
-
- ///////////////////////////////////////
- myInt i, j;
- InitGraph(InitialGraph->nVertices);
- for (i = 0; i < nVertices; i++) {
- for (j = 0; j < nVertices; j++) {
- Edge[i][j] = InitialGraph->Edge[i][j];
- }
- }
-
- return;
-}
-
-Graph::~Graph() {
- myInt i; // cout << "-> "; fflush(stdout);
-
- for (i = 0; i < nVertices; i++) {
- delete[] Edge[i];
- Edge[i] = NULL;
- }
- delete[] Edge;
- Edge = NULL;
-
- delete[] Labels;
- Labels = NULL;
-
- delete[] Cliques[0]; // for(i=0; inVertices;
- for (i = 0; i < n; i++) {
- for (j = 0; j < n; j++) {
- Edge[i][j] = G->Edge[i][j];
- }
- }
- nLabels = G->nLabels;
- for (i = 0; i < n; i++) {
- Labels[i] = G->Labels[i];
- }
- nCliques = G->nCliques;
- for (i = 0; i < n; i++) {
- CliquesDimens[i] = G->CliquesDimens[i];
- for (j = 0; j < n; j++) {
- Cliques[i][j] = G->Cliques[i][j];
- }
- }
- nTreeEdges = G->nTreeEdges;
- for (i = 0; i < n; i++) {
- TreeEdgeA[i] = G->TreeEdgeA[i];
- TreeEdgeB[i] = G->TreeEdgeB[i];
- }
- nSeparators = G->nSeparators;
- for (i = 0; i < n; i++) {
- SeparatorsDimens[i] = G->SeparatorsDimens[i];
- for (j = 0; j < n; j++) {
- Separators[i][j] = G->Separators[i][j];
- }
- }
- for (i = 0; i < n; i++) {
- localord[i] = G->localord[i];
- }
-}
-
-void Graph::GenerateCliques(myInt label) {
- myInt i, j, k, p, r;
- myInt n = nVertices;
- myInt *clique = new myInt[nVertices];
- memset(clique, 0, nVertices * sizeof(myInt));
-
- myInt countA, countB;
- myInt *listA = new myInt[n];
- myInt *listB2 = new myInt[n];
-
- // clean memory
- memset(localord, 0, nVertices * sizeof(myInt));
-
- myInt v, vk;
- myInt PrevCard = 0;
- myInt NewCard;
- myInt s = nCliques - 1; // cout << "s = " << s << endl;
-
- countA = 0;
- for (i = 0; i < n; i++) {
- if (Labels[i] == label) {
- listA[countA] = i;
- countA++;
- }
- };
- countB = 0;
-
- for (i = n; i >= 0; i--) {
- NewCard = -1;
-
- // choose a vertex v...
- for (j = 0; j < countA; j++) {
- myInt maxj = 0;
- for (r = 0; r < countB; r++) {
- if (Edge[listA[j]][listB2[r]]) {
- maxj++;
- }
- }
- if (maxj > NewCard) {
- v = listA[j];
- NewCard = maxj;
- }
- }
-
- // printf("i=%d, NewCard=%d, PrevCard=%d countA=%d,
- // countB=%d\n",i,NewCard,PrevCard,countA,countB);
-
- if (NewCard == -1) {
- break;
- }
-
- localord[v] = i;
- if (NewCard <= PrevCard) { // begin new clique
- s++;
-
- for (r = 0; r < countB; r++) {
- if (Edge[v][listB2[r]]) {
- Cliques[s][CliquesDimens[s]] = listB2[r];
- CliquesDimens[s]++;
- }
- }
-
- if (NewCard != 0) { // get edge to parent
- vk = Cliques[s][0];
- k = localord[vk]; // cout << "(" << clique[Cliques[s][0]] << " ";
- for (r = 1; r < CliquesDimens[s]; r++) {
- if (localord[Cliques[s][r]] < k) {
- vk = Cliques[s][r];
- k = localord[vk];
- }; // cout << clique[Cliques[s][r]] << " ";
- }
- // cout << "| ";
- // for(r=0; rIsDecomposable())
- return;
-
- myInt v1 = gfill->SearchVertex();
- // add edges to Def(Adj(x)) so that Adj(x) becomes a clique
- for (u = 0; u < gfill->nVertices; u++) {
- if (gfill->Edge[v1][u] == 1) {
- for (v = u + 1; v < gfill->nVertices; v++) {
- if ((gfill->Edge[v1][v] == 1) && (gfill->Edge[u][v] == 0)) {
- gfill->Edge[v][u] = gfill->Edge[u][v] = 1;
- }
- }
- }
- }
- LPEliminationGraph egraph = new EliminationGraph(graph, v1);
- for (i = 1; i < graph->nVertices - 1; i++) {
- v1 = egraph->SearchVertex();
- for (u = 0; u < egraph->nVertices; u++) {
- if (egraph->Eliminated[u])
- continue;
- if (egraph->Edge[v1][u] == 1) {
- for (v = u + 1; v < egraph->nVertices; v++) {
- if (egraph->Eliminated[v])
- continue;
- if ((egraph->Edge[v1][v] == 1) && (egraph->Edge[u][v] == 0)) {
- gfill->Edge[v][u] = gfill->Edge[u][v] = 1;
- // these are the edges that are added to the initial graph
- }
- }
- }
- }
- egraph->EliminateVertex(v1);
- }
- delete egraph;
- return;
-}
-
-// -----------------------------------------------------------------------------------------------
-// Based on Scott & Carvalho 2008
-
-bool Graph::CanAddEdge(myInt a, myInt b) {
- myInt i, j, k;
- bool canadd = 0;
- int n = (int)nVertices;
- myInt nS = 0, pR, pT;
- myInt *R = new myInt[nTreeEdges];
- myInt *T = new myInt[nTreeEdges];
- myInt *S = new myInt[n];
- myInt common_parent = 0;
- myInt contain_a, contain_b, contain_S;
-
- if (Labels[a] != Labels[b]) {
- canadd = 1;
- } // else ..
- else {
-
- nS = 0;
- for (j = 0; j < n; j++) {
- if (Edge[a][j] && Edge[b][j]) {
- S[nS] = j;
- nS++;
- }
- }
- if (nS == 0) {
- canadd = 0;
- } // else ....
- else { // HERE;
-
- // find higest-indexed cliques containing (a & S) and (b & S)
- myInt aSi = -1, bSi = -1;
- for (i = 0; i < nCliques; i++) {
- contain_a = 0;
- contain_b = 0;
- contain_S = 0;
- for (j = 0; j < CliquesDimens[i]; j++) {
- if (Cliques[i][j] == a) {
- contain_a = 1;
- };
- if (Cliques[i][j] == b) {
- contain_b = 1;
- }
- for (k = 0; k < nS; k++) {
- if (Cliques[i][j] == S[k]) {
- contain_S++;
- }
- }
- }
- if (contain_a && (contain_S == nS)) {
- aSi = i;
- }
- if (contain_b && (contain_S == nS)) {
- bSi = i;
- }
- }
-
- // find the path from aSi to root
- R[0] = -1;
- pR = -1;
- for (i = 0; i < nTreeEdges; i++) {
- if (TreeEdgeA[i] == aSi) {
- R[0] = i;
- pR = 0;
- break;
- }
- }
- for (i = R[0]; i >= 0; i--) {
- if (TreeEdgeA[i] == TreeEdgeB[R[pR]]) {
- pR++;
- R[pR] = i;
- }
- }
-
- // find the path from bSi to root
- T[0] = -1;
- pT = -1;
- for (i = 0; i < nTreeEdges; i++) {
- if (TreeEdgeA[i] == bSi) {
- T[0] = i;
- pT = 0;
- break;
- }
- }
- for (i = T[0]; i >= 0; i--) {
- if (TreeEdgeA[i] == TreeEdgeB[T[pT]]) {
- pT++;
- T[pT] = i;
- }
- }
-
- // find the branching point
- k = (pR < pT) ? pR : pT; // min(pR,pT);
- for (i = 0; i <= k; i++) {
- if (TreeEdgeB[R[pR - i]] == TreeEdgeB[T[pT - i]]) {
- common_parent = i;
- } else {
- break;
- }
- }
- if (k != -1) {
- if (TreeEdgeA[R[pR - common_parent]] ==
- TreeEdgeA[T[pT - common_parent]]) {
- common_parent++;
- }
- }
-
- // check if S is in the path from R[0] to common_parent
- for (i = 0; i <= (pR - common_parent); i++) {
- if (SeparatorsDimens[R[i]] == nS) {
- contain_S = 0;
- for (j = 0; j < SeparatorsDimens[R[i]]; j++) {
- for (k = 0; k < nS; k++) {
- if (Separators[R[i]][j] == S[k]) {
- contain_S++;
- }
- }
- }
- if (contain_S == nS) {
- canadd = 1;
- }
- }
- if (canadd) {
- break;
- }
- }
-
- // else: check if S is in the path from T[0] to common_parent
- if (!canadd) {
- for (i = 0; i <= (pT - common_parent); i++) {
- if (SeparatorsDimens[T[i]] == nS) {
- contain_S = 0;
- for (j = 0; j < SeparatorsDimens[T[i]]; j++) {
- for (k = 0; k < nS; k++) {
- if (Separators[T[i]][j] == S[k]) {
- contain_S++;
- }
- }
- }
- if (contain_S == nS) {
- canadd = 1;
- }
- }
- if (canadd) {
- break;
- }
- }
- }
- }
- }
-
- delete[] R;
- delete[] T;
- delete[] S;
-
- return (canadd);
-}
-
-Real Graph::ScoreAddEdge(myInt a, myInt b, Real *D_prior, Real *D_post,
- myInt delta, myInt n_sub, Real score, int nEdges) {
- myInt j;
- int n = (int)nVertices;
- myInt nS = 0;
- myInt *S = new myInt[n];
-
- nS = 0;
- for (j = 0; j < n; j++) {
- if (Edge[a][j] && Edge[b][j]) {
- S[nS] = j;
- nS++;
- }
- }
-
- {
- myInt *C = new myInt[n];
- myInt nC;
- Real *sub_D = new Real[n * (n + 1) / 2];
- Real cScore;
-
- nC = nS + 2;
- for (j = 0; j < nS; j++) {
- C[j] = S[j];
- };
- C[nS] = a;
- C[nS + 1] = b;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(n, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score += cScore;
-
- nC = nS + 1;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(n, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score -= cScore;
-
- nC = nS + 1;
- C[nS] = b;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(n, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score -= cScore;
-
- nC = nS;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(n, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score += cScore; // cout << nS << "->";
-
- delete[] S;
- delete[] C;
- delete[] sub_D;
- return score;
- }
-}
-
-myInt Graph::CanDeleteEdge(myInt a, myInt b) {
- myInt i, j;
- myBool contain_a, contain_b;
- myInt count = 0, which_ab;
-
- for (i = 0; i < nCliques; i++) {
- contain_a = 0;
- contain_b = 0;
- for (j = 0; j < CliquesDimens[i]; j++) {
- if (Cliques[i][j] == a) {
- contain_a = 1;
- };
- if (Cliques[i][j] == b) {
- contain_b = 1;
- }
- }
- if (contain_a && contain_b) {
- which_ab = i;
- count++;
- }
- }
-
- if (count == 1) {
- return (which_ab);
- } else {
- return (-1);
- }
-}
-
-Real Graph::ScoreDeleteEdge(myInt a, myInt b, myInt which_ab, Real *D_prior,
- Real *D_post, myInt delta, myInt n_sub, Real score,
- int nEdges) {
- myInt j;
-
- {
- myInt p = nVertices;
- myInt *C = new myInt[p];
- myInt nC;
- Real *sub_D = new Real[p * (p + 1) / 2];
- Real cScore;
-
- nC = 0;
- for (j = 0; j < CliquesDimens[which_ab]; j++) {
- if ((Cliques[which_ab][j] != a) && (Cliques[which_ab][j] != b)) {
- C[nC] = Cliques[which_ab][j];
- nC++;
- }
- };
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(p, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score -= cScore;
-
- nC = CliquesDimens[which_ab] - 1;
- C[nC - 1] = a;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(p, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score += cScore;
-
- nC = CliquesDimens[which_ab] - 1;
- C[nC - 1] = b;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(p, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score += cScore;
-
- nC = CliquesDimens[which_ab];
- C[nC - 2] = a;
- C[nC - 1] = b;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(p, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score -= cScore;
-
- delete[] C;
- delete[] sub_D;
- return score;
- }
-}
\ No newline at end of file
diff --git a/src/sss-hip/graph.h b/src/sss-hip/graph.h
deleted file mode 100644
index f5cda3b912..0000000000
--- a/src/sss-hip/graph.h
+++ /dev/null
@@ -1,117 +0,0 @@
-#define GRAPH_H
-
-typedef class Graph *LPGraph;
-
-class Graph {
-public:
- // data members
- myInt nVertices;
- myInt *d_nVertices; // number of vertices in the graph
- myInt **Edge;
- int *d_EdgeField; // matrix containing the edges of the graph
-
- myInt *Labels;
- myInt *d_Labels; // identifies the connected components of the graph
- myInt nLabels;
- myInt *d_nLabels; // number of labels or connected components
- myInt **Cliques; // storage for cliques
- myInt *CliquesDimens; // number of vertices in each clique
- myInt nCliques;
- myInt *d_nCliques; // number of cliques
-
-public:
- myInt *TreeEdgeA; // edges of the clique tree
- myInt *TreeEdgeB;
- myInt nTreeEdges; // number of edges in the generated clique tree
-public:
- myInt **Separators; // storage for separators
- myInt *SeparatorsDimens;
- myInt nSeparators;
- // private:
- myInt *localord;
-
- // methods
-public:
- Graph(); // constructor
- Graph(LPGraph InitialGraph); // constructor
- ~Graph(); // destructor
-public:
- myInt SearchVertex(); // identifies the next vertex to be eliminated
- void FlipEdge(myInt which); // flips an edge on a graph which is a myInteger
- // between 0 and the total number of edges
-
-public:
- // the MSS (Minimal Sufficient Statistics) are the maximal cliques for our
- // graph
- void InitGraph(myInt n);
- void CopyGraph(LPGraph G);
- void GenerateCliques(myInt label);
- myInt CheckCliques(myInt start,
- myInt end); // checks whether each generated component is
- // complete in the given graph
- myInt IsClique(
- myInt *vect,
- myInt nvect); // checks if the vertices in vect form a clique in our graph
-
- void GenerateSeparators();
- void AttachLabel(myInt v, myInt label);
- void GenerateLabels();
- myInt GenerateAllCliques();
- myInt IsDecomposable();
- myInt IfDecomposable();
-
- myInt CanDeleteEdge(myInt a, myInt b);
- bool CanAddEdge(myInt a, myInt b);
-
- Real ScoreDeleteEdge(myInt a, myInt b, myInt which_ab, Real *D_prior,
- Real *D_post, myInt delta, myInt n_sub, Real score,
- int nEdges);
- Real ScoreAddEdge(myInt a, myInt b, Real *D_prior, Real *D_post, myInt delta,
- myInt n_sub, Real score, int nEdges);
-};
-
-//////////////////////////////////////////////////////////////////////
-
-typedef class SectionGraph *LPSectionGraph;
-
-class SectionGraph : public Graph {
-public:
- myInt *
- Eliminated; // shows which vertices were eliminated from the initial graph
- myInt nEliminated; // number of vertices we eliminated
-
- // methods
-public:
- SectionGraph(LPGraph InitialGraph, myInt *velim); // constructor
- ~SectionGraph(); // destructor
-
-public:
- myInt IsChain(myInt u, myInt v); // see if there is a chain between u and v
- // or, equivalently, checks if u and v are in
- // the same connected component
-};
-
-////////////////////////////////////////////////////////////////////////
-
-typedef class EliminationGraph *LPEliminationGraph;
-
-class EliminationGraph : public Graph {
-public:
- myInt *
- Eliminated; // shows which vertices were eliminated from the initial graph
- myInt nEliminated; // number of vertices we eliminated
-
- // methods
-public:
- EliminationGraph(LPGraph InitialGraph, myInt vertex); // constructor
- ~EliminationGraph(); // destructor
-public:
- myInt SearchVertex(); // identify a vertex to be eliminated
-public:
- void EliminateVertex(myInt x); // eliminates an extra vertex
-};
-
-//////////////////////////////////////////////////////////////////////////
-
-// constructs the minimum fill-in graph for a nondecomposable graph
-void TurnFillInGraph(LPGraph graph);
\ No newline at end of file
diff --git a/src/sss-hip/gwish.cpp b/src/sss-hip/gwish.cpp
deleted file mode 100644
index 343adf69ef..0000000000
--- a/src/sss-hip/gwish.cpp
+++ /dev/null
@@ -1,167 +0,0 @@
-// gwish.cpp: This is a collection of device functions that manipulate graph
-// objects according to the G-Wishart distribution. Note: this depends on the
-// graph.cpp library Chiranjit Mukherjee : chiranjit@soe.ucsc.edu -- based on
-// Alex Lenkoski : lenkoski@stat.washington.edu
-
-#define GWISH_CPP
-#ifndef GRAPH_H
-#include "graph.h"
-#endif
-
-void log_det(int p, Real *A, Real *result) {
- int i, j, k;
- Real temp;
- *result = 0;
- Real br;
-
- for (i = 0; i < p; i++) {
- for (j = i; j < p; j++) {
- temp = A[j * (j + 1) / 2 + i];
- for (k = i; k > 0; k--) {
- temp = temp - A[j * (j + 1) / 2 + k - 1] * A[i * (i + 1) / 2 + k - 1];
- };
- if (i == j) {
- if (temp <= 0.0) {
- *result = NEG_INF;
- return;
- } else {
- br = sqrt(temp);
- }
- }
- A[j * (j + 1) / 2 + i] = temp / br;
- }
- }
-
-#ifdef ISFLOAT
- for (i = 0; i < p; i++) {
- *result += logf(A[i * (i + 1) / 2 + i]);
- };
- *result = 2 * (*result);
-#else
- for (i = 0; i < p; i++) {
- *result += log(A[i * (i + 1) / 2 + i]);
- };
- *result = 2 * (*result);
-#endif
-}
-
-// Computes the normalizing constant of a G-Wishart distribution for a full
-// p-dimensional graph with parameters delta and D : br p (reusable)
-Real gwish_nc_complete(myInt delta, int p, Real *D, bool flag) {
- Real d, c, a, g;
- Real dblDelta; // Recasting the inputs makes life easier below
- Real dblP;
- myInt i;
- c = 0.0;
- a = 0.0;
- g = 0.0;
- d = 0.0;
- dblDelta = delta;
- dblP = p;
-
- if (flag) {
- log_det(p, D, &d);
- }
- a = (dblDelta + dblP - 1) / 2.0;
- d = a * d;
- c = dblP * a * log_2;
- g = dblP * (dblP - 1) * log_pi_over_4;
-
- int signp;
-#ifdef ISFLOAT
- for (i = 0; i < p; i++)
- g += lgammaf_r(a - (Real)i / 2.0, &signp);
-#else
- for (i = 0; i < p; i++)
- g += lgamma_r(a - (Real)i / 2.0, &signp);
-#endif
-
- return (-d + c + g);
-}
-
-// Utility function for making submatrices
-void make_sub_mat_dbl(int p, int p_sub, myInt *sub, Real *A, Real *B) {
- int i, j;
- for (i = 0; i < p_sub; i++) {
- for (j = 0; j <= i; j++) {
- B[i * (i + 1) / 2 + j] =
- ((sub[i] >= sub[j]) ? A[sub[i] * (sub[i] + 1) / 2 + sub[j]]
- : A[sub[j] * (sub[j] + 1) / 2 + sub[i]]);
- }
- }
-}
-
-// Utility function for making a mean vector and a covariance matrix over a
-// subset of the dataset indicated by the vector sub.
-void make_sub_means_and_cov(Real *X, myInt *sub, myInt sub_match, int p,
- myInt n, myInt n_sub, Real *means, Real *D) {
- int i, j, k, ii;
- for (i = 0; i < p; i++) {
- means[i] = 0.0;
- };
- for (i = 0; i < p * (p + 1) / 2; i++) {
- D[i] = 0.0;
- }
-
- if (n_sub == 0) {
- return;
- }
- for (i = 0; i < p; i++) {
- for (k = 0; k < n; k++) {
- if (sub[k] == sub_match) {
- means[i] += (X[k * p + i] / (Real)n_sub);
- }
- }
- }
- if (n_sub < 2) {
- return;
- }
-
- for (i = 0; i < p; i++) {
- ii = i * (i + 1) / 2;
- for (j = 0; j <= i; j++) {
- for (k = 0; k < n; k++) {
- if (sub[k] == sub_match) {
- D[ii + j] += (X[k * p + i] - means[i]) * (X[k * p + j] - means[j]);
- }
- }
- }
- }
-
- return;
-}
-
-Real j_g_decomposable(LPGraph graph, Real *D_prior, Real *D_post, myInt delta,
- myInt n, bool flag) {
- Real mypost = 0;
- int p = graph->nVertices;
- myInt i;
- myInt sub_p;
- Real *sub_D = new Real[2 * p * p];
-
- //----- First loop through all the prime components (cliques since we're
- //decomposable) ---
- for (i = 0; i < graph->nCliques; i++) {
- sub_p = graph->CliquesDimens[i];
- if (flag) {
- make_sub_mat_dbl(p, sub_p, graph->Cliques[i], D_prior, sub_D);
- };
- mypost -= gwish_nc_complete(delta, sub_p, sub_D, flag);
- make_sub_mat_dbl(p, sub_p, graph->Cliques[i], D_post, sub_D);
- mypost += gwish_nc_complete(delta + n, sub_p, sub_D, 1);
- }
-
- //------- Now subtract off the separators -----------------------------------
- for (i = 0; i < graph->nSeparators; i++) {
- sub_p = graph->SeparatorsDimens[i];
- if (flag) {
- make_sub_mat_dbl(p, sub_p, graph->Separators[i], D_prior, sub_D);
- };
- mypost += gwish_nc_complete(delta, sub_p, sub_D, flag);
- make_sub_mat_dbl(p, sub_p, graph->Separators[i], D_post, sub_D);
- mypost -= gwish_nc_complete(delta + n, sub_p, sub_D, 1);
- }
-
- delete[] sub_D;
- return (mypost);
-}
\ No newline at end of file
diff --git a/src/sss-hip/kernels.cu b/src/sss-hip/kernels.cu
deleted file mode 100644
index e1218b4da0..0000000000
--- a/src/sss-hip/kernels.cu
+++ /dev/null
@@ -1,272 +0,0 @@
-#define KERNELS
-
-__global__ void CanDeleteEdge(myInt *d_in_delete, myInt *isDecomposable) {
- myInt tid = threadIdx.x;
- int bid = blockIdx.x;
- myInt bdim = blockDim.x;
-
- int n = *d_in_delete;
- myInt nCliques = *(d_in_delete + 1);
- myInt *CliquesDimens = d_in_delete + 2;
- myInt *Cliques = CliquesDimens + nCliques;
- // myInt nTasks = *(Cliques+nCliques*n);
- myInt *d_a = Cliques + nCliques * n + 1;
- myInt *d_b = d_a + *(Cliques + nCliques * n); // nTasks;
-
- myInt i, j, k;
- int ii; // myInt n = *d_n;
- myInt a = d_a[bid];
- myInt b = d_b[bid];
- __shared__ myInt count;
- if (tid == 0) {
- count = 0;
- }
- __shared__ myInt contain_a, contain_b, which_ab;
-
- for (i = 0; i < nCliques; i++) {
- ii = i * n;
- if (tid == 0) {
- contain_a = 0;
- contain_b = 0;
- }
- for (j = tid; j < CliquesDimens[i]; j += bdim) {
- k = Cliques[ii + j];
- if (k == a) {
- contain_a = 1;
- };
- if (k == b) {
- contain_b = 1;
- }
- }
- if (tid == 0) {
- if (contain_a && contain_b) {
- count++;
- which_ab = i;
- }
- }
- if (count > 1) {
- break;
- }
- }
-
- if (tid == 0) {
- if (count == 1) {
- isDecomposable[bid] = which_ab;
- } else {
- isDecomposable[bid] = -1;
- }
- }
-}
-
-// shared myInt demand: p+2*nTreeEdges+BLOCK_SIZE
-__global__ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add,
- myInt *isDecomposable) {
- myInt tid = threadIdx.x;
- int bid = blockIdx.x;
- myInt bdim = blockDim.x;
- extern __shared__ myInt shmem[];
- myInt *bi = shmem;
- int i, j, k, c, t;
-
- int n = (int)*d_in_delete;
- myInt nCliques = *(d_in_delete + 1);
- myInt *CliquesDimens = d_in_delete + 2;
- myInt *Cliques = CliquesDimens + nCliques;
-
- myInt *d_Labels = d_in_add;
- myInt nSeparators = *(d_in_add + n);
- myInt *SeparatorsDimens = d_in_add + n + 1;
- myInt *Separators = SeparatorsDimens + nSeparators;
- myInt nTreeEdges = *(Separators + n * nSeparators);
- myInt *TreeEdgeA = Separators + n * nSeparators + 1;
- myInt *TreeEdgeB = TreeEdgeA + nTreeEdges;
- myInt *d_Edge = TreeEdgeB + nTreeEdges;
- // myInt nTasks = *(d_Edge + n*n);
- myInt *d_a = d_Edge + (n * n) + 1;
- myInt *d_b = d_a + *(d_Edge + (n * n)); // nTasks;
-
- __shared__ myInt nS, pR, pT, contain_a, contain_b, aSi, bSi, common_parent, a,
- b;
- __shared__ myBool flag;
- myInt *R;
- myInt *T;
- myInt *S;
- myInt *contain_S;
- R = bi;
- bi += nTreeEdges;
- T = bi;
- bi += nTreeEdges;
- S = bi;
- bi += n;
- contain_S = bi; // bi += BLOCK_SIZE;
-
- if (tid == 0) {
- flag = 0;
- aSi = -1;
- bSi = -1;
- common_parent = 0;
- a = d_a[bid];
- b = d_b[bid];
- isDecomposable[bid] = 0;
-
- if (d_Labels[a] != d_Labels[b]) {
- flag = 1;
- isDecomposable[bid] = 1;
- }
- };
- SYNC;
-
- if (flag) {
- return;
- } // else ..
-
- if (tid == 0) {
- nS = 0;
- for (j = 0; j < n; j++) {
- if (d_Edge[a * n + j] && d_Edge[b * n + j]) {
- S[nS] = j;
- nS++;
- }
- };
- if (nS == 0) {
- flag = 1;
- }
- };
- SYNC;
-
- if (flag) {
- return;
- } // else ..
-
- // find higest-indexed cliques containing (a & S) and (b & S)
- for (i = 0; i < nCliques; i++) {
- if (tid == 0) {
- contain_a = 0;
- contain_b = 0;
- };
- contain_S[tid] = 0;
- t = i * n;
- for (j = tid; j < CliquesDimens[i]; j += bdim) {
- c = Cliques[t + j];
- if (c == a) {
- contain_a = 1;
- };
- if (c == b) {
- contain_b = 1;
- }
- for (k = 0; k < nS; k++) {
- if (c == S[k]) {
- contain_S[tid]++;
- }
- }
- }
- if (tid == 0) {
- k = 0;
- for (j = 0; j < BLOCK_SIZE; j++) {
- k += contain_S[j];
- }
- if (contain_a && (k == nS)) {
- aSi = i;
- };
- if (contain_b && (k == nS)) {
- bSi = i;
- }
- }
- }
-
- if (tid == 0) { // find the path from aSi to root
- R[0] = -1;
- pR = -1;
- for (i = 0; i < nTreeEdges; i++) {
- if (TreeEdgeA[i] == aSi) {
- R[0] = i;
- pR = 0;
- break;
- }
- }
- for (i = R[0]; i >= 0; i--) {
- if (TreeEdgeA[i] == TreeEdgeB[R[pR]]) {
- pR++;
- R[pR] = i;
- }
- }
- } else if (tid == 1) {
- // find the path from bSi to root
- T[0] = -1;
- pT = -1;
- for (i = 0; i < nTreeEdges; i++) {
- if (TreeEdgeA[i] == bSi) {
- T[0] = i;
- pT = 0;
- break;
- }
- }
- for (i = T[0]; i >= 0; i--) {
- if (TreeEdgeA[i] == TreeEdgeB[T[pT]]) {
- pT++;
- T[pT] = i;
- }
- }
- }
-
- if (tid == 0) {
- // find the branching point
- t = ((pR <= pT) ? pR : pT);
- for (i = 0; i <= t; i++) {
- if (TreeEdgeB[R[pR - i]] == TreeEdgeB[T[pT - i]]) {
- common_parent = i;
- } else {
- break;
- }
- }
- if (t != -1) {
- if (TreeEdgeA[R[pR - common_parent]] ==
- TreeEdgeA[T[pT - common_parent]]) {
- common_parent++;
- }
- }
- }
- SYNC;
-
- // check if S is in the path from R[0] to common_parent
- for (i = tid; i <= ((pR - common_parent) + (pT - common_parent) + 1);
- i += bdim) {
- if (i <= (pR - common_parent)) {
- if (SeparatorsDimens[R[i]] == nS) {
- contain_S[tid] = 0;
- t = R[i] * n;
- for (j = 0; j < nS; j++) {
- for (k = 0; k < nS; k++) {
- if (Separators[t + j] == S[k]) {
- contain_S[tid]++;
- }
- }
- }
- if (contain_S[tid] == nS) {
- flag = 1;
- isDecomposable[bid] = 1;
- }
- }
- } else {
- c = i - (pR - common_parent) - 1;
- if (SeparatorsDimens[T[c]] == nS) {
- contain_S[tid] = 0;
- t = T[c] * n;
- for (j = 0; j < nS; j++) {
- for (k = 0; k < nS; k++) {
- if (Separators[t + j] == S[k]) {
- contain_S[tid]++;
- }
- }
- }
- if (contain_S[tid] == nS) {
- flag = 1;
- isDecomposable[bid] = 1;
- }
- }
- }
- }
- SYNC;
-
- return;
-}
diff --git a/src/sss-hip/main.cu b/src/sss-hip/main.cu
index 014b5e8b28..a7ec485634 100644
--- a/src/sss-hip/main.cu
+++ b/src/sss-hip/main.cu
@@ -66,13 +66,11 @@
#define ISFLOAT 10
using namespace std;
-#include
-#include
+#include "gsl_compat.h"
#define GSL_INTEGRATION_GRIDSIZE 1000
gsl_integration_workspace *w;
gsl_function F;
-#include
#define RANDOMSEED 314159265
// Define hyperparameters for the prior distribution of (mu, K | G)
diff --git a/src/sss-hip/utilities.cpp b/src/sss-hip/utilities.cpp
deleted file mode 100644
index 00280ff271..0000000000
--- a/src/sss-hip/utilities.cpp
+++ /dev/null
@@ -1,74 +0,0 @@
-#define UTILITIES_CPP
-
-#ifndef GRAPH_CPP
-#include "graph.cpp"
-#endif
-
-// function to return a random myInteger with probability according to a given
-// (normalized!) weights
-myInt rand_myInt_weighted(myInt n, Real *weights) {
- myInt i;
- Real r = gsl_ran_flat(rnd, 0.0, 1.0);
- for (i = 0; i < n; i++) {
- if (r < weights[i]) {
- break;
- };
- r -= weights[i];
- };
- return i;
-}
-
-// returns a random myInterger between 0 and n - 1
-myInt rand_myInt(myInt n) {
- Real alpha = gsl_ran_flat(rnd, 0.0, 1.0);
- myInt value = (myInt)(n * alpha) / 1;
- return (value);
-}
-
-// function to return a random myInteger with probability according to a given
-// (normalized!) weights
-int rand_int_weighted(int n, Real *weights) {
- int i;
- Real r = gsl_ran_flat(rnd, 0.0, 1.0);
- for (i = 0; i < n; i++) {
- if (r < weights[i]) {
- break;
- };
- r -= weights[i];
- };
- return i;
-}
-
-// returns a random myInterger between 0 and n - 1
-int rand_int(int n) {
- Real alpha = gsl_ran_flat(rnd, 0.0, 1.0);
- int value = (int)(n * alpha) / 1;
- return (value);
-}
-
-// This samples an object from a subset of 0 to n - 1
-int sample_from(int *s, int n, int n_sub) {
- int i;
- int k = 0;
- int sample = rand_int(n_sub);
- for (i = 0; i < n; i++) {
- if (s[i]) {
- if (k == sample) {
- return (i);
- };
- k++;
- }
- }
-
- return (-1);
-}
-
-void shuffle(int *randomorder, int start, int end, int size) {
- int i, j, k;
- for (i = start; i < end; i++) {
- j = rand_int(size);
- k = randomorder[j];
- randomorder[j] = randomorder[i];
- randomorder[i] = k;
- }
-}
\ No newline at end of file
diff --git a/src/sss-sycl/CMakeLists.txt b/src/sss-sycl/CMakeLists.txt
index ce3b2ab3f2..7d13060bf0 100644
--- a/src/sss-sycl/CMakeLists.txt
+++ b/src/sss-sycl/CMakeLists.txt
@@ -4,6 +4,6 @@ add_hecbench_benchmark(
NAME sss
MODEL sycl
SOURCES main.cpp
+ INCLUDE_DIRS ${CMAKE_CURRENT_LIST_DIR}/../sss-cuda
CATEGORIES algorithms
- LINK_LIBRARIES gsl
)
diff --git a/src/sss-sycl/DPmixGGM.cpp b/src/sss-sycl/DPmixGGM.cpp
deleted file mode 100644
index 8dccb876e7..0000000000
--- a/src/sss-sycl/DPmixGGM.cpp
+++ /dev/null
@@ -1,516 +0,0 @@
-#include
-#define DPMIXGGM_CPP
-#ifndef GRAPH_CPP
-#include "graph.cpp"
-#endif
-#ifndef GWISH_CPP
-#include "gwish.cpp"
-#endif
-
-typedef class DPmixGGM *State;
-
-class DPmixGGM {
- //-------- Variables ---------------
-public:
- int n;
- int p;
- Real *X;
- myInt *xi;
- myInt L;
- LPGraph *graphlist;
- Real plp; // log-prior of the partition
- Real *pll; // array of log-likelihood for the cluster
- Real alpha;
-
- //--------- Functions --------------
-public:
- DPmixGGM(Real *data, myInt L_start, myInt n_obs, myInt p_model,
- Real edgeInclusionProb, ifstream &initfile);
- DPmixGGM(State a);
- ~DPmixGGM();
-
- void RandomStartAllXi(myInt L);
- void RandomStartAllG(myInt L, Real edgeInclusionProb);
-
- void RandomStart(myInt L, Real edgeInclusionProb);
- void InformedStart(ifstream &initfile);
-
- void ReadState(ifstream &initfile);
- void WriteState(ofstream &out, int itr);
- void CopyState(State a);
-
- Real partitionlogPrior(myInt theL, myInt *thexi, Real thealpha);
- Real lA(myInt the_n, myInt m, myInt k);
-
- Real cluster_k_loglikelihood(myInt k, myInt *thexi, LPGraph thegraph);
- Real predictiveDistribution(myInt k, myInt l, myInt *thexi, LPGraph thegraph);
-};
-
-//--------- This is the initializer ----------------------
-DPmixGGM::DPmixGGM(Real *data, myInt L_start, myInt n_obs, myInt p_model,
- Real edgeInclusionProb, ifstream &initfile) {
- myInt i;
- X = data;
- n = n_obs;
- p = p_model;
- L = L_start;
- xi = new myInt[n];
- graphlist = new LPGraph[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- };
- pll = new Real[L];
-
- // alpha = 20;
- alpha = 1;
-
-#ifdef RANDOMSTART
- RandomStartAllXi(L);
- RandomStartAllG(L, edgeInclusionProb);
-#else
- ReadState(initfile);
-#endif
-
- plp = partitionlogPrior(L, xi, alpha);
- for (i = 0; i < L; i++) {
- pll[i] = cluster_k_loglikelihood(i, xi, graphlist[i]);
- }
-}
-
-DPmixGGM::DPmixGGM(State a) {
- myInt i;
-
- n = a->n;
- p = a->p;
- L = a->L;
- X = a->X;
- alpha = a->alpha;
- plp = a->plp;
- xi = new myInt[n];
- for (i = 0; i < n; i++) {
- xi[i] = a->xi[i];
- }
- graphlist = new LPGraph[L];
- pll = new Real[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- graphlist[i]->CopyGraph(a->graphlist[i]);
- pll[i] = a->pll[i];
- }
-}
-
-//-------- Destructor ------------------------------
-DPmixGGM::~DPmixGGM() {
- myInt i;
- delete[] xi;
- for (i = 0; i < L; i++) {
- delete graphlist[i];
- };
- delete[] graphlist;
- delete[] pll;
-}
-
-void DPmixGGM::RandomStartAllXi(myInt L_start) {
- for (myInt i = 0; i < n; i++) {
- xi[i] = (myInt)(L_start * gsl_ran_flat(rnd, 0.0, 1.0)) / 1;
- }
-}
-
-void DPmixGGM::RandomStartAllG(myInt L_start, Real edgeInclusionProb) {
- myInt i, k, l;
- bool temp;
-
- for (i = 0; i < L_start; i++) {
- for (k = 0; k < p - 1; k++) {
- graphlist[i]->Edge[k][k] = 0;
- for (l = k + 1; l < p; l++) {
- temp = (gsl_ran_flat(rnd, 0.0, 1.0) < edgeInclusionProb);
- graphlist[i]->Edge[k][l] = temp;
- graphlist[i]->Edge[l][k] = temp;
- }
- }
- graphlist[i]->Edge[p - 1][p - 1] = 0;
- TurnFillInGraph(graphlist[i]);
- graphlist[i]->GenerateAllCliques();
- }
-}
-
-void DPmixGGM::RandomStart(myInt L_start, Real edgeInclusionProb) {
- myInt i;
- for (i = 0; i < L; i++) {
- delete graphlist[i];
- };
- delete graphlist;
- delete pll;
- L = L_start;
- graphlist = new LPGraph[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- };
- pll = new Real[L];
-
- RandomStartAllXi(L);
- RandomStartAllG(L, edgeInclusionProb);
-
- plp = partitionlogPrior(L, xi, alpha);
- for (i = 0; i < L; i++) {
- pll[i] = cluster_k_loglikelihood(i, xi, graphlist[i]);
- }
-}
-
-void DPmixGGM::InformedStart(ifstream &initfile) {
- ReadState(initfile);
- plp = partitionlogPrior(L, xi, alpha);
- for (myInt i = 0; i < L; i++) {
- pll[i] = cluster_k_loglikelihood(i, xi, graphlist[i]);
- }
-}
-
-void DPmixGGM::ReadState(ifstream &initfile) {
- myInt i, j, k, l;
-
- for (l = 0; l < L; l++) {
- delete graphlist[l];
- }
- delete[] graphlist;
- delete[] pll;
-
- int supern, superp;
- initfile >> supern;
- initfile >> superp;
- initfile >> L; // cout << "supern = " << supern << " superp = " << superp << "
- // L = " << L << endl; fflush(stdout);
- for (i = 0; i < n; i++) {
- initfile >> xi[i];
- xi[i]--;
- };
- for (i = n; i < supern; i++) {
- initfile >> j;
- }
- graphlist = new LPGraph[L];
- pll = new Real[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- for (k = 0; k < p; k++) {
- for (l = 0; l < p; l++) {
- initfile >> graphlist[i]->Edge[k][l];
- };
- for (l = p; l < superp; l++) {
- initfile >> j;
- }
- }
- for (k = p; k < superp; k++) {
- for (l = 0; l < superp; l++) {
- initfile >> j;
- }
- }
- TurnFillInGraph(graphlist[i]);
- graphlist[i]->GenerateAllCliques();
- }
-
- // cout << "end of ReadState" << endl; fflush(stdout);
-}
-
-void DPmixGGM::WriteState(ofstream &out, int itr) {
- myInt i, j, l, q, r;
- Real score = plp;
- for (l = 0; l < L; l++) {
- score += pll[l];
- }
-
- out << L << " " << score << " " << itr << " ";
- for (i = 0; i < n; i++) {
- out << xi[i] << " ";
- }
- for (l = 0; l < L; l++) {
- for (q = 0; q < p - 1; q++) {
- for (r = q + 1; r < p; r++) {
- out << graphlist[l]->Edge[q][r] << " ";
- }
- }
- }
- for (l = 0; l < L; l++) {
- out << graphlist[l]->nCliques << " ";
- for (i = 0; i < (graphlist[l]->nCliques); i++) {
- out << graphlist[l]->CliquesDimens[i] << " ";
- for (j = 0; j < (graphlist[l]->CliquesDimens[i]); j++) {
- out << graphlist[l]->Cliques[i][j] << " ";
- }
- }
-
- out << graphlist[l]->nTreeEdges << " ";
- for (i = 0; i < (graphlist[l]->nTreeEdges); i++) {
- out << graphlist[l]->TreeEdgeA[i] << " " << graphlist[l]->TreeEdgeB[i]
- << " ";
- }
-
- out << (graphlist[l]->nSeparators) << " ";
- for (i = 0; i < (graphlist[l]->nSeparators); i++) {
- out << graphlist[l]->SeparatorsDimens[i] << " ";
- for (j = 0; j < (graphlist[l]->SeparatorsDimens[i]); j++) {
- out << graphlist[l]->Separators[i][j] << " ";
- }
- }
- }
- out << endl;
-}
-
-void DPmixGGM::CopyState(State a) {
- myInt i;
- myInt oldL = L;
-
- n = a->n;
- p = a->p;
- X = a->X;
- alpha = a->alpha;
- plp = a->plp;
- for (i = 0; i < n; i++) {
- xi[i] = a->xi[i];
- };
- if (L != a->L) {
- L = a->L;
- pll = new Real[L];
- for (i = 0; i < oldL; i++) {
- delete graphlist[i];
- };
- delete[] graphlist;
- graphlist = new LPGraph[L];
- for (i = 0; i < L; i++) {
- graphlist[i] = new Graph();
- graphlist[i]->InitGraph(p);
- }
- }
- for (i = 0; i < L; i++) {
- graphlist[i]->CopyGraph(a->graphlist[i]);
- pll[i] = a->pll[i];
- }
-}
-
-#ifdef JEFFREYS_PRIOR
-struct f_params {
- int n;
- int m;
- int k;
-};
-double f(double beta, void *params) {
- struct f_params *iparams = (f_params *)params;
- int k = iparams->k;
- int n = iparams->n;
- int m = iparams->m;
- double s;
- double sum = 0;
- for (int j = 1; j < m; j++) {
- s = beta + j;
- sum += ((Real)j) / (s * s);
- }
- return exp(lgamma(beta) - lgamma(beta + n) + lgamma(n + 1.0) +
- (k - 0.5) * log(beta) + log(sqrt(sum)));
-}
-
-Real DPmixGGM::lA(myInt the_n, myInt m, myInt k) {
- struct f_params params = {the_n, m, k};
- F.function = &f;
- F.params = ¶ms;
- double result, error;
- gsl_integration_qagiu(&F, 0.0, 1e-7, 1e-7, GSL_INTEGRATION_GRIDSIZE, w,
- &result, &error);
- return log(result);
-}
-#endif
-
-//--------------------------------------------------------------------
-// This returns the partition prior probability -- normalising constant ignored
-// assuming alpha is fixed maxL is maximum number of clusters, effective number
-// of clusters can be smaller
-Real DPmixGGM::partitionlogPrior(myInt maxL, myInt *thexi, Real thealpha) {
- myInt i, j;
- Real siz;
- Real pri = 0;
- myInt effectiveL = 0;
- for (i = 0; i < maxL; i++) {
- siz = 0.0;
- for (j = 0; j < n; j++) {
- if (thexi[j] == i)
- siz = siz + 1;
- };
- if (siz > 0) {
- pri += lgamma(siz);
- effectiveL++;
- }
- }
-
-#ifdef JEFFREYS_PRIOR
- struct f_params params = {n, n, effectiveL};
- F.function = &f;
- F.params = ¶ms;
- double result, error;
- gsl_integration_qagiu(&F, 0.0, 1e-7, 1e-7, GSL_INTEGRATION_GRIDSIZE, w,
- &result, &error);
- return pri + log(result);
-#else
- return pri + effectiveL * log(thealpha);
-#endif
-}
-
-// This returns the likelihood of the observations in cluster k
-Real DPmixGGM::cluster_k_loglikelihood(myInt k, myInt *thexi,
- LPGraph thegraph) {
- //-------- Variables, follows paper ---------------
- int i, j, ii;
- Real n0 = N0;
- Real *mu0 = new Real[p];
- for (i = 0; i < p; i++)
- mu0[i] = 0;
- Real *xbar = new Real[p];
- Real *mu_bar = new Real[p];
- Real *D_prior = new Real[p * (p + 1) / 2];
- Real *D_post = new Real[p * (p + 1) / 2];
- myInt n_sub = 0;
- Real J_G;
- Real Norm_terms;
-
- for (i = 0; i < n; i++) {
- if (thexi[i] == k)
- n_sub++;
- } // count the number in cluster k
- if (n_sub == 0) {
- delete[] D_post;
- delete[] D_prior;
- delete[] xbar;
- delete[] mu_bar;
- delete[] mu0;
- return (0);
- }
- // cout << "l = " << k << ", n_sub = " << n_sub << endl; fflush(stdout);
-
- //----------- Form the sufficient statistics ----------
- make_sub_means_and_cov(X, thexi, k, p, n, n_sub, xbar, D_post);
-
- //------------------ Update these ---------------------
- for (i = 0; i < p; i++)
- mu_bar[i] = (n_sub * xbar[i] + n0 * mu0[i]) / (n_sub + n0);
- for (i = 0; i < p * (p + 1) / 2; i++)
- D_prior[i] = 0;
- for (i = 0; i < p; i++)
- D_prior[i * (i + 1) / 2 + i] = 1;
- for (i = 0; i < p; i++)
- D_post[i * (i + 1) / 2 + i] += 1;
-
- //----------------- Factor in the mean for the matrix D_posterior
- //-------------
- for (i = 0; i < p; i++) {
- ii = i * (i + 1) / 2;
- for (j = 0; j <= i; j++) {
- D_post[ii + j] += (-(n_sub + n0) * mu_bar[i] * mu_bar[j] +
- n_sub * xbar[i] * xbar[j] + n0 * mu0[i] * mu0[j]);
- }
- }
-
- //------------ Calculate the score
- //---------------------------------------------
- J_G = j_g_decomposable(thegraph, D_prior, D_post, DELTA0, n_sub, 0);
- Norm_terms =
- -(Real(n_sub * p) / 2) * log_2_pi + Real(p) / 2 * log(n0 / (n_sub + n0));
-
- delete[] D_post;
- delete[] D_prior;
- delete[] xbar;
- delete[] mu_bar;
- delete[] mu0;
- return (Norm_terms + J_G);
-}
-
-//--------------------------------------------------------------------
-// This returns the predictive distribution score for observation k, if it were
-// to belong to cluster l thegraph is the graph associated with cluster l
-Real DPmixGGM::predictiveDistribution(myInt k, myInt l, myInt *thexi,
- LPGraph thegraph) {
- int i, j, ii;
- Real a;
- Real n0 = N0;
- Real *mu0 = new Real[p];
- Real *xbar = new Real[p];
- Real *mu_bar = new Real[p];
- Real *mu_tilde = new Real[p];
- for (i = 0; i < p; i++) {
- mu0[i] = 0;
- }
- Real *D_prior = new Real[p * (p + 1) / 2];
- Real *D_post = new Real[p * (p + 1) / 2];
- myInt n_sub = 0;
- Real coef;
- Real J_G;
- Real Norm_terms;
-
- for (i = 0; i < n; i++) {
- if (thexi[i] == l) {
- n_sub++;
- }
- } // count the number in this cluster
-#ifdef JEFFREYS_PRIOR
- if (l == L) {
- coef = lA(n, n, L + 1) - lA(n - 1, n, L);
- } else if (n_sub > 0) {
- coef = log(n_sub) + lA(n, n, L) - lA(n - 1, n, L);
- } else {
- coef = 0;
- }
-#else
- if (l == L) {
- coef = log(alpha);
- } else if (n_sub > 0) {
- coef = log(n_sub);
- } else {
- coef = 0;
- }
-#endif
-
- //----------- Form the sufficient statistics ----------
- make_sub_means_and_cov(X, thexi, l, p, n, n_sub, xbar, D_prior);
-
- //------------------ Update these ---------------------
- for (i = 0; i < p; i++)
- mu_bar[i] = (n_sub * xbar[i] + n0 * mu0[i]) / (n_sub + n0);
- for (i = 0; i < p; i++)
- D_prior[i * (i + 1) / 2 + i] += 1;
- //----------------- Factor in the mean for the matrix D_prior -------------
- for (i = 0; i < p; i++) {
- ii = i * (i + 1) / 2;
- for (j = 0; j <= i; j++) {
- D_prior[ii + j] += -(n_sub + n0) * mu_bar[i] * mu_bar[j] +
- n_sub * xbar[i] * xbar[j] + n0 * mu0[i] * mu0[j];
- }
- }
-
- //-------------- Get the posterior information ----------------------------
- for (i = 0; i < p * (p + 1) / 2; i++)
- D_post[i] = D_prior[i];
- for (i = 0; i < p; i++)
- mu_tilde[i] = (X[k * p + i] + (n_sub + n0) * mu_bar[i]) / (n_sub + n0 + 1);
- for (i = 0; i < p; i++) {
- ii = i * (i + 1) / 2;
- for (j = 0; j <= i; j++) {
- a = -(n_sub + 1 + n0) * mu_tilde[i] * mu_tilde[j] +
- X[k * p + i] * X[k * p + j] + (n_sub + n0) * mu_bar[i] * mu_bar[j];
- D_post[ii + j] += a;
- }
- }
-
- //------------ Calculate the score
- //---------------------------------------------
- J_G = j_g_decomposable(thegraph, D_prior, D_post, DELTA0 + n_sub, 1, 1);
- Norm_terms =
- -(p / 2) * log_2_pi + (p / 2) * log((n_sub + n0) / (n_sub + 1 + n0));
-
- delete[] D_post;
- delete[] D_prior;
- delete[] xbar;
- delete[] mu_bar;
- delete[] mu_tilde;
- delete[] mu0;
- return (coef + Norm_terms + J_G);
-}
\ No newline at end of file
diff --git a/src/sss-sycl/DPmixGGM_Lists.cpp b/src/sss-sycl/DPmixGGM_Lists.cpp
deleted file mode 100644
index 191ced7e5e..0000000000
--- a/src/sss-sycl/DPmixGGM_Lists.cpp
+++ /dev/null
@@ -1,413 +0,0 @@
-#define LISTS_CPP
-#ifndef GRAPH_CPP
-#include "graph.cpp"
-#endif
-#ifndef GWISH_CPP
-#include "gwish.cpp"
-#endif
-#ifndef DPMIXGGM_CPP
-#include "DPmixGGM.cpp"
-#endif
-
-typedef class DPmixGGMlist *List;
-
-class DPmixGGMlist {
- // variables
-public:
- int n;
- int p;
- myInt M;
-
- myInt *L_list;
- myInt *xi_list;
- myInt *edge_list;
- Real *score_list;
-
- myInt ind_minOfListScores;
- Real minOfListScores;
-
- // functions
-public:
- DPmixGGMlist(myInt size, int n, int p);
- // DPmixGGMlist (myInt size, State a);
- ~DPmixGGMlist();
-
- void FlushList(State a);
- void WriteList(ofstream &out);
- void ReadFromList(State a, myInt m);
-
- void UpdateList(myInt L, myInt *xi, LPGraph *graphlist, Real score);
- void UpdateList(State a);
-
- LPGraph ProposeGraph(myInt NN_xi, myInt *which_xi, Real sFactor);
-};
-
-DPmixGGMlist::DPmixGGMlist(myInt size, int in_n, int in_p) {
- int i;
- n = in_n;
- p = in_p;
- M = size;
- int ee = p * (p - 1) / 2;
-
- L_list = new myInt[M];
- for (i = 0; i < M; i++) {
- L_list[i] = -1;
- }
- xi_list = new myInt[M * n];
- for (i = 0; i < M * n; i++) {
- xi_list[i] = -1;
- }
- edge_list = new myInt[M * n * ee];
- for (i = 0; i < (M * n * ee); i++) {
- edge_list[i] = -1;
- }
- score_list = new Real[M];
- for (i = 0; i < M; i++) {
- score_list[i] = NEG_INF;
- }
-
- ind_minOfListScores = 0;
- minOfListScores = NEG_INF;
-}
-
-DPmixGGMlist::~DPmixGGMlist() {
- delete[] L_list;
- delete[] xi_list;
- delete[] edge_list;
- delete[] score_list;
-}
-
-void DPmixGGMlist::FlushList(State a) {
- int i;
- int ee = p * (p - 1) / 2;
-
- for (i = 0; i < M; i++) {
- L_list[i] = -1;
- }
- for (i = 0; i < M * n; i++) {
- xi_list[i] = -1;
- }
- for (i = 0; i < (M * n * ee); i++) {
- edge_list[i] = -1;
- }
- for (i = 0; i < M; i++) {
- score_list[i] = NEG_INF;
- }
-
- int q, r, l;
- L_list[0] = a->L;
- for (i = 0; i < n; i++) {
- xi_list[i] = a->xi[i];
- }
- for (i = 0; i < n; i++) {
- l = 0;
- for (q = 0; q < p - 1; q++) {
- for (r = q + 1; r < p; r++) {
- edge_list[i * ee + l] = a->graphlist[a->xi[i]]->Edge[q][r];
- l++;
- }
- }
- }
- score_list[0] = a->plp;
- for (i = 0; i < (a->L); i++) {
- score_list[0] += a->pll[i];
- }
-
- ind_minOfListScores = 1;
- minOfListScores = NEG_INF;
-}
-
-void DPmixGGMlist::WriteList(ofstream &out) {
- int i, j, l, L, q, r, t;
- int ee = p * (p - 1) / 2;
-
- for (t = 0; t < M; t++) {
- if (L_list[t] == -1) {
- continue;
- }
-
- L = L_list[t];
- out << L << " " << score_list[t] << " ";
- for (i = 0; i < n; i++) {
- out << xi_list[t * n + i] << " ";
- }
-
- LPGraph *graphlist = new LPGraph[L];
- for (l = 0; l < L; l++) {
- graphlist[l] = new Graph();
- graphlist[l]->InitGraph(p);
- }
- for (l = 0; l < L; l++) {
- for (i = 0; i < n; i++) {
- if (xi_list[t * n + i] == l)
- break;
- }
- j = 0;
- for (q = 0; q < p - 1; q++) {
- graphlist[l]->Edge[q][q] = 0;
- for (r = q + 1; r < p; r++) {
- graphlist[l]->Edge[q][r] = edge_list[t * n * ee + i * ee + j];
- graphlist[l]->Edge[r][q] = graphlist[l]->Edge[q][r];
- out << graphlist[l]->Edge[q][r] << " ";
- j++;
- }
- }
- graphlist[l]->Edge[p - 1][p - 1] = 0;
- graphlist[l]->GenerateAllCliques();
- }
-
- for (l = 0; l < L; l++) {
- out << graphlist[l]->nCliques << " ";
- for (i = 0; i < (graphlist[l]->nCliques); i++) {
- out << graphlist[l]->CliquesDimens[i] << " ";
- for (j = 0; j < (graphlist[l]->CliquesDimens[i]); j++) {
- out << graphlist[l]->Cliques[i][j] << " ";
- }
- }
-
- out << graphlist[l]->nTreeEdges << " ";
- for (i = 0; i < (graphlist[l]->nTreeEdges); i++) {
- out << graphlist[l]->TreeEdgeA[i] << " " << graphlist[l]->TreeEdgeB[i]
- << " ";
- }
-
- out << graphlist[l]->nSeparators << " ";
- for (i = 0; i < (graphlist[l]->nSeparators); i++) {
- out << graphlist[l]->SeparatorsDimens[i] << " ";
- for (j = 0; j < (graphlist[l]->SeparatorsDimens[i]); j++) {
- out << graphlist[l]->Separators[i][j] << " ";
- }
- }
- }
-
- out << endl;
- for (l = 0; l < L; l++) {
- delete graphlist[l];
- };
- delete[] graphlist;
- }
-}
-
-void DPmixGGMlist::ReadFromList(State a, myInt m) {
- int i, j, k, l, L, t;
- int ee = p * (p - 1) / 2;
-
- for (l = 0; l < (a->L); l++) {
- delete a->graphlist[l];
- };
- delete[] a->graphlist;
- delete[] a->pll; // delete previous memory spaces
- L = L_list[m];
- a->L = L;
- a->pll = new Real[L];
- a->graphlist = new LPGraph[L];
- for (l = 0; l < L; l++) {
- a->graphlist[l] = new Graph();
- a->graphlist[l]->InitGraph(p);
- }
- for (i = 0; i < n; i++) {
- a->xi[i] = xi_list[m * n + i];
- }
- for (l = 0; l < L; l++) {
- for (k = 0; k < n; k++) {
- if (a->xi[k] == l) {
- break;
- }
- };
- t = 0;
- for (i = 0; i < (p - 1); i++) {
- for (j = (i + 1); j < p; j++) {
- a->graphlist[l]->Edge[i][j] = edge_list[m * n * ee + k * ee + t];
- a->graphlist[l]->Edge[j][i] = a->graphlist[l]->Edge[i][j];
- t++;
- }
- }
- a->graphlist[l]->GenerateAllCliques();
- }
-}
-
-// OVERWRITES A GIVEN LIST OF MODELS WITH A GIVEN MODEL IF IT HAS HIGHER SCORE
-void DPmixGGMlist::UpdateList(myInt L, myInt *xi, LPGraph *graphlist,
- Real score) {
- int i, j, k, l, q, r, t;
- int ee = p * (p - 1) / 2;
-
- // check if the candidate model deserve to be in the list of models -- if so,
- // do the necessary
- bool flag;
- int count;
- myInt *xi_new = new myInt[n];
- myInt *index_set = new myInt[L];
- for (i = 0; i < n; i++) {
- xi_new[i] = -1;
- }
-
- if (score > minOfListScores) {
- // reindexing (lexicographically) the incoming model -- if there is a
- // redundant index, getting rid of it
- myInt xi_temp = xi[0];
- index_set[0] = xi_temp;
- for (i = 0; i < n; i++) {
- if (xi[i] == xi_temp) {
- xi_new[i] = 0;
- }
- }
-
- j = 1;
- l = 0;
- while (j < n) {
- k = j;
- l++;
- while (k < n) // finding a new cluster index
- {
- flag = 1;
- for (t = 0; t < l; t++) {
- if (xi[k] == index_set[t]) {
- flag = 0;
- break;
- }
- }
- if (flag) {
- xi_temp = xi[k];
- index_set[l] = xi_temp;
- for (i = 0; i < n; i++) {
- if (xi[i] == xi_temp) {
- xi_new[i] = l;
- }
- }
- break;
- } else {
- k++;
- }
- }
- j = k + 1;
- }
- k = l;
-
- // check if this model is already in the list of models
- count = 0;
- for (t = 0; t < M; t++) {
- flag = 0;
- if (L_list[t] != k) {
- flag = 1;
- }
-
- if (flag == 0) {
- for (i = 0; i < n; i++) {
- if (xi_list[t * n + i] != xi_new[i]) {
- flag = 1;
- break;
- }
- }
- };
-
- if (flag == 0) {
- for (l = 0; l < k; l++) {
- for (i = 0; i < n; i++) {
- if (xi_list[t * n + i] == l)
- break;
- } // finiding an observation with cluster index l
-
- j = 0;
- for (q = 0; q < p - 1; q++) {
- for (r = q + 1; r < p; r++) {
- if (edge_list[t * n * ee + i * ee + j] !=
- graphlist[xi[i]]->Edge[q][r]) {
- flag = 1;
- break;
- };
- j++;
- }
- if (flag)
- break;
- }
- if (flag)
- break;
- }
- }
-
- if (flag)
- count++;
- }
-
- // if the model is not in list, substitute the lowest scoring model with the
- // candidate
- if (count == M) {
- L_list[ind_minOfListScores] = k; // cout << "effective L = " << k << endl;
- for (i = 0; i < n; i++) {
- xi_list[ind_minOfListScores * n + i] = xi_new[i];
- t = 0;
- for (q = 0; q < p - 1; q++) {
- for (r = q + 1; r < p; r++) {
- edge_list[ind_minOfListScores * n * ee + i * ee + t] =
- graphlist[xi[i]]->Edge[q][r];
- t++;
- }
- }
- }
- score_list[ind_minOfListScores] = score;
-
- // searching for minimum scoring model in the list
- ind_minOfListScores = 0;
- minOfListScores = score_list[0];
- for (t = 1; t < M; t++) {
- if (score_list[t] < minOfListScores) {
- minOfListScores = score_list[t];
- ind_minOfListScores = t;
- }
- }
- }
- }
-
- delete[] xi_new;
- delete[] index_set;
-}
-
-void DPmixGGMlist::UpdateList(State a) {
- Real score = a->plp;
- for (myInt i = 0; i < a->L; i++) {
- score += a->pll[i];
- }
- UpdateList(a->L, a->xi, a->graphlist, score);
-}
-
-// PROPOSES A DECOMPOSABLE GRAPH
-LPGraph DPmixGGMlist::ProposeGraph(myInt NN_xi, myInt *which_xi, Real sFactor) {
- int i, j, k, l, m, t;
- Real a, s;
- bool temp;
- int ee = p * (p - 1) / 2;
-
- LPGraph newgraph = new Graph();
- newgraph->InitGraph(p);
-
- // feature extraction
- l = 0;
- for (i = 0; i < p - 1; i++) {
- newgraph->Edge[i][i] = 0;
- for (j = i + 1; j < p; j++) {
- s = NN_xi * M * sFactor;
- a = s;
-
- for (m = 0; m < M; m++) {
- for (k = 0; k < NN_xi; k++) {
- t = edge_list[m * n * ee + which_xi[k] * ee + l];
- if (t != -1) {
- a += Real((bool)t);
- s += 1.0;
- }
- }
- }
- l++;
-
- temp = (gsl_ran_flat(rnd, 0.0, 1.0) < a / s);
- newgraph->Edge[i][j] = temp;
- newgraph->Edge[j][i] = temp;
- }
- }
- newgraph->Edge[p - 1][p - 1] = 0;
- TurnFillInGraph(newgraph);
- newgraph->GenerateAllCliques();
-
- return newgraph;
-}
\ No newline at end of file
diff --git a/src/sss-sycl/Makefile b/src/sss-sycl/Makefile
index 97e0da3d8a..886dba6c8b 100644
--- a/src/sss-sycl/Makefile
+++ b/src/sss-sycl/Makefile
@@ -6,7 +6,6 @@
CC = clang++
OPTIMIZE = yes
DEBUG = no
-GSL = /path/to/GSL
LAUNCHER =
CUDA = no
@@ -30,8 +29,7 @@ obj = $(source:.cpp=.o)
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -fsycl \
- -I$(GSL)/include \
+CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -fsycl -I../sss-cuda \
--gcc-toolchain=$(GCC_TOOLCHAIN)
ifeq ($(VENDOR), AdaptiveCpp)
@@ -40,7 +38,7 @@ ifeq ($(VENDOR), AdaptiveCpp)
endif
# Linker Flags
-LDFLAGS = -L$(GSL)/lib -lgsl
+LDFLAGS =
ifeq ($(CUDA), yes)
CFLAGS += -fsycl-targets=nvptx64-nvidia-cuda \
@@ -70,12 +68,16 @@ endif
$(program): $(obj)
$(CC) $(CFLAGS) $(obj) -o $@ $(LDFLAGS)
-%.o: %.cpp DPmixGGM.cpp DPmixGGM_Lists.cpp DPmixGGM_SSSmoves.cpp \
- kernels.cpp graph.cpp gwish.cpp utilities.cpp
+SHARED_SRC = ../sss-cuda/DPmixGGM.cpp ../sss-cuda/DPmixGGM_Lists.cpp \
+ ../sss-cuda/graph.cpp ../sss-cuda/gwish.cpp \
+ ../sss-cuda/utilities.cpp ../sss-cuda/graph.h \
+ ../sss-cuda/gsl_compat.h
+
+%.o: %.cpp DPmixGGM_SSSmoves.cpp kernels.cpp $(SHARED_SRC)
$(CC) $(CFLAGS) -c $< -o $@
clean:
rm -rf $(program) $(obj) RES/f9_n150_p50_modes_GPU.txt
run: $(program)
- $(LAUNCHER) GSL_RNG_SEED=123 ./$(program) f9_n150_p50
+ $(LAUNCHER) ./$(program) f9_n150_p50
diff --git a/src/sss-sycl/README.md b/src/sss-sycl/README.md
index d020e7e2f9..d313bd8d90 100644
--- a/src/sss-sycl/README.md
+++ b/src/sss-sycl/README.md
@@ -1,22 +1,28 @@
DPmixGGM
========
-This folder contains source codes for the "GPU-powered Stochastic Shotgun Search for Dirichlet proces mixtures of Gaussian Graphical Models"
+This folder contains source codes for the "GPU-powered Stochastic Shotgun Search for Dirichlet proces mixtures of Gaussian Graphical Models"
by Chiranjit Mukherjee and Abel Rodriguez.
-
-The "DPmixGGM_SSS_main.cpp" file contains tuning parameters for the algorithm, as elaborated below:
+
+The "DPmixGGM_SSS_main.cpp" file contains tuning parameters for the algorithm, as elaborated below:
1. Run the SSS
-2. Run GPU/CPU versions of the SSS by enabling / disabling the macro CUDA.
-3. Specify maximum number of mixture components that the model should accommodate (for pre-allocation of memory).
-4. Set SSS runtime parameters C, D, R, S, M, g, h, f, t.
-5. Set SSS number of chain parameters. User needs to provide at least one initial point.
-7. Set hyperparameters of for the prior on (mu, K | G) with N0, DELTA0.
-
+2. Run GPU/CPU versions of the SSS by enabling / disabling the macro CUDA.
+3. Specify maximum number of mixture components that the model should accommodate (for pre-allocation of memory).
+4. Set SSS runtime parameters C, D, R, S, M, g, h, f, t.
+5. Set SSS number of chain parameters. User needs to provide at least one initial point.
+7. Set hyperparameters of for the prior on (mu, K | G) with N0, DELTA0.
+
Complie source codes using the "make" command and run with "./main f9_n150_p50" command.
-
-The program expects an input-data file (e.g. f9_n150_p50) in the DATA/ folder and at least one initialization point (e.g. f9_n150_p50_init1).
-The input-data file should specify n and p in the first row and then provide n rows of length p. The initial point data-file should specify n, p
-and L of the initial model configuration in the first row and xi-indices of the initial point in the second row. Subsequent L rows specify
-G_l (l=1:L).
-
+
+The program expects an input-data file (e.g. f9_n150_p50) in the DATA/ folder and at least one initialization point (e.g. f9_n150_p50_init1).
+The input-data file should specify n and p in the first row and then provide n rows of length p. The initial point data-file should specify n, p
+and L of the initial model configuration in the first row and xi-indices of the initial point in the second row. Subsequent L rows specify
+G_l (l=1:L).
+
A list of highest-score models is stored in folder RES/.
+
+Dependencies
+------------
+None beyond a C++17 compiler. The uniform random number generator (MT19937) and
+the adaptive Gauss-Kronrod quadrature the sampler needs are implemented in
+`gsl_compat.h` in this folder.
diff --git a/src/sss-sycl/graph.cpp b/src/sss-sycl/graph.cpp
deleted file mode 100644
index 55e75deb02..0000000000
--- a/src/sss-sycl/graph.cpp
+++ /dev/null
@@ -1,973 +0,0 @@
-#include
-#include
-#include
-#include
-
-#define GRAPH_CPP
-#ifndef GRAPH_H
-#include "graph.h"
-#endif
-#ifndef GWISH_CPP
-#include "gwish.cpp"
-#endif
-
-// class Graph::Begins
-Graph::Graph() {
- nVertices = 0;
- Edge = NULL;
- Labels = NULL;
- nLabels = 0;
- Cliques = NULL;
- CliquesDimens = NULL;
- nCliques = 0;
- TreeEdgeA = NULL;
- TreeEdgeB = NULL;
- nTreeEdges = 0;
- Separators = NULL;
- SeparatorsDimens = NULL;
- nSeparators = 0;
- localord = NULL;
- return;
-}
-
-Graph::Graph(LPGraph InitialGraph) {
- nVertices = 0;
- Edge = NULL;
- Labels = NULL;
- nLabels = 0;
- Cliques = NULL;
- CliquesDimens = NULL;
- nCliques = 0;
- TreeEdgeA = NULL;
- TreeEdgeB = NULL;
- nTreeEdges = 0;
- Separators = NULL;
- SeparatorsDimens = NULL;
- nSeparators = 0;
- localord = NULL;
-
- ///////////////////////////////////////
- myInt i, j;
- InitGraph(InitialGraph->nVertices);
- for (i = 0; i < nVertices; i++) {
- for (j = 0; j < nVertices; j++) {
- Edge[i][j] = InitialGraph->Edge[i][j];
- }
- }
-
- return;
-}
-
-Graph::~Graph() {
- myInt i; // cout << "-> "; fflush(stdout);
-
- for (i = 0; i < nVertices; i++) {
- delete[] Edge[i];
- Edge[i] = NULL;
- }
- delete[] Edge;
- Edge = NULL;
-
- delete[] Labels;
- Labels = NULL;
-
- delete[] Cliques[0]; // for(i=0; inVertices;
- for (i = 0; i < n; i++) {
- for (j = 0; j < n; j++) {
- Edge[i][j] = G->Edge[i][j];
- }
- }
- nLabels = G->nLabels;
- for (i = 0; i < n; i++) {
- Labels[i] = G->Labels[i];
- }
- nCliques = G->nCliques;
- for (i = 0; i < n; i++) {
- CliquesDimens[i] = G->CliquesDimens[i];
- for (j = 0; j < n; j++) {
- Cliques[i][j] = G->Cliques[i][j];
- }
- }
- nTreeEdges = G->nTreeEdges;
- for (i = 0; i < n; i++) {
- TreeEdgeA[i] = G->TreeEdgeA[i];
- TreeEdgeB[i] = G->TreeEdgeB[i];
- }
- nSeparators = G->nSeparators;
- for (i = 0; i < n; i++) {
- SeparatorsDimens[i] = G->SeparatorsDimens[i];
- for (j = 0; j < n; j++) {
- Separators[i][j] = G->Separators[i][j];
- }
- }
- for (i = 0; i < n; i++) {
- localord[i] = G->localord[i];
- }
-}
-
-void Graph::GenerateCliques(myInt label) {
- myInt i, j, k, p, r;
- myInt n = nVertices;
- myInt *clique = new myInt[nVertices];
- memset(clique, 0, nVertices * sizeof(myInt));
-
- myInt countA, countB;
- myInt *listA = new myInt[n];
- myInt *listB2 = new myInt[n];
-
- // clean memory
- memset(localord, 0, nVertices * sizeof(myInt));
-
- myInt v, vk;
- myInt PrevCard = 0;
- myInt NewCard;
- myInt s = nCliques - 1; // cout << "s = " << s << endl;
-
- countA = 0;
- for (i = 0; i < n; i++) {
- if (Labels[i] == label) {
- listA[countA] = i;
- countA++;
- }
- };
- countB = 0;
-
- for (i = n; i >= 0; i--) {
- NewCard = -1;
-
- // choose a vertex v...
- for (j = 0; j < countA; j++) {
- myInt maxj = 0;
- for (r = 0; r < countB; r++) {
- if (Edge[listA[j]][listB2[r]]) {
- maxj++;
- }
- }
- if (maxj > NewCard) {
- v = listA[j];
- NewCard = maxj;
- }
- }
-
- // printf("i=%d, NewCard=%d, PrevCard=%d countA=%d,
- // countB=%d\n",i,NewCard,PrevCard,countA,countB);
-
- if (NewCard == -1) {
- break;
- }
-
- localord[v] = i;
- if (NewCard <= PrevCard) { // begin new clique
- s++;
-
- for (r = 0; r < countB; r++) {
- if (Edge[v][listB2[r]]) {
- Cliques[s][CliquesDimens[s]] = listB2[r];
- CliquesDimens[s]++;
- }
- }
-
- if (NewCard != 0) { // get edge to parent
- vk = Cliques[s][0];
- k = localord[vk]; // cout << "(" << clique[Cliques[s][0]] << " ";
- for (r = 1; r < CliquesDimens[s]; r++) {
- if (localord[Cliques[s][r]] < k) {
- vk = Cliques[s][r];
- k = localord[vk];
- }; // cout << clique[Cliques[s][r]] << " ";
- }
- // cout << "| ";
- // for(r=0; rIsDecomposable())
- return;
-
- myInt v1 = gfill->SearchVertex();
- // add edges to Def(Adj(x)) so that Adj(x) becomes a clique
- for (u = 0; u < gfill->nVertices; u++) {
- if (gfill->Edge[v1][u] == 1) {
- for (v = u + 1; v < gfill->nVertices; v++) {
- if ((gfill->Edge[v1][v] == 1) && (gfill->Edge[u][v] == 0)) {
- gfill->Edge[v][u] = gfill->Edge[u][v] = 1;
- }
- }
- }
- }
- LPEliminationGraph egraph = new EliminationGraph(graph, v1);
- for (i = 1; i < graph->nVertices - 1; i++) {
- v1 = egraph->SearchVertex();
- for (u = 0; u < egraph->nVertices; u++) {
- if (egraph->Eliminated[u])
- continue;
- if (egraph->Edge[v1][u] == 1) {
- for (v = u + 1; v < egraph->nVertices; v++) {
- if (egraph->Eliminated[v])
- continue;
- if ((egraph->Edge[v1][v] == 1) && (egraph->Edge[u][v] == 0)) {
- gfill->Edge[v][u] = gfill->Edge[u][v] = 1;
- // these are the edges that are added to the initial graph
- }
- }
- }
- }
- egraph->EliminateVertex(v1);
- }
- delete egraph;
- return;
-}
-
-// -----------------------------------------------------------------------------------------------
-// Based on Scott & Carvalho 2008
-
-bool Graph::CanAddEdge(myInt a, myInt b) {
- myInt i, j, k;
- bool canadd = 0;
- int n = (int)nVertices;
- myInt nS = 0, pR, pT;
- myInt *R = new myInt[nTreeEdges];
- myInt *T = new myInt[nTreeEdges];
- myInt *S = new myInt[n];
- myInt common_parent = 0;
- myInt contain_a, contain_b, contain_S;
-
- if (Labels[a] != Labels[b]) {
- canadd = 1;
- } // else ..
- else {
-
- nS = 0;
- for (j = 0; j < n; j++) {
- if (Edge[a][j] && Edge[b][j]) {
- S[nS] = j;
- nS++;
- }
- }
- if (nS == 0) {
- canadd = 0;
- } // else ....
- else { // HERE;
-
- // find higest-indexed cliques containing (a & S) and (b & S)
- myInt aSi = -1, bSi = -1;
- for (i = 0; i < nCliques; i++) {
- contain_a = 0;
- contain_b = 0;
- contain_S = 0;
- for (j = 0; j < CliquesDimens[i]; j++) {
- if (Cliques[i][j] == a) {
- contain_a = 1;
- };
- if (Cliques[i][j] == b) {
- contain_b = 1;
- }
- for (k = 0; k < nS; k++) {
- if (Cliques[i][j] == S[k]) {
- contain_S++;
- }
- }
- }
- if (contain_a && (contain_S == nS)) {
- aSi = i;
- }
- if (contain_b && (contain_S == nS)) {
- bSi = i;
- }
- }
-
- // find the path from aSi to root
- R[0] = -1;
- pR = -1;
- for (i = 0; i < nTreeEdges; i++) {
- if (TreeEdgeA[i] == aSi) {
- R[0] = i;
- pR = 0;
- break;
- }
- }
- for (i = R[0]; i >= 0; i--) {
- if (TreeEdgeA[i] == TreeEdgeB[R[pR]]) {
- pR++;
- R[pR] = i;
- }
- }
-
- // find the path from bSi to root
- T[0] = -1;
- pT = -1;
- for (i = 0; i < nTreeEdges; i++) {
- if (TreeEdgeA[i] == bSi) {
- T[0] = i;
- pT = 0;
- break;
- }
- }
- for (i = T[0]; i >= 0; i--) {
- if (TreeEdgeA[i] == TreeEdgeB[T[pT]]) {
- pT++;
- T[pT] = i;
- }
- }
-
- // find the branching point
- k = (pR < pT) ? pR : pT; // min(pR,pT);
- for (i = 0; i <= k; i++) {
- if (TreeEdgeB[R[pR - i]] == TreeEdgeB[T[pT - i]]) {
- common_parent = i;
- } else {
- break;
- }
- }
- if (k != -1) {
- if (TreeEdgeA[R[pR - common_parent]] ==
- TreeEdgeA[T[pT - common_parent]]) {
- common_parent++;
- }
- }
-
- // check if S is in the path from R[0] to common_parent
- for (i = 0; i <= (pR - common_parent); i++) {
- if (SeparatorsDimens[R[i]] == nS) {
- contain_S = 0;
- for (j = 0; j < SeparatorsDimens[R[i]]; j++) {
- for (k = 0; k < nS; k++) {
- if (Separators[R[i]][j] == S[k]) {
- contain_S++;
- }
- }
- }
- if (contain_S == nS) {
- canadd = 1;
- }
- }
- if (canadd) {
- break;
- }
- }
-
- // else: check if S is in the path from T[0] to common_parent
- if (!canadd) {
- for (i = 0; i <= (pT - common_parent); i++) {
- if (SeparatorsDimens[T[i]] == nS) {
- contain_S = 0;
- for (j = 0; j < SeparatorsDimens[T[i]]; j++) {
- for (k = 0; k < nS; k++) {
- if (Separators[T[i]][j] == S[k]) {
- contain_S++;
- }
- }
- }
- if (contain_S == nS) {
- canadd = 1;
- }
- }
- if (canadd) {
- break;
- }
- }
- }
- }
- }
-
- delete[] R;
- delete[] T;
- delete[] S;
-
- return (canadd);
-}
-
-Real Graph::ScoreAddEdge(myInt a, myInt b, Real *D_prior, Real *D_post,
- myInt delta, myInt n_sub, Real score, int nEdges) {
- myInt j;
- int n = (int)nVertices;
- myInt nS = 0;
- myInt *S = new myInt[n];
-
- nS = 0;
- for (j = 0; j < n; j++) {
- if (Edge[a][j] && Edge[b][j]) {
- S[nS] = j;
- nS++;
- }
- }
-
- {
- myInt *C = new myInt[n];
- myInt nC;
- Real *sub_D = new Real[n * (n + 1) / 2];
- Real cScore;
-
- nC = nS + 2;
- for (j = 0; j < nS; j++) {
- C[j] = S[j];
- };
- C[nS] = a;
- C[nS + 1] = b;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(n, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score += cScore;
-
- nC = nS + 1;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(n, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score -= cScore;
-
- nC = nS + 1;
- C[nS] = b;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(n, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score -= cScore;
-
- nC = nS;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(n, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score += cScore; // cout << nS << "->";
-
- delete[] S;
- delete[] C;
- delete[] sub_D;
- return score;
- }
-}
-
-myInt Graph::CanDeleteEdge(myInt a, myInt b) {
- myInt i, j;
- myBool contain_a, contain_b;
- myInt count = 0, which_ab;
-
- for (i = 0; i < nCliques; i++) {
- contain_a = 0;
- contain_b = 0;
- for (j = 0; j < CliquesDimens[i]; j++) {
- if (Cliques[i][j] == a) {
- contain_a = 1;
- };
- if (Cliques[i][j] == b) {
- contain_b = 1;
- }
- }
- if (contain_a && contain_b) {
- which_ab = i;
- count++;
- }
- }
-
- if (count == 1) {
- return (which_ab);
- } else {
- return (-1);
- }
-}
-
-Real Graph::ScoreDeleteEdge(myInt a, myInt b, myInt which_ab, Real *D_prior,
- Real *D_post, myInt delta, myInt n_sub, Real score,
- int nEdges) {
- myInt j;
-
- {
- myInt p = nVertices;
- myInt *C = new myInt[p];
- myInt nC;
- Real *sub_D = new Real[p * (p + 1) / 2];
- Real cScore;
-
- nC = 0;
- for (j = 0; j < CliquesDimens[which_ab]; j++) {
- if ((Cliques[which_ab][j] != a) && (Cliques[which_ab][j] != b)) {
- C[nC] = Cliques[which_ab][j];
- nC++;
- }
- };
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(p, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score -= cScore;
-
- nC = CliquesDimens[which_ab] - 1;
- C[nC - 1] = a;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(p, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score += cScore;
-
- nC = CliquesDimens[which_ab] - 1;
- C[nC - 1] = b;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(p, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score += cScore;
-
- nC = CliquesDimens[which_ab];
- C[nC - 2] = a;
- C[nC - 1] = b;
- cScore = 0;
- cScore -= gwish_nc_complete(delta, nC, sub_D, 0);
- make_sub_mat_dbl(p, nC, C, D_post, sub_D);
- cScore += gwish_nc_complete(delta + n_sub, nC, sub_D, 1);
- score -= cScore;
-
- delete[] C;
- delete[] sub_D;
- return score;
- }
-}
\ No newline at end of file
diff --git a/src/sss-sycl/graph.h b/src/sss-sycl/graph.h
deleted file mode 100644
index f5cda3b912..0000000000
--- a/src/sss-sycl/graph.h
+++ /dev/null
@@ -1,117 +0,0 @@
-#define GRAPH_H
-
-typedef class Graph *LPGraph;
-
-class Graph {
-public:
- // data members
- myInt nVertices;
- myInt *d_nVertices; // number of vertices in the graph
- myInt **Edge;
- int *d_EdgeField; // matrix containing the edges of the graph
-
- myInt *Labels;
- myInt *d_Labels; // identifies the connected components of the graph
- myInt nLabels;
- myInt *d_nLabels; // number of labels or connected components
- myInt **Cliques; // storage for cliques
- myInt *CliquesDimens; // number of vertices in each clique
- myInt nCliques;
- myInt *d_nCliques; // number of cliques
-
-public:
- myInt *TreeEdgeA; // edges of the clique tree
- myInt *TreeEdgeB;
- myInt nTreeEdges; // number of edges in the generated clique tree
-public:
- myInt **Separators; // storage for separators
- myInt *SeparatorsDimens;
- myInt nSeparators;
- // private:
- myInt *localord;
-
- // methods
-public:
- Graph(); // constructor
- Graph(LPGraph InitialGraph); // constructor
- ~Graph(); // destructor
-public:
- myInt SearchVertex(); // identifies the next vertex to be eliminated
- void FlipEdge(myInt which); // flips an edge on a graph which is a myInteger
- // between 0 and the total number of edges
-
-public:
- // the MSS (Minimal Sufficient Statistics) are the maximal cliques for our
- // graph
- void InitGraph(myInt n);
- void CopyGraph(LPGraph G);
- void GenerateCliques(myInt label);
- myInt CheckCliques(myInt start,
- myInt end); // checks whether each generated component is
- // complete in the given graph
- myInt IsClique(
- myInt *vect,
- myInt nvect); // checks if the vertices in vect form a clique in our graph
-
- void GenerateSeparators();
- void AttachLabel(myInt v, myInt label);
- void GenerateLabels();
- myInt GenerateAllCliques();
- myInt IsDecomposable();
- myInt IfDecomposable();
-
- myInt CanDeleteEdge(myInt a, myInt b);
- bool CanAddEdge(myInt a, myInt b);
-
- Real ScoreDeleteEdge(myInt a, myInt b, myInt which_ab, Real *D_prior,
- Real *D_post, myInt delta, myInt n_sub, Real score,
- int nEdges);
- Real ScoreAddEdge(myInt a, myInt b, Real *D_prior, Real *D_post, myInt delta,
- myInt n_sub, Real score, int nEdges);
-};
-
-//////////////////////////////////////////////////////////////////////
-
-typedef class SectionGraph *LPSectionGraph;
-
-class SectionGraph : public Graph {
-public:
- myInt *
- Eliminated; // shows which vertices were eliminated from the initial graph
- myInt nEliminated; // number of vertices we eliminated
-
- // methods
-public:
- SectionGraph(LPGraph InitialGraph, myInt *velim); // constructor
- ~SectionGraph(); // destructor
-
-public:
- myInt IsChain(myInt u, myInt v); // see if there is a chain between u and v
- // or, equivalently, checks if u and v are in
- // the same connected component
-};
-
-////////////////////////////////////////////////////////////////////////
-
-typedef class EliminationGraph *LPEliminationGraph;
-
-class EliminationGraph : public Graph {
-public:
- myInt *
- Eliminated; // shows which vertices were eliminated from the initial graph
- myInt nEliminated; // number of vertices we eliminated
-
- // methods
-public:
- EliminationGraph(LPGraph InitialGraph, myInt vertex); // constructor
- ~EliminationGraph(); // destructor
-public:
- myInt SearchVertex(); // identify a vertex to be eliminated
-public:
- void EliminateVertex(myInt x); // eliminates an extra vertex
-};
-
-//////////////////////////////////////////////////////////////////////////
-
-// constructs the minimum fill-in graph for a nondecomposable graph
-void TurnFillInGraph(LPGraph graph);
\ No newline at end of file
diff --git a/src/sss-sycl/gwish.cpp b/src/sss-sycl/gwish.cpp
deleted file mode 100644
index 03c44cbafe..0000000000
--- a/src/sss-sycl/gwish.cpp
+++ /dev/null
@@ -1,168 +0,0 @@
-#include
-// gwish.cpp: This is a collection of device functions that manipulate graph
-// objects according to the G-Wishart distribution. Note: this depends on the
-// graph.cpp library Chiranjit Mukherjee : chiranjit@soe.ucsc.edu -- based on
-// Alex Lenkoski : lenkoski@stat.washington.edu
-
-#define GWISH_CPP
-#ifndef GRAPH_H
-#include "graph.h"
-#endif
-
-void log_det(int p, Real *A, Real *result) {
- int i, j, k;
- Real temp;
- *result = 0;
- Real br;
-
- for (i = 0; i < p; i++) {
- for (j = i; j < p; j++) {
- temp = A[j * (j + 1) / 2 + i];
- for (k = i; k > 0; k--) {
- temp = temp - A[j * (j + 1) / 2 + k - 1] * A[i * (i + 1) / 2 + k - 1];
- };
- if (i == j) {
- if (temp <= 0.0) {
- *result = NEG_INF;
- return;
- } else {
- br = sqrt(temp);
- }
- }
- A[j * (j + 1) / 2 + i] = temp / br;
- }
- }
-
-#ifdef ISFLOAT
- for (i = 0; i < p; i++) {
- *result += logf(A[i * (i + 1) / 2 + i]);
- };
- *result = 2 * (*result);
-#else
- for (i = 0; i < p; i++) {
- *result += log(A[i * (i + 1) / 2 + i]);
- };
- *result = 2 * (*result);
-#endif
-}
-
-// Computes the normalizing constant of a G-Wishart distribution for a full
-// p-dimensional graph with parameters delta and D : br p (reusable)
-Real gwish_nc_complete(myInt delta, int p, Real *D, bool flag) {
- Real d, c, a, g;
- Real dblDelta; // Recasting the inputs makes life easier below
- Real dblP;
- myInt i;
- c = 0.0;
- a = 0.0;
- g = 0.0;
- d = 0.0;
- dblDelta = delta;
- dblP = p;
-
- if (flag) {
- log_det(p, D, &d);
- }
- a = (dblDelta + dblP - 1) / 2.0;
- d = a * d;
- c = dblP * a * log_2;
- g = dblP * (dblP - 1) * log_pi_over_4;
-
- int signp;
-#ifdef ISFLOAT
- for (i = 0; i < p; i++)
- g += lgammaf_r(a - (Real)i / 2.0, &signp);
-#else
- for (i = 0; i < p; i++)
- g += lgamma_r(a - (Real)i / 2.0, &signp);
-#endif
-
- return (-d + c + g);
-}
-
-// Utility function for making submatrices
-void make_sub_mat_dbl(int p, int p_sub, myInt *sub, Real *A, Real *B) {
- int i, j;
- for (i = 0; i < p_sub; i++) {
- for (j = 0; j <= i; j++) {
- B[i * (i + 1) / 2 + j] =
- ((sub[i] >= sub[j]) ? A[sub[i] * (sub[i] + 1) / 2 + sub[j]]
- : A[sub[j] * (sub[j] + 1) / 2 + sub[i]]);
- }
- }
-}
-
-// Utility function for making a mean vector and a covariance matrix over a
-// subset of the dataset indicated by the vector sub.
-void make_sub_means_and_cov(Real *X, myInt *sub, myInt sub_match, int p,
- myInt n, myInt n_sub, Real *means, Real *D) {
- int i, j, k, ii;
- for (i = 0; i < p; i++) {
- means[i] = 0.0;
- };
- for (i = 0; i < p * (p + 1) / 2; i++) {
- D[i] = 0.0;
- }
-
- if (n_sub == 0) {
- return;
- }
- for (i = 0; i < p; i++) {
- for (k = 0; k < n; k++) {
- if (sub[k] == sub_match) {
- means[i] += (X[k * p + i] / (Real)n_sub);
- }
- }
- }
- if (n_sub < 2) {
- return;
- }
-
- for (i = 0; i < p; i++) {
- ii = i * (i + 1) / 2;
- for (j = 0; j <= i; j++) {
- for (k = 0; k < n; k++) {
- if (sub[k] == sub_match) {
- D[ii + j] += (X[k * p + i] - means[i]) * (X[k * p + j] - means[j]);
- }
- }
- }
- }
-
- return;
-}
-
-Real j_g_decomposable(LPGraph graph, Real *D_prior, Real *D_post, myInt delta,
- myInt n, bool flag) {
- Real mypost = 0;
- int p = graph->nVertices;
- myInt i;
- myInt sub_p;
- Real *sub_D = new Real[2 * p * p];
-
- //----- First loop through all the prime components (cliques since we're
- //decomposable) ---
- for (i = 0; i < graph->nCliques; i++) {
- sub_p = graph->CliquesDimens[i];
- if (flag) {
- make_sub_mat_dbl(p, sub_p, graph->Cliques[i], D_prior, sub_D);
- };
- mypost -= gwish_nc_complete(delta, sub_p, sub_D, flag);
- make_sub_mat_dbl(p, sub_p, graph->Cliques[i], D_post, sub_D);
- mypost += gwish_nc_complete(delta + n, sub_p, sub_D, 1);
- }
-
- //------- Now subtract off the separators -----------------------------------
- for (i = 0; i < graph->nSeparators; i++) {
- sub_p = graph->SeparatorsDimens[i];
- if (flag) {
- make_sub_mat_dbl(p, sub_p, graph->Separators[i], D_prior, sub_D);
- };
- mypost += gwish_nc_complete(delta, sub_p, sub_D, flag);
- make_sub_mat_dbl(p, sub_p, graph->Separators[i], D_post, sub_D);
- mypost -= gwish_nc_complete(delta + n, sub_p, sub_D, 1);
- }
-
- delete[] sub_D;
- return (mypost);
-}
\ No newline at end of file
diff --git a/src/sss-sycl/kernels.cpp b/src/sss-sycl/kernels.cpp
index 9ece307d00..0bce608b00 100644
--- a/src/sss-sycl/kernels.cpp
+++ b/src/sss-sycl/kernels.cpp
@@ -24,12 +24,15 @@ void CanDeleteEdge(myInt *d_in_delete, myInt *isDecomposable,
count = 0;
}
+ SYNC;
+
for (i = 0; i < nCliques; i++) {
ii = i * n;
if (tid == 0) {
contain_a = 0;
contain_b = 0;
}
+ SYNC;
for (j = tid; j < CliquesDimens[i]; j += bdim) {
k = Cliques[ii + j];
if (k == a) {
@@ -39,12 +42,16 @@ void CanDeleteEdge(myInt *d_in_delete, myInt *isDecomposable,
contain_b = 1;
}
}
+ SYNC;
if (tid == 0) {
if (contain_a && contain_b) {
count++;
which_ab = i;
}
}
+ // count is updated by thread 0 only, so it must be visible to every thread
+ // before the loop-exit test below, or threads leave the loop divergently.
+ SYNC;
if (count > 1) {
break;
}
@@ -145,6 +152,7 @@ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add, myInt *isDecomposable,
contain_b = 0;
};
contain_S[tid] = 0;
+ SYNC;
t = i * n;
for (j = tid; j < CliquesDimens[i]; j += bdim) {
c = Cliques[t + j];
@@ -160,6 +168,9 @@ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add, myInt *isDecomposable,
}
}
}
+ // every thread contributes to contain_S[] and contain_a/contain_b above,
+ // so those writes must be visible before thread 0 reduces them.
+ SYNC;
if (tid == 0) {
k = 0;
for (j = 0; j < BLOCK_SIZE; j++) {
@@ -172,6 +183,7 @@ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add, myInt *isDecomposable,
bSi = i;
}
}
+ SYNC;
}
if (tid == 0) { // find the path from aSi to root
@@ -208,6 +220,8 @@ void CanAddEdge(myInt *d_in_delete, myInt *d_in_add, myInt *isDecomposable,
}
}
}
+ // R/pR are produced by thread 0 and T/pT by thread 1; both are read below.
+ SYNC;
if (tid == 0) {
// find the branching point
diff --git a/src/sss-sycl/main.cpp b/src/sss-sycl/main.cpp
index b703b7236a..5607230a0c 100644
--- a/src/sss-sycl/main.cpp
+++ b/src/sss-sycl/main.cpp
@@ -66,13 +66,11 @@
#define ISFLOAT 10
using namespace std;
-#include
-#include
+#include "gsl_compat.h"
#define GSL_INTEGRATION_GRIDSIZE 1000
gsl_integration_workspace *w;
gsl_function F;
-#include
#define RANDOMSEED 314159265
// Define hyperparameters for the prior distribution of (mu, K | G)
@@ -271,6 +269,7 @@ int main(int argc, char *argv[]) {
for (l = 0; l < L; l++) {
score += state->pll[l];
}
+ k = 0;
printf("initial: k=%ld L=%d score=%.4f localBestScore=%.4f globalBestScore=%.4f "
"nmodes=%d num_cases=%d num_allModels=%ld\n",
k, state->L, score, localBestScore, globalBestScore, nmodes,
@@ -278,7 +277,6 @@ int main(int argc, char *argv[]) {
// start the stopwatch
auto start = std::chrono::steady_clock::now();
- k = 0;
while (nmodes <= maxNmodes) {
k++;
num_cases = 0;
@@ -352,7 +350,6 @@ int main(int argc, char *argv[]) {
modesList->UpdateList(globalBestState);
nmodes++;
gsl_rng_set(rnd, seedset[nmodes - 1]);
- start = std::chrono::steady_clock::now();
k = 0;
localBestScore = NEG_INF;
globalBestScore = NEG_INF;
diff --git a/src/sss-sycl/utilities.cpp b/src/sss-sycl/utilities.cpp
deleted file mode 100644
index 00280ff271..0000000000
--- a/src/sss-sycl/utilities.cpp
+++ /dev/null
@@ -1,74 +0,0 @@
-#define UTILITIES_CPP
-
-#ifndef GRAPH_CPP
-#include "graph.cpp"
-#endif
-
-// function to return a random myInteger with probability according to a given
-// (normalized!) weights
-myInt rand_myInt_weighted(myInt n, Real *weights) {
- myInt i;
- Real r = gsl_ran_flat(rnd, 0.0, 1.0);
- for (i = 0; i < n; i++) {
- if (r < weights[i]) {
- break;
- };
- r -= weights[i];
- };
- return i;
-}
-
-// returns a random myInterger between 0 and n - 1
-myInt rand_myInt(myInt n) {
- Real alpha = gsl_ran_flat(rnd, 0.0, 1.0);
- myInt value = (myInt)(n * alpha) / 1;
- return (value);
-}
-
-// function to return a random myInteger with probability according to a given
-// (normalized!) weights
-int rand_int_weighted(int n, Real *weights) {
- int i;
- Real r = gsl_ran_flat(rnd, 0.0, 1.0);
- for (i = 0; i < n; i++) {
- if (r < weights[i]) {
- break;
- };
- r -= weights[i];
- };
- return i;
-}
-
-// returns a random myInterger between 0 and n - 1
-int rand_int(int n) {
- Real alpha = gsl_ran_flat(rnd, 0.0, 1.0);
- int value = (int)(n * alpha) / 1;
- return (value);
-}
-
-// This samples an object from a subset of 0 to n - 1
-int sample_from(int *s, int n, int n_sub) {
- int i;
- int k = 0;
- int sample = rand_int(n_sub);
- for (i = 0; i < n; i++) {
- if (s[i]) {
- if (k == sample) {
- return (i);
- };
- k++;
- }
- }
-
- return (-1);
-}
-
-void shuffle(int *randomorder, int start, int end, int size) {
- int i, j, k;
- for (i = start; i < end; i++) {
- j = rand_int(size);
- k = randomorder[j];
- randomorder[j] = randomorder[i];
- randomorder[i] = k;
- }
-}
\ No newline at end of file
diff --git a/src/xlqc-cuda/CMakeLists.txt b/src/xlqc-cuda/CMakeLists.txt
index 122c738971..7f42adbab8 100644
--- a/src/xlqc-cuda/CMakeLists.txt
+++ b/src/xlqc-cuda/CMakeLists.txt
@@ -11,8 +11,8 @@ add_hecbench_benchmark(
basis.cc
int_lib/crys.cc
int_lib/cints.cc
- INCLUDE_DIRS ${SRC_DIR}/int_lib
+ INCLUDE_DIRS ${SRC_DIR} ${SRC_DIR}/int_lib
CATEGORIES algorithms
COMPILE_OPTIONS -dc
- LINK_LIBRARIES gsl gslcblas
+ LINK_LIBRARIES Eigen3::Eigen
)
diff --git a/src/xlqc-cuda/Makefile b/src/xlqc-cuda/Makefile
index 6c909a712a..73fff78817 100644
--- a/src/xlqc-cuda/Makefile
+++ b/src/xlqc-cuda/Makefile
@@ -7,10 +7,14 @@ CC = nvcc
OPTIMIZE = yes
DEBUG = no
ARCH = sm_60
-GSL_INC =-I/path/to/gsl/include
-GSL_LIB =-L/path/to/gsl/lib -lgsl -lgslcblas
+EIGEN_INC ?= $(shell pkg-config --cflags eigen3 2>/dev/null)
LAUNCHER =
+ifneq ($(MAKECMDGOALS),clean)
+ifeq ($(EIGEN_INC),)
+$(error Eigen not found; install libeigen3-dev or set EIGEN_INC=-I/path/to/eigen)
+endif
+endif
#===============================================================================
# Program name & source code list
@@ -26,10 +30,11 @@ obj=basis.o scf.o main.o crys.o cints.o cuda_rys_sp.o cuda_rys_dp.o
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Xcompiler -Wall $(GSL_INC) -arch=$(ARCH)
+CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Xcompiler -Wall $(EIGEN_INC) -arch=$(ARCH) \
+ --expt-relaxed-constexpr
# Linker Flags
-LDFLAGS = $(GSL_LIB)
+LDFLAGS =
# Debug Flags
ifeq ($(DEBUG),yes)
diff --git a/src/xlqc-cuda/README.md b/src/xlqc-cuda/README.md
index 7872cf0727..bac7b26a50 100644
--- a/src/xlqc-cuda/README.md
+++ b/src/xlqc-cuda/README.md
@@ -9,3 +9,14 @@ http://pyquante.sourceforge.net/
Refereces for Quantum Chemistry on GPU:
a) ERI evaluation. http://pubs.acs.org/doi/abs/10.1021/ct700268q
b) Direct SCF. http://pubs.acs.org/doi/abs/10.1021/ct800526s
+
+Dependencies
+------------
+The Hartree-Fock SCF driver needs a dense linear algebra library for the
+symmetric eigensolver, matrix products and the DIIS linear solve. These are
+provided by Eigen (header-only) through the thin `gsl_compat.h` wrapper in
+`xlqc-cuda/`. The Makefile discovers Eigen with `pkg-config --cflags eigen3`,
+so an Eigen installed by the system package manager (`libeigen3-dev` on Debian
+and Ubuntu, `eigen3-devel` on Fedora) needs no configuration. For an Eigen in a
+non-standard prefix, either add it to `PKG_CONFIG_PATH` or override the include
+path directly, for example `make EIGEN_INC=-I/opt/eigen-3.4.0/include/eigen3`.
diff --git a/src/xlqc-cuda/basis.cc b/src/xlqc-cuda/basis.cc
index 1ea31b674b..c794318cbc 100644
--- a/src/xlqc-cuda/basis.cc
+++ b/src/xlqc-cuda/basis.cc
@@ -24,9 +24,6 @@ of this software, even if advised of the possibility of such damage.
#include
-#include
-#include
-
#include "typedef.h"
#include "int_lib/cints.h"
#include "int_lib/crys.h"
diff --git a/src/xlqc-cuda/gsl_compat.h b/src/xlqc-cuda/gsl_compat.h
new file mode 100644
index 0000000000..5245c6509a
--- /dev/null
+++ b/src/xlqc-cuda/gsl_compat.h
@@ -0,0 +1,301 @@
+/*****************************************************************************
+ Eigen-backed replacement for the subset of the GNU Scientific Library used by
+ the XLQC benchmark, added so that the benchmark no longer links against the
+ GPL-licensed GSL (see HeCBench issue #319).
+
+ The names and signatures of the GSL entry points are kept so that the SCF code
+ reads unchanged; the dense linear algebra is delegated to Eigen (MPL2):
+
+ gsl_blas_dgemm -> Eigen matrix product
+ gsl_eigen_symmv -> Eigen::SelfAdjointEigenSolver
+ gsl_linalg_LU_decomp -> Eigen::PartialPivLU (partial pivoting, as in GSL)
+ gsl_linalg_LU_solve -> Eigen::PartialPivLU::solve
+
+ Only what XLQC needs is implemented. Two deliberate deviations from GSL, both
+ invisible to this benchmark: matrices are zero-filled on allocation rather
+ than left uninitialised, and gsl_eigen_symmv leaves its input intact instead
+ of destroying it.
+ *****************************************************************************/
+
+#ifndef GSL_COMPAT_H
+#define GSL_COMPAT_H
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+#define GSL_SUCCESS 0
+#define GSL_EFAILED 5
+#define GSL_EBADLEN 19
+
+//===============================
+// matrices and vectors
+//===============================
+
+typedef struct gsl_matrix_struct {
+ size_t size1;
+ size_t size2;
+ Eigen::MatrixXd m;
+ // Filled in by gsl_linalg_LU_decomp and consumed by gsl_linalg_LU_solve.
+ // GSL keeps the factorisation in the matrix itself, so it lives here too.
+ Eigen::PartialPivLU *lu;
+} gsl_matrix;
+
+typedef struct gsl_vector_struct {
+ size_t size;
+ Eigen::VectorXd v;
+} gsl_vector;
+
+inline gsl_matrix *gsl_matrix_alloc(const size_t n1, const size_t n2) {
+ gsl_matrix *a = new gsl_matrix;
+ a->size1 = n1;
+ a->size2 = n2;
+ a->m = Eigen::MatrixXd::Zero(n1, n2);
+ a->lu = NULL;
+ return a;
+}
+
+inline gsl_matrix *gsl_matrix_calloc(const size_t n1, const size_t n2) {
+ return gsl_matrix_alloc(n1, n2);
+}
+
+inline void gsl_matrix_free(gsl_matrix *a) {
+ if (a == NULL) return;
+ delete a->lu;
+ delete a;
+}
+
+inline double gsl_matrix_get(const gsl_matrix *a, const size_t i,
+ const size_t j) {
+ return a->m(i, j);
+}
+
+inline void gsl_matrix_set(gsl_matrix *a, const size_t i, const size_t j,
+ const double x) {
+ a->m(i, j) = x;
+}
+
+inline void gsl_matrix_set_zero(gsl_matrix *a) { a->m.setZero(); }
+
+inline int gsl_matrix_memcpy(gsl_matrix *dest, const gsl_matrix *src) {
+ if (dest->size1 != src->size1 || dest->size2 != src->size2)
+ return GSL_EBADLEN;
+ dest->m = src->m;
+ return GSL_SUCCESS;
+}
+
+inline gsl_vector *gsl_vector_alloc(const size_t n) {
+ gsl_vector *x = new gsl_vector;
+ x->size = n;
+ x->v = Eigen::VectorXd::Zero(n);
+ return x;
+}
+
+inline gsl_vector *gsl_vector_calloc(const size_t n) {
+ return gsl_vector_alloc(n);
+}
+
+inline void gsl_vector_free(gsl_vector *x) { delete x; }
+
+inline double gsl_vector_get(const gsl_vector *x, const size_t i) {
+ return x->v(i);
+}
+
+inline void gsl_vector_set(gsl_vector *x, const size_t i, const double y) {
+ x->v(i) = y;
+}
+
+inline void gsl_vector_set_zero(gsl_vector *x) { x->v.setZero(); }
+
+//===============================
+// BLAS level 3
+//===============================
+
+typedef enum {
+ CblasNoTrans = 111,
+ CblasTrans = 112,
+ CblasConjTrans = 113
+} CBLAS_TRANSPOSE_t;
+
+// C = alpha * op(A) * op(B) + beta * C
+inline int gsl_blas_dgemm(const CBLAS_TRANSPOSE_t TransA,
+ const CBLAS_TRANSPOSE_t TransB, const double alpha,
+ const gsl_matrix *A, const gsl_matrix *B,
+ const double beta, gsl_matrix *C) {
+ const bool ta = (TransA != CblasNoTrans);
+ const bool tb = (TransB != CblasNoTrans);
+
+ const size_t ma = ta ? A->size2 : A->size1;
+ const size_t na = ta ? A->size1 : A->size2;
+ const size_t mb = tb ? B->size2 : B->size1;
+ const size_t nb = tb ? B->size1 : B->size2;
+
+ if (na != mb || ma != C->size1 || nb != C->size2) return GSL_EBADLEN;
+
+ // Evaluated into a temporary so that C may alias A or B
+ Eigen::MatrixXd prod(ma, nb);
+ if (ta) {
+ if (tb)
+ prod.noalias() = A->m.transpose() * B->m.transpose();
+ else
+ prod.noalias() = A->m.transpose() * B->m;
+ } else {
+ if (tb)
+ prod.noalias() = A->m * B->m.transpose();
+ else
+ prod.noalias() = A->m * B->m;
+ }
+
+ if (beta == 0.0)
+ C->m = alpha * prod;
+ else
+ C->m = alpha * prod + beta * C->m;
+
+ return GSL_SUCCESS;
+}
+
+//===============================
+// symmetric eigenproblem
+//===============================
+
+typedef struct {
+ size_t size;
+} gsl_eigen_symmv_workspace;
+
+inline gsl_eigen_symmv_workspace *gsl_eigen_symmv_alloc(const size_t n) {
+ gsl_eigen_symmv_workspace *w = new gsl_eigen_symmv_workspace;
+ w->size = n;
+ return w;
+}
+
+inline void gsl_eigen_symmv_free(gsl_eigen_symmv_workspace *w) { delete w; }
+
+// Eigenvalues in eval, corresponding eigenvectors in the columns of evec.
+// Eigen already returns them in ascending order of eigenvalue.
+inline int gsl_eigen_symmv(gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
+ gsl_eigen_symmv_workspace *w) {
+ (void)w;
+
+ if (A->size1 != A->size2) return GSL_EBADLEN;
+ if (eval->size != A->size1 || evec->size1 != A->size1 ||
+ evec->size2 != A->size2)
+ return GSL_EBADLEN;
+
+ Eigen::SelfAdjointEigenSolver es(A->m);
+ if (es.info() != Eigen::Success) return GSL_EFAILED;
+
+ eval->v = es.eigenvalues();
+ evec->m = es.eigenvectors();
+
+ return GSL_SUCCESS;
+}
+
+typedef enum {
+ GSL_EIGEN_SORT_VAL_ASC,
+ GSL_EIGEN_SORT_VAL_DESC,
+ GSL_EIGEN_SORT_ABS_ASC,
+ GSL_EIGEN_SORT_ABS_DESC
+} gsl_eigen_sort_t;
+
+inline int gsl_eigen_symmv_sort(gsl_vector *eval, gsl_matrix *evec,
+ const gsl_eigen_sort_t sort_type) {
+ const size_t n = eval->size;
+ if (evec->size2 != n) return GSL_EBADLEN;
+
+ std::vector idx(n);
+ std::iota(idx.begin(), idx.end(), (size_t)0);
+
+ const Eigen::VectorXd &e = eval->v;
+ std::stable_sort(idx.begin(), idx.end(), [&](size_t i, size_t j) {
+ switch (sort_type) {
+ case GSL_EIGEN_SORT_VAL_ASC: return e(i) < e(j);
+ case GSL_EIGEN_SORT_VAL_DESC: return e(i) > e(j);
+ case GSL_EIGEN_SORT_ABS_ASC: return std::abs(e(i)) < std::abs(e(j));
+ default: return std::abs(e(i)) > std::abs(e(j));
+ }
+ });
+
+ Eigen::VectorXd eval_sorted(n);
+ Eigen::MatrixXd evec_sorted(evec->size1, n);
+ for (size_t k = 0; k < n; k++) {
+ eval_sorted(k) = eval->v(idx[k]);
+ evec_sorted.col(k) = evec->m.col(idx[k]);
+ }
+
+ eval->v = eval_sorted;
+ evec->m = evec_sorted;
+
+ return GSL_SUCCESS;
+}
+
+//===============================
+// LU decomposition
+//===============================
+
+typedef struct {
+ size_t size;
+ std::vector data;
+} gsl_permutation;
+
+inline gsl_permutation *gsl_permutation_alloc(const size_t n) {
+ gsl_permutation *p = new gsl_permutation;
+ p->size = n;
+ p->data.resize(n);
+ std::iota(p->data.begin(), p->data.end(), (size_t)0);
+ return p;
+}
+
+inline gsl_permutation *gsl_permutation_calloc(const size_t n) {
+ return gsl_permutation_alloc(n);
+}
+
+inline void gsl_permutation_free(gsl_permutation *p) { delete p; }
+
+inline size_t gsl_permutation_get(const gsl_permutation *p, const size_t i) {
+ return p->data[i];
+}
+
+// As in GSL, A is replaced by its LU factors and p by the row permutation.
+inline int gsl_linalg_LU_decomp(gsl_matrix *A, gsl_permutation *p, int *signum) {
+ *signum = 0;
+
+ if (A->size1 != A->size2 || p->size != A->size1) return GSL_EBADLEN;
+
+ delete A->lu;
+ A->lu = new Eigen::PartialPivLU(A->m);
+
+ const Eigen::PermutationMatrix &P =
+ A->lu->permutationP();
+
+ for (size_t i = 0; i < p->size; i++)
+ p->data[i] = (size_t)P.indices()(i);
+
+ *signum = (int)P.determinant();
+
+ A->m = A->lu->matrixLU();
+
+ return GSL_SUCCESS;
+}
+
+inline int gsl_linalg_LU_solve(const gsl_matrix *LU, const gsl_permutation *p,
+ const gsl_vector *b, gsl_vector *x) {
+ (void)p;
+
+ if (LU->lu == NULL) {
+ fprintf(stderr, "Error: gsl_linalg_LU_solve called before LU_decomp\n");
+ return GSL_EFAILED;
+ }
+ if (b->size != LU->size1 || x->size != LU->size2) return GSL_EBADLEN;
+
+ x->v = LU->lu->solve(b->v);
+
+ return GSL_SUCCESS;
+}
+
+#endif // GSL_COMPAT_H
diff --git a/src/xlqc-cuda/main.cu b/src/xlqc-cuda/main.cu
index b31adf2a4a..2361118b70 100644
--- a/src/xlqc-cuda/main.cu
+++ b/src/xlqc-cuda/main.cu
@@ -2,7 +2,6 @@
This file is part of the XLQC program.
Copyright (C) 2015 Xin Li
- Filename: main.cu
License: BSD 3-Clause License
This software is provided by the copyright holders and contributors "as is"
@@ -18,7 +17,6 @@
*****************************************************************************/
#include
-#include
#include
#include
#include
@@ -26,11 +24,7 @@
#include
#include
-#include
-#include
-#include
-#include
-#include
+#include "gsl_compat.h"
#include "int_lib/cints.h"
#include "int_lib/crys.h"
@@ -48,6 +42,8 @@ int main(int argc, char* argv[])
int use_5d = 1;
// use double precision?
int use_dp = 1;
+ // status is set on failure of result check
+ int status = 0;
if (argc > 1) {
for (int i = 1; i < argc; ++ i) {
@@ -590,29 +586,27 @@ int main(int argc, char* argv[])
const double err = fabs(ene_total - ref_energy);
if (err < tol) {
fprintf(stdout, "PASS: E_total error %.2e within tolerance %.0e\n", err, tol);
+ // print MO information
+ fprintf(stdout, "%5s %10s %15s %12s\n", "MO", "State", "E(Eh)", "E(eV)");
+ for (ibasis = 0; ibasis < p_basis->num; ++ ibasis)
+ {
+ char occ[10];
+ if (ibasis < n_occ) { strcpy(occ, "occ."); }
+ else { strcpy(occ, "virt."); }
+
+ double ener = gsl_vector_get(emo, ibasis);
+ fprintf(stdout, "%5d %10s %15.5f %12.2f\n",
+ ibasis + 1, occ, ener, ener * HARTREE2EV);
+ }
} else {
fprintf(stderr, "FAIL: E_total = %.10f, expected %.10f (error %.2e > tol %.0e)\n",
ene_total, ref_energy, err, tol);
- return 1;
- }
-
- // print MO information
- start = std::chrono::steady_clock::now();
-
- fprintf(stdout, "%5s %10s %15s %12s\n", "MO", "State", "E(Eh)", "E(eV)");
- for (ibasis = 0; ibasis < p_basis->num; ++ ibasis)
- {
- char occ[10];
- if (ibasis < n_occ) { strcpy(occ, "occ."); }
- else { strcpy(occ, "virt."); }
-
- double ener = gsl_vector_get(emo, ibasis);
- fprintf(stdout, "%5d %10s %15.5f %12.2f\n",
- ibasis + 1, occ, ener, ener * HARTREE2EV);
+ status = 1;
}
//====== free device memories ========
+ start = std::chrono::steady_clock::now();
cudaFree(dev_pbf_xlec);
cudaFree(dev_pbf_to_cbf);
@@ -716,5 +710,5 @@ int main(int argc, char* argv[])
//====== the end of program ========
- return 0;
+ return status;
}
diff --git a/src/xlqc-cuda/scf.cc b/src/xlqc-cuda/scf.cc
index 31730f901b..ca367d95a1 100644
--- a/src/xlqc-cuda/scf.cc
+++ b/src/xlqc-cuda/scf.cc
@@ -22,11 +22,7 @@
#include
#include
-#include
-#include
-#include
-#include
-#include
+#include "gsl_compat.h"
#include "typedef.h"
#include "basis.h"
diff --git a/src/xlqc-hip/CMakeLists.txt b/src/xlqc-hip/CMakeLists.txt
index 48a9cc16ac..54300b62f5 100644
--- a/src/xlqc-hip/CMakeLists.txt
+++ b/src/xlqc-hip/CMakeLists.txt
@@ -13,8 +13,8 @@ add_hecbench_benchmark(
${INC_DIR}/basis.cc
${INC_DIR}/int_lib/crys.cc
${INC_DIR}/int_lib/cints.cc
- INCLUDE_DIRS ${GSL_INCLUDE_DIRS} ${INC_DIR} ${INC_DIR}/int_lib
+ INCLUDE_DIRS ${INC_DIR} ${INC_DIR}/int_lib
CATEGORIES algorithms
COMPILE_OPTIONS -fgpu-rdc
- LINK_LIBRARIES -fgpu-rdc --hip-link gsl gslcblas
+ LINK_LIBRARIES -fgpu-rdc --hip-link Eigen3::Eigen
)
diff --git a/src/xlqc-hip/Makefile b/src/xlqc-hip/Makefile
index 878030f10a..b6dce1e676 100644
--- a/src/xlqc-hip/Makefile
+++ b/src/xlqc-hip/Makefile
@@ -6,10 +6,14 @@
CC = hipcc
OPTIMIZE = yes
DEBUG = no
-GSL_INC =-I/path/to/gsl/include
-GSL_LIB =-L/path/to/gsl/lib -lgsl -lgslcblas
+EIGEN_INC ?= $(shell pkg-config --cflags eigen3 2>/dev/null)
LAUNCHER =
+ifneq ($(MAKECMDGOALS),clean)
+ifeq ($(EIGEN_INC),)
+$(error Eigen not found; install libeigen3-dev or set EIGEN_INC=-I/path/to/eigen)
+endif
+endif
#===============================================================================
# Program name & source code list
@@ -25,10 +29,10 @@ obj=basis.o scf.o main.o crys.o cints.o cuda_rys_sp.o cuda_rys_dp.o
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../xlqc-cuda -I../xlqc-cuda/int_lib $(GSL_INC)
+CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../xlqc-cuda -I../xlqc-cuda/int_lib $(EIGEN_INC)
# Linker Flags
-LDFLAGS = $(GSL_LIB) -fgpu-rdc --hip-link
+LDFLAGS = -fgpu-rdc --hip-link
# Debug Flags
ifeq ($(DEBUG),yes)
@@ -53,16 +57,16 @@ $(program): $(obj)
$(CC) -fgpu-rdc -c -o $@ $< $(CFLAGS)
scf.o: ../xlqc-cuda/scf.cc
- $(CC) -c -o $@ $< $(CFLAGS)
+ $(CC) -fgpu-rdc -c -o $@ $< $(CFLAGS)
basis.o: ../xlqc-cuda/basis.cc
- $(CC) -c -o $@ $< $(CFLAGS)
+ $(CC) -fgpu-rdc -c -o $@ $< $(CFLAGS)
crys.o: ../xlqc-cuda/int_lib/crys.cc
- $(CC) -c -o $@ $< $(CFLAGS)
+ $(CC) -fgpu-rdc -c -o $@ $< $(CFLAGS)
cints.o: ../xlqc-cuda/int_lib/cints.cc
- $(CC) -c -o $@ $< $(CFLAGS)
+ $(CC) -fgpu-rdc -c -o $@ $< $(CFLAGS)
clean:
rm -rf $(program) $(obj)
diff --git a/src/xlqc-hip/README.md b/src/xlqc-hip/README.md
index 7872cf0727..bac7b26a50 100644
--- a/src/xlqc-hip/README.md
+++ b/src/xlqc-hip/README.md
@@ -9,3 +9,14 @@ http://pyquante.sourceforge.net/
Refereces for Quantum Chemistry on GPU:
a) ERI evaluation. http://pubs.acs.org/doi/abs/10.1021/ct700268q
b) Direct SCF. http://pubs.acs.org/doi/abs/10.1021/ct800526s
+
+Dependencies
+------------
+The Hartree-Fock SCF driver needs a dense linear algebra library for the
+symmetric eigensolver, matrix products and the DIIS linear solve. These are
+provided by Eigen (header-only) through the thin `gsl_compat.h` wrapper in
+`xlqc-cuda/`. The Makefile discovers Eigen with `pkg-config --cflags eigen3`,
+so an Eigen installed by the system package manager (`libeigen3-dev` on Debian
+and Ubuntu, `eigen3-devel` on Fedora) needs no configuration. For an Eigen in a
+non-standard prefix, either add it to `PKG_CONFIG_PATH` or override the include
+path directly, for example `make EIGEN_INC=-I/opt/eigen-3.4.0/include/eigen3`.
diff --git a/src/xlqc-hip/main.cu b/src/xlqc-hip/main.cu
index 5d692e4990..37fab0a1f7 100644
--- a/src/xlqc-hip/main.cu
+++ b/src/xlqc-hip/main.cu
@@ -2,7 +2,6 @@
This file is part of the XLQC program.
Copyright (C) 2015 Xin Li
- Filename: main.cu
License: BSD 3-Clause License
This software is provided by the copyright holders and contributors "as is"
@@ -18,7 +17,6 @@
*****************************************************************************/
#include
-#include
#include
#include
#include
@@ -26,11 +24,7 @@
#include
#include
-#include
-#include
-#include
-#include
-#include
+#include "gsl_compat.h"
#include "int_lib/cints.h"
#include "int_lib/crys.h"
@@ -48,6 +42,8 @@ int main(int argc, char* argv[])
int use_5d = 1;
// use double precision?
int use_dp = 1;
+ // status is set on failure of result check
+ int status = 0;
if (argc > 1) {
for (int i = 1; i < argc; ++ i) {
@@ -586,29 +582,26 @@ int main(int argc, char* argv[])
const double err = fabs(ene_total - ref_energy);
if (err < tol) {
fprintf(stdout, "PASS: E_total error %.2e within tolerance %.0e\n", err, tol);
+ // print MO information
+ fprintf(stdout, "%5s %10s %15s %12s\n", "MO", "State", "E(Eh)", "E(eV)");
+ for (ibasis = 0; ibasis < p_basis->num; ++ ibasis)
+ {
+ char occ[10];
+ if (ibasis < n_occ) { strcpy(occ, "occ."); }
+ else { strcpy(occ, "virt."); }
+
+ double ener = gsl_vector_get(emo, ibasis);
+ fprintf(stdout, "%5d %10s %15.5f %12.2f\n",
+ ibasis + 1, occ, ener, ener * HARTREE2EV);
+ }
} else {
fprintf(stderr, "FAIL: E_total = %.10f, expected %.10f (error %.2e > tol %.0e)\n",
ene_total, ref_energy, err, tol);
- return 1;
+ status = 1;
}
- // print MO information
- start = std::chrono::steady_clock::now();
-
- fprintf(stdout, "%5s %10s %15s %12s\n", "MO", "State", "E(Eh)", "E(eV)");
- for (ibasis = 0; ibasis < p_basis->num; ++ ibasis)
- {
- char occ[10];
- if (ibasis < n_occ) { strcpy(occ, "occ."); }
- else { strcpy(occ, "virt."); }
-
- double ener = gsl_vector_get(emo, ibasis);
- fprintf(stdout, "%5d %10s %15.5f %12.2f\n",
- ibasis + 1, occ, ener, ener * HARTREE2EV);
- }
-
-
//====== free device memories ========
+ start = std::chrono::steady_clock::now();
hipFree(dev_pbf_xlec);
hipFree(dev_pbf_to_cbf);
@@ -712,5 +705,5 @@ int main(int argc, char* argv[])
//====== the end of program ========
- return 0;
+ return status;
}
diff --git a/src/xlqc-omp/CMakeLists.txt b/src/xlqc-omp/CMakeLists.txt
index d7fc40129e..455ef10ace 100644
--- a/src/xlqc-omp/CMakeLists.txt
+++ b/src/xlqc-omp/CMakeLists.txt
@@ -1,8 +1,20 @@
# xlqc-omp/CMakeLists.txt
+set(SRC_DIR "${CMAKE_CURRENT_LIST_DIR}")
+UNZIPFILE("${SRC_DIR}/example.tar.gz")
+
+set(INC_DIR "${SRC_DIR}/../xlqc-cuda")
+
+# cuda_rys_sp.cpp and cuda_rys_dp.cpp are #included by main.cpp
add_hecbench_benchmark(
NAME xlqc
MODEL omp
- SOURCES cuda_rys_dp.cpp cuda_rys_sp.cpp main.cpp
+ SOURCES main.cpp
+ ${INC_DIR}/scf.cc
+ ${INC_DIR}/basis.cc
+ ${INC_DIR}/int_lib/crys.cc
+ ${INC_DIR}/int_lib/cints.cc
+ INCLUDE_DIRS ${INC_DIR} ${INC_DIR}/int_lib
CATEGORIES algorithms
+ LINK_LIBRARIES Eigen3::Eigen
)
diff --git a/src/xlqc-omp/Makefile b/src/xlqc-omp/Makefile
index d472ee8fd2..59811ab234 100644
--- a/src/xlqc-omp/Makefile
+++ b/src/xlqc-omp/Makefile
@@ -7,10 +7,14 @@ CC = icpx
OPTIMIZE = yes
DEBUG = no
DEVICE = gpu
-GSL_INC =-I/path/to/gsl/include
-GSL_LIB =-L/path/to/gsl/lib -lgsl -lgslcblas
+EIGEN_INC ?= $(shell pkg-config --cflags eigen3 2>/dev/null)
LAUNCHER =
+ifneq ($(MAKECMDGOALS),clean)
+ifeq ($(EIGEN_INC),)
+$(error Eigen not found; install libeigen3-dev or set EIGEN_INC=-I/path/to/eigen)
+endif
+endif
#===============================================================================
# Program name & source code list
@@ -26,10 +30,11 @@ obj=basis.o scf.o main.o crys.o cints.o
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../xlqc-cuda -I../xlqc-cuda/int_lib $(GSL_INC)
+CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall \
+ -I../xlqc-cuda -I../xlqc-cuda/int_lib $(EIGEN_INC)
# Linker Flags
-LDFLAGS = $(GSL_LIB)
+LDFLAGS =
# Debug Flags
ifeq ($(DEBUG),yes)
diff --git a/src/xlqc-omp/Makefile.aomp b/src/xlqc-omp/Makefile.aomp
index 68a31576e9..1a5356f3df 100644
--- a/src/xlqc-omp/Makefile.aomp
+++ b/src/xlqc-omp/Makefile.aomp
@@ -8,10 +8,14 @@ OPTIMIZE = yes
DEBUG = no
DEVICE = gpu
ARCH = gfx906
-GSL_INC =-I/path/to/gsl/include
-GSL_LIB =-L/path/to/gsl/lib -lgsl -lgslcblas
+EIGEN_INC ?= $(shell pkg-config --cflags eigen3 2>/dev/null)
LAUNCHER =
+ifneq ($(MAKECMDGOALS),clean)
+ifeq ($(EIGEN_INC),)
+$(error Eigen not found; install libeigen3-dev or set EIGEN_INC=-I/path/to/eigen)
+endif
+endif
#===============================================================================
# Program name & source code list
@@ -27,10 +31,11 @@ obj=basis.o scf.o main.o crys.o cints.o
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall -I../xlqc-cuda -I../xlqc-cuda/int_lib $(GSL_INC)
+CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Wall \
+ -I../xlqc-cuda -I../xlqc-cuda/int_lib $(EIGEN_INC)
# Linker Flags
-LDFLAGS = $(GSL_LIB)
+LDFLAGS =
# Debug Flags
ifeq ($(DEBUG),yes)
diff --git a/src/xlqc-omp/Makefile.nvc b/src/xlqc-omp/Makefile.nvc
index 943c302800..54697060e6 100644
--- a/src/xlqc-omp/Makefile.nvc
+++ b/src/xlqc-omp/Makefile.nvc
@@ -8,10 +8,14 @@ OPTIMIZE = yes
DEBUG = no
DEVICE = gpu
SM = cc70
-GSL_INC =-I/path/to/gsl/include
-GSL_LIB =-L/path/to/gsl/lib -lgsl -lgslcblas
+EIGEN_INC ?= $(shell pkg-config --cflags eigen3 2>/dev/null)
LAUNCHER =
+ifneq ($(MAKECMDGOALS),clean)
+ifeq ($(EIGEN_INC),)
+$(error Eigen not found; install libeigen3-dev or set EIGEN_INC=-I/path/to/eigen)
+endif
+endif
#===============================================================================
# Program name & source code list
@@ -27,10 +31,11 @@ obj=basis.o scf.o main.o crys.o cints.o
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -std=c++14 -Wall -I../xlqc-cuda -I../xlqc-cuda/int_lib $(GSL_INC)
+CFLAGS := $(EXTRA_CFLAGS) -std=c++14 -Wall \
+ -I../xlqc-cuda -I../xlqc-cuda/int_lib $(EIGEN_INC)
# Linker Flags
-LDFLAGS = $(GSL_LIB)
+LDFLAGS =
# Debug Flags
ifeq ($(DEBUG),yes)
diff --git a/src/xlqc-omp/README.md b/src/xlqc-omp/README.md
index 7872cf0727..bac7b26a50 100644
--- a/src/xlqc-omp/README.md
+++ b/src/xlqc-omp/README.md
@@ -9,3 +9,14 @@ http://pyquante.sourceforge.net/
Refereces for Quantum Chemistry on GPU:
a) ERI evaluation. http://pubs.acs.org/doi/abs/10.1021/ct700268q
b) Direct SCF. http://pubs.acs.org/doi/abs/10.1021/ct800526s
+
+Dependencies
+------------
+The Hartree-Fock SCF driver needs a dense linear algebra library for the
+symmetric eigensolver, matrix products and the DIIS linear solve. These are
+provided by Eigen (header-only) through the thin `gsl_compat.h` wrapper in
+`xlqc-cuda/`. The Makefile discovers Eigen with `pkg-config --cflags eigen3`,
+so an Eigen installed by the system package manager (`libeigen3-dev` on Debian
+and Ubuntu, `eigen3-devel` on Fedora) needs no configuration. For an Eigen in a
+non-standard prefix, either add it to `PKG_CONFIG_PATH` or override the include
+path directly, for example `make EIGEN_INC=-I/opt/eigen-3.4.0/include/eigen3`.
diff --git a/src/xlqc-omp/main.cpp b/src/xlqc-omp/main.cpp
index 7c3813ebf2..ac62aa6bd7 100644
--- a/src/xlqc-omp/main.cpp
+++ b/src/xlqc-omp/main.cpp
@@ -2,7 +2,6 @@
This file is part of the XLQC program.
Copyright (C) 2015 Xin Li
- Filename: main.cu
License: BSD 3-Clause License
This software is provided by the copyright holders and contributors "as is"
@@ -18,22 +17,16 @@
*****************************************************************************/
#include
-#include
#include
#include
#include
#include
#include
#include
-
-#include
-#include
-#include
-#include
-#include
-
#include
+#include "gsl_compat.h"
+
#include "int_lib/cints.h"
#include "int_lib/crys.h"
@@ -51,6 +44,8 @@ int main(int argc, char* argv[])
int use_5d = 1;
// use double precision?
int use_dp = 1;
+ // status is set on failure of result check
+ int status = 0;
if (argc > 1) {
for (int i = 1; i < argc; ++ i) {
@@ -571,27 +566,26 @@ int main(int argc, char* argv[])
const double err = fabs(ene_total - ref_energy);
if (err < tol) {
fprintf(stdout, "PASS: E_total error %.2e within tolerance %.0e\n", err, tol);
+ // print MO information
+ fprintf(stdout, "%5s %10s %15s %12s\n", "MO", "State", "E(Eh)", "E(eV)");
+ for (ibasis = 0; ibasis < p_basis->num; ++ ibasis)
+ {
+ char occ[10];
+ if (ibasis < n_occ) { strcpy(occ, "occ."); }
+ else { strcpy(occ, "virt."); }
+
+ double ener = gsl_vector_get(emo, ibasis);
+ fprintf(stdout, "%5d %10s %15.5f %12.2f\n",
+ ibasis + 1, occ, ener, ener * HARTREE2EV);
+ }
} else {
fprintf(stderr, "FAIL: E_total = %.10f, expected %.10f (error %.2e > tol %.0e)\n",
ene_total, ref_energy, err, tol);
- return 1;
+ status = 1;
}
- // print MO information
start = std::chrono::steady_clock::now();
- fprintf(stdout, "%5s %10s %15s %12s\n", "MO", "State", "E(Eh)", "E(eV)");
- for (ibasis = 0; ibasis < p_basis->num; ++ ibasis)
- {
- char occ[10];
- if (ibasis < n_occ) { strcpy(occ, "occ."); }
- else { strcpy(occ, "virt."); }
-
- double ener = gsl_vector_get(emo, ibasis);
- fprintf(stdout, "%5d %10s %15.5f %12.2f\n",
- ibasis + 1, occ, ener, ener * HARTREE2EV);
- }
-
} // omp target
//====== free host memories ========
@@ -688,5 +682,5 @@ int main(int argc, char* argv[])
//====== the end of program ========
- return 0;
+ return status;
}
diff --git a/src/xlqc-sycl/CMakeLists.txt b/src/xlqc-sycl/CMakeLists.txt
index 26e3ed6265..a32f920aaa 100644
--- a/src/xlqc-sycl/CMakeLists.txt
+++ b/src/xlqc-sycl/CMakeLists.txt
@@ -13,7 +13,7 @@ add_hecbench_benchmark(
${INC_DIR}/basis.cc
${INC_DIR}/int_lib/crys.cc
${INC_DIR}/int_lib/cints.cc
- INCLUDE_DIRS ${GSL_INCLUDE_DIRS} ${INC_DIR} ${INC_DIR}/int_lib
+ INCLUDE_DIRS ${INC_DIR} ${INC_DIR}/int_lib
CATEGORIES algorithms
- LINK_LIBRARIES gsl gslcblas
+ LINK_LIBRARIES Eigen3::Eigen
)
diff --git a/src/xlqc-sycl/Makefile b/src/xlqc-sycl/Makefile
index e34874aaa1..2dd36f61f6 100644
--- a/src/xlqc-sycl/Makefile
+++ b/src/xlqc-sycl/Makefile
@@ -6,8 +6,7 @@
CC = clang++
OPTIMIZE = yes
DEBUG = no
-GSL_INC =-I/path/to/gsl/include
-GSL_LIB =-L/path/to/gsl/lib -lgsl -lgslcblas
+EIGEN_INC ?= $(shell pkg-config --cflags eigen3 2>/dev/null)
LAUNCHER =
GPU = yes
@@ -17,6 +16,12 @@ HIP = no
HIP_ARCH = gfx908
#GCC_TOOLCHAIN = "/auto/software/gcc/x86_64/gcc-9.1.0/"
+ifneq ($(MAKECMDGOALS),clean)
+ifeq ($(EIGEN_INC),)
+$(error Eigen not found; install libeigen3-dev or set EIGEN_INC=-I/path/to/eigen)
+endif
+endif
+
#===============================================================================
# Program name & source code list
#===============================================================================
@@ -31,7 +36,7 @@ obj=basis.o scf.o main.o crys.o cints.o
#===============================================================================
# Standard Flags
-CFLAGS := $(EXTRA_CFLAGS) -I../xlqc-cuda -I../xlqc-cuda/int_lib $(GSL_INC) \
+CFLAGS := $(EXTRA_CFLAGS) -I../xlqc-cuda -I../xlqc-cuda/int_lib $(EIGEN_INC) \
--gcc-toolchain=$(GCC_TOOLCHAIN) \
-std=c++17 -Wall -fsycl
@@ -41,7 +46,7 @@ ifeq ($(VENDOR), AdaptiveCpp)
endif
# Linker Flags
-LDFLAGS = $(GSL_LIB)
+LDFLAGS =
ifeq ($(CUDA), yes)
CFLAGS += -fsycl-targets=nvptx64-nvidia-cuda \
@@ -92,6 +97,8 @@ cints.o: ../xlqc-cuda/int_lib/cints.cc
clean:
rm -rf $(program) $(obj)
+# ONEAPI_DEVICE_SELECTOR=level_zero:gpu
+# ONEAPI_DEVICE_SELECTOR=level_zero:gpu OverrideDefaultFP64Settings=1 IGC_EnableDPEmulation=1
run: $(program)
$(LAUNCHER) ./$(program) sp
$(LAUNCHER) ./$(program) dp
diff --git a/src/xlqc-sycl/README.md b/src/xlqc-sycl/README.md
index 7872cf0727..bac7b26a50 100644
--- a/src/xlqc-sycl/README.md
+++ b/src/xlqc-sycl/README.md
@@ -9,3 +9,14 @@ http://pyquante.sourceforge.net/
Refereces for Quantum Chemistry on GPU:
a) ERI evaluation. http://pubs.acs.org/doi/abs/10.1021/ct700268q
b) Direct SCF. http://pubs.acs.org/doi/abs/10.1021/ct800526s
+
+Dependencies
+------------
+The Hartree-Fock SCF driver needs a dense linear algebra library for the
+symmetric eigensolver, matrix products and the DIIS linear solve. These are
+provided by Eigen (header-only) through the thin `gsl_compat.h` wrapper in
+`xlqc-cuda/`. The Makefile discovers Eigen with `pkg-config --cflags eigen3`,
+so an Eigen installed by the system package manager (`libeigen3-dev` on Debian
+and Ubuntu, `eigen3-devel` on Fedora) needs no configuration. For an Eigen in a
+non-standard prefix, either add it to `PKG_CONFIG_PATH` or override the include
+path directly, for example `make EIGEN_INC=-I/opt/eigen-3.4.0/include/eigen3`.
diff --git a/src/xlqc-sycl/main.cpp b/src/xlqc-sycl/main.cpp
index a151334254..c88ca6b165 100644
--- a/src/xlqc-sycl/main.cpp
+++ b/src/xlqc-sycl/main.cpp
@@ -17,7 +17,6 @@ of this software, even if advised of the possibility of such damage.
*****************************************************************************/
#include
-#include
#include
#include
#include
@@ -25,12 +24,9 @@ of this software, even if advised of the possibility of such damage.
#include
#include
+#include
-#include
-#include
-#include
-#include
-#include
+#include "gsl_compat.h"
#include "int_lib/cints.h"
#include "int_lib/crys.h"
@@ -39,7 +35,6 @@ of this software, even if advised of the possibility of such damage.
#include "basis.h"
#include "scf.h"
-#include
#include "rys.h"
#include "cuda_rys_sp.cpp"
#include "cuda_rys_dp.cpp"
@@ -50,6 +45,8 @@ int main(int argc, char* argv[])
int use_5d = 1;
// use double precision?
int use_dp = 1;
+ // status is set on failure of result check
+ int status = 0;
if (argc > 1) {
for (int i = 1; i < argc; ++ i) {
@@ -638,28 +635,26 @@ int main(int argc, char* argv[])
const double err = fabs(ene_total - ref_energy);
if (err < tol) {
fprintf(stdout, "PASS: E_total error %.2e within tolerance %.0e\n", err, tol);
+ // print MO information
+ fprintf(stdout, "%5s %10s %15s %12s\n", "MO", "State", "E(Eh)", "E(eV)");
+ for (ibasis = 0; ibasis < p_basis->num; ++ ibasis)
+ {
+ char occ[10];
+ if (ibasis < n_occ) { strcpy(occ, "occ."); }
+ else { strcpy(occ, "virt."); }
+
+ double ener = gsl_vector_get(emo, ibasis);
+ fprintf(stdout, "%5d %10s %15.5f %12.2f\n",
+ ibasis + 1, occ, ener, ener * HARTREE2EV);
+ }
} else {
fprintf(stderr, "FAIL: E_total = %.10f, expected %.10f (error %.2e > tol %.0e)\n",
ene_total, ref_energy, err, tol);
- return 1;
- }
-
- // print MO information
- start = std::chrono::steady_clock::now();
-
- fprintf(stdout, "%5s %10s %15s %12s\n", "MO", "State", "E(Eh)", "E(eV)");
- for (ibasis = 0; ibasis < p_basis->num; ++ ibasis)
- {
- char occ[10];
- if (ibasis < n_occ) { strcpy(occ, "occ."); }
- else { strcpy(occ, "virt."); }
-
- double ener = gsl_vector_get(emo, ibasis);
- fprintf(stdout, "%5d %10s %15.5f %12.2f\n",
- ibasis + 1, occ, ener, ener * HARTREE2EV);
+ status = 1;
}
//====== free device memories ========
+ start = std::chrono::steady_clock::now();
sycl::free(d_pbf_xlec, q);
sycl::free(d_pbf_to_cbf, q);
@@ -762,5 +757,5 @@ int main(int argc, char* argv[])
//====== the end of program ========
- return 0;
+ return status;
}