From 24a04d755b3bdfe04d2499baf9f25f3ded24fca9 Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Tue, 11 Aug 2026 16:05:10 -0700 Subject: [PATCH 1/4] [xlqc] replace GSL with Eigen and discover Eigen instead of hardcoding its path GSL is GPL, so a redistributed binary that links it together with the CUDA or HIP runtime cannot satisfy both licenses. Reported in https://github.com/ORNL/HeCBench/issues/319 The SCF driver only needs a handful of GSL entry points: matrix and vector allocation and element access, dgemm, the symmetric eigensolver, and an LU solve for DIIS. gsl_compat.h maps exactly those onto Eigen (MPL2, header-only), keeping the existing gsl_* call sites in main, scf.cc and basis.cc unchanged. The header lives in xlqc-cuda and is shared by the hip, omp and sycl variants, which already compile those sources through -I../xlqc-cuda. Eigen's SelfAdjointEigenSolver returns ascending eigenvalues like gsl_eigen_symmv_sort, and PartialPivLU matches the pivoting of gsl_linalg_LU_decomp, so the SCF trajectory is preserved. The Makefiles carried the same "/path/to/..." placeholder that https://github.com/ORNL/HeCBench/pull/305 replaces for GSL, so apply that fix to Eigen: EIGEN_INC defaults to `pkg-config --cflags eigen3` and a missing Eigen stops the build with an actionable message instead of a bare missing-header error. Eigen is header-only, so there is no library counterpart to GSL_LIB and LDFLAGS becomes empty. The check is skipped for `clean` so the tree can be cleaned without Eigen installed. CMake gains find_package(Eigen3) and guards xlqc on the Eigen3::Eigen imported target rather than a _FOUND variable, since that target is what carries the include path. xlqc-omp additionally never compiled the shared scf.cc, basis.cc and int_lib sources or unpacked its example data, so it could not link; both are fixed here. Co-authored-by: Cursor --- CMakeLists.txt | 1 + README.md | 5 +- cmake/modules/BenchmarkMacros.cmake | 13 +- src/xlqc-cuda/CMakeLists.txt | 4 +- src/xlqc-cuda/Makefile | 12 +- src/xlqc-cuda/README.md | 11 + src/xlqc-cuda/basis.cc | 3 - src/xlqc-cuda/gsl_compat.h | 300 ++++++++++++++++++++++++++++ src/xlqc-cuda/main.cu | 35 ++-- src/xlqc-cuda/scf.cc | 6 +- src/xlqc-hip/CMakeLists.txt | 4 +- src/xlqc-hip/Makefile | 12 +- src/xlqc-hip/README.md | 11 + src/xlqc-hip/main.cu | 36 ++-- src/xlqc-omp/CMakeLists.txt | 14 +- src/xlqc-omp/Makefile | 13 +- src/xlqc-omp/Makefile.aomp | 13 +- src/xlqc-omp/Makefile.nvc | 13 +- src/xlqc-omp/README.md | 11 + src/xlqc-omp/main.cpp | 32 ++- src/xlqc-sycl/CMakeLists.txt | 4 +- src/xlqc-sycl/Makefile | 13 +- src/xlqc-sycl/README.md | 11 + src/xlqc-sycl/main.cpp | 35 ++-- 24 files changed, 487 insertions(+), 125 deletions(-) create mode 100644 src/xlqc-cuda/gsl_compat.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 0ef653d8cb..16688073f0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,6 +25,7 @@ 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..d2d35f9098 100644 --- a/README.md +++ b/README.md @@ -10,11 +10,12 @@ 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/), [GSL](https://www.gnu.org/software/gsl), [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
+GSL: sss
+Eigen: xlqc
GDAL: stsg
BZip2: gerbil diff --git a/cmake/modules/BenchmarkMacros.cmake b/cmake/modules/BenchmarkMacros.cmake index 4f6d18bb2d..491167e5ea 100644 --- a/cmake/modules/BenchmarkMacros.cmake +++ b/cmake/modules/BenchmarkMacros.cmake @@ -12,7 +12,10 @@ set(DEPEND_ON_BOOST "hbc" "ge-spmm" "mmcsf" "warpsort" "gerbil") 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") +set(DEPEND_ON_GSL "sss") + +# Global list for benchmarks that require Eigen +set(DEPEND_ON_EIGEN "xlqc") # Global list for benchmarks that require GDAL set(DEPEND_ON_GDAL "stsg") @@ -123,6 +126,14 @@ function(add_hecbench_benchmark) endif() endif() + 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() + if(${BENCH_NAME} IN_LIST DEPEND_ON_GDAL) if(NOT GDAL_FOUND) message(STATUS "Skipping ${BENCH_NAME}-${BENCH_MODEL_LOWER} (GDAL not found)") 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..655e30f642 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,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 -Xcompiler -Wall $(GSL_INC) -arch=$(ARCH) +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Xcompiler -Wall $(EIGEN_INC) -arch=$(ARCH) # 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..fc9928c055 --- /dev/null +++ b/src/xlqc-cuda/gsl_compat.h @@ -0,0 +1,300 @@ +/***************************************************************************** + 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 + +#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..e5ddf8d181 100644 --- a/src/xlqc-cuda/main.cu +++ b/src/xlqc-cuda/main.cu @@ -26,11 +26,7 @@ #include #include -#include -#include -#include -#include -#include +#include "gsl_compat.h" #include "int_lib/cints.h" #include "int_lib/crys.h" @@ -590,29 +586,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); } //====== free device memories ======== + start = std::chrono::steady_clock::now(); cudaFree(dev_pbf_xlec); cudaFree(dev_pbf_to_cbf); 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..d3418ebd91 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) 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..282eb02312 100644 --- a/src/xlqc-hip/main.cu +++ b/src/xlqc-hip/main.cu @@ -26,11 +26,7 @@ #include #include -#include -#include -#include -#include -#include +#include "gsl_compat.h" #include "int_lib/cints.h" #include "int_lib/crys.h" @@ -586,29 +582,25 @@ 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); - } - - //====== free device memories ======== + start = std::chrono::steady_clock::now(); hipFree(dev_pbf_xlec); hipFree(dev_pbf_to_cbf); 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..9e5855f2ae 100644 --- a/src/xlqc-omp/main.cpp +++ b/src/xlqc-omp/main.cpp @@ -26,11 +26,7 @@ #include #include -#include -#include -#include -#include -#include +#include "gsl_compat.h" #include @@ -571,27 +567,25 @@ 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); - } - } // omp target //====== free host memories ======== 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..e20161a577 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 \ 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..8c8e14b372 100644 --- a/src/xlqc-sycl/main.cpp +++ b/src/xlqc-sycl/main.cpp @@ -26,11 +26,7 @@ of this software, even if advised of the possibility of such damage. #include #include -#include -#include -#include -#include -#include +#include "gsl_compat.h" #include "int_lib/cints.h" #include "int_lib/crys.h" @@ -638,28 +634,25 @@ 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); } //====== free device memories ======== + start = std::chrono::steady_clock::now(); sycl::free(d_pbf_xlec, q); sycl::free(d_pbf_to_cbf, q); From 2a98a7d61e1389e67c32ab1dacb7a68b5f910ecb Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Tue, 11 Aug 2026 19:10:43 -0700 Subject: [PATCH 2/4] [sss] replace GSL with a self-contained header, share sources, and fix kernel data races GSL is GPL, so distributing sss binaries built from CUDA/HIP sources conflicts with the CUDA EULA (https://github.com/ORNL/HeCBench/issues/319). sss only used GSL for a Mersenne Twister stream and one semi-infinite integral, so replace it with a self-contained gsl_compat.h that reproduces mt19937 with GSL's 2002 Knuth seeding and implements adaptive Gauss-Kronrod (QK15) quadrature from the public domain QUADPACK algorithm. The RNG stream is bit-identical to GSL 2.7.1 and the integral agrees to ~1e-12, so results are unchanged. sss now builds with no external dependency, so drop the GSL discovery and the guards that skipped the benchmark when GSL was absent. Deduplicate the variants: the files that carry no GPU API calls now live only in sss-cuda and are pulled in through an include path, removing about 4500 lines of copies. Only main, kernels, and DPmixGGM_SSSmoves stay per-variant, and sss-cuda/kernels.cu is now shared with sss-hip since the two were identical and it uses no CUDA-only names. Add the missing barriers in CanDeleteEdge and CanAddEdge. Both kernels let every thread write shared state and then had thread 0 read all of it without a barrier, and CanDeleteEdge tested a thread-0-only counter for its loop exit, so threads could leave the loop divergently. This was benign on NVIDIA, where a 32-thread block is one warp in lockstep, but on gfx908 thread 0 read stale values: sss-hip scored -5721.97 and reported num_allModels=2731275 against -5675.31 and 12504300 elsewhere, converging on a worse optimum after 4.6x fewer iterations and so also reporting a misleadingly short run time. All three variants now reproduce the CPU reference path iteration for iteration and write identical output. Also fix two sss-sycl reporting bugs: the stopwatch was restarted on every search restart, so wall_time only covered the final restart (0.0036s instead of 35s), and the initial progress line printed k before it was assigned. Co-authored-by: Cursor --- CMakeLists.txt | 1 - README.md | 3 +- cmake/modules/BenchmarkMacros.cmake | 10 - src/sss-cuda/CMakeLists.txt | 2 - src/sss-cuda/DPmixGGM.cpp | 1 + src/sss-cuda/Makefile | 10 +- src/sss-cuda/README.md | 36 +- src/sss-cuda/gsl_compat.h | 872 +++++++++++++++++++++++++ src/sss-cuda/gwish.cpp | 1 + src/sss-cuda/kernels.cu | 14 + src/sss-cuda/main.cu | 4 +- src/sss-hip/CMakeLists.txt | 3 +- src/sss-hip/DPmixGGM.cpp | 515 --------------- src/sss-hip/DPmixGGM_Lists.cpp | 413 ------------ src/sss-hip/Makefile | 15 +- src/sss-hip/README.md | 36 +- src/sss-hip/graph.cpp | 973 ---------------------------- src/sss-hip/graph.h | 117 ---- src/sss-hip/gwish.cpp | 167 ----- src/sss-hip/kernels.cu | 272 -------- src/sss-hip/main.cu | 4 +- src/sss-hip/utilities.cpp | 74 --- src/sss-sycl/CMakeLists.txt | 2 +- src/sss-sycl/DPmixGGM.cpp | 516 --------------- src/sss-sycl/DPmixGGM_Lists.cpp | 413 ------------ src/sss-sycl/Makefile | 16 +- src/sss-sycl/README.md | 36 +- src/sss-sycl/graph.cpp | 973 ---------------------------- src/sss-sycl/graph.h | 117 ---- src/sss-sycl/gwish.cpp | 168 ----- src/sss-sycl/kernels.cpp | 14 + src/sss-sycl/main.cpp | 7 +- src/sss-sycl/utilities.cpp | 74 --- 33 files changed, 994 insertions(+), 4885 deletions(-) create mode 100644 src/sss-cuda/gsl_compat.h delete mode 100644 src/sss-hip/DPmixGGM.cpp delete mode 100644 src/sss-hip/DPmixGGM_Lists.cpp delete mode 100644 src/sss-hip/graph.cpp delete mode 100644 src/sss-hip/graph.h delete mode 100644 src/sss-hip/gwish.cpp delete mode 100644 src/sss-hip/kernels.cu delete mode 100644 src/sss-hip/utilities.cpp delete mode 100644 src/sss-sycl/DPmixGGM.cpp delete mode 100644 src/sss-sycl/DPmixGGM_Lists.cpp delete mode 100644 src/sss-sycl/graph.cpp delete mode 100644 src/sss-sycl/graph.h delete mode 100644 src/sss-sycl/gwish.cpp delete mode 100644 src/sss-sycl/utilities.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 16688073f0..bcea2f625c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,7 +24,6 @@ set(CMAKE_CXX_EXTENSIONS OFF) find_package(Boost) find_package(MPI) find_package(GDAL CONFIG) -find_package(GSL) find_package(Eigen3) find_package(BZip2) diff --git a/README.md b/README.md index d2d35f9098..85472bc473 100644 --- a/README.md +++ b/README.md @@ -10,11 +10,10 @@ 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), [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).
+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
Eigen: xlqc
GDAL: stsg
BZip2: gerbil diff --git a/cmake/modules/BenchmarkMacros.cmake b/cmake/modules/BenchmarkMacros.cmake index 491167e5ea..63528e0b6e 100644 --- a/cmake/modules/BenchmarkMacros.cmake +++ b/cmake/modules/BenchmarkMacros.cmake @@ -11,9 +11,6 @@ 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") - # Global list for benchmarks that require Eigen set(DEPEND_ON_EIGEN "xlqc") @@ -119,13 +116,6 @@ 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)") - return() - endif() - endif() - if(${BENCH_NAME} IN_LIST DEPEND_ON_EIGEN) # Eigen is header-only; the imported target carries the include path if(NOT TARGET Eigen3::Eigen) 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..7b7ca4780e --- /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--; + } + + size_t top; + if (last < (limit / 2 + 2)) + top = last; + else + top = limit - last + 1; + + size_t i = 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]; + + size_t 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 From fc44fe218109328dfeabb221d6ff57035a761846 Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Thu, 13 Aug 2026 08:03:24 -0700 Subject: [PATCH 3/4] [xlqc] double check the source and build files --- src/xlqc-cuda/Makefile | 3 ++- src/xlqc-cuda/gsl_compat.h | 1 + src/xlqc-cuda/main.cu | 7 ++++--- src/xlqc-hip/Makefile | 8 ++++---- src/xlqc-hip/main.cu | 7 ++++--- src/xlqc-omp/main.cpp | 10 +++++----- src/xlqc-sycl/Makefile | 2 ++ src/xlqc-sycl/main.cpp | 8 +++++--- 8 files changed, 27 insertions(+), 19 deletions(-) diff --git a/src/xlqc-cuda/Makefile b/src/xlqc-cuda/Makefile index 655e30f642..73fff78817 100644 --- a/src/xlqc-cuda/Makefile +++ b/src/xlqc-cuda/Makefile @@ -30,7 +30,8 @@ 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 $(EIGEN_INC) -arch=$(ARCH) +CFLAGS := $(EXTRA_CFLAGS) -std=c++17 -Xcompiler -Wall $(EIGEN_INC) -arch=$(ARCH) \ + --expt-relaxed-constexpr # Linker Flags LDFLAGS = diff --git a/src/xlqc-cuda/gsl_compat.h b/src/xlqc-cuda/gsl_compat.h index fc9928c055..5245c6509a 100644 --- a/src/xlqc-cuda/gsl_compat.h +++ b/src/xlqc-cuda/gsl_compat.h @@ -21,6 +21,7 @@ #define GSL_COMPAT_H #include +#include #include #include #include diff --git a/src/xlqc-cuda/main.cu b/src/xlqc-cuda/main.cu index e5ddf8d181..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 @@ -44,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) { @@ -601,6 +601,7 @@ int main(int argc, char* argv[]) } else { fprintf(stderr, "FAIL: E_total = %.10f, expected %.10f (error %.2e > tol %.0e)\n", ene_total, ref_energy, err, tol); + status = 1; } @@ -709,5 +710,5 @@ int main(int argc, char* argv[]) //====== the end of program ======== - return 0; + return status; } diff --git a/src/xlqc-hip/Makefile b/src/xlqc-hip/Makefile index d3418ebd91..b6dce1e676 100644 --- a/src/xlqc-hip/Makefile +++ b/src/xlqc-hip/Makefile @@ -57,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/main.cu b/src/xlqc-hip/main.cu index 282eb02312..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 @@ -44,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) { @@ -597,6 +597,7 @@ int main(int argc, char* argv[]) } else { fprintf(stderr, "FAIL: E_total = %.10f, expected %.10f (error %.2e > tol %.0e)\n", ene_total, ref_energy, err, tol); + status = 1; } //====== free device memories ======== @@ -704,5 +705,5 @@ int main(int argc, char* argv[]) //====== the end of program ======== - return 0; + return status; } diff --git a/src/xlqc-omp/main.cpp b/src/xlqc-omp/main.cpp index 9e5855f2ae..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,18 +17,16 @@ *****************************************************************************/ #include -#include #include #include #include #include #include #include +#include #include "gsl_compat.h" -#include - #include "int_lib/cints.h" #include "int_lib/crys.h" @@ -47,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) { @@ -582,6 +581,7 @@ int main(int argc, char* argv[]) } else { fprintf(stderr, "FAIL: E_total = %.10f, expected %.10f (error %.2e > tol %.0e)\n", ene_total, ref_energy, err, tol); + status = 1; } start = std::chrono::steady_clock::now(); @@ -682,5 +682,5 @@ int main(int argc, char* argv[]) //====== the end of program ======== - return 0; + return status; } diff --git a/src/xlqc-sycl/Makefile b/src/xlqc-sycl/Makefile index e20161a577..2dd36f61f6 100644 --- a/src/xlqc-sycl/Makefile +++ b/src/xlqc-sycl/Makefile @@ -97,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/main.cpp b/src/xlqc-sycl/main.cpp index 8c8e14b372..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,6 +24,7 @@ of this software, even if advised of the possibility of such damage. #include #include +#include #include "gsl_compat.h" @@ -35,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" @@ -46,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) { @@ -649,6 +650,7 @@ int main(int argc, char* argv[]) } else { fprintf(stderr, "FAIL: E_total = %.10f, expected %.10f (error %.2e > tol %.0e)\n", ene_total, ref_energy, err, tol); + status = 1; } //====== free device memories ======== @@ -755,5 +757,5 @@ int main(int argc, char* argv[]) //====== the end of program ======== - return 0; + return status; } From d2fc4a2b04c5a3cad4e948590c37e23fbfaf8b62 Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Thu, 13 Aug 2026 08:09:37 -0700 Subject: [PATCH 4/4] [sss] double check gsl_compat.h --- src/sss-cuda/gsl_compat.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/sss-cuda/gsl_compat.h b/src/sss-cuda/gsl_compat.h index 7b7ca4780e..4369bd91ae 100644 --- a/src/sss-cuda/gsl_compat.h +++ b/src/sss-cuda/gsl_compat.h @@ -362,13 +362,13 @@ inline void gsl_compat_qpsrt(gsl_integration_workspace *w) { i_nrmax--; } - size_t top; + int top; if (last < (limit / 2 + 2)) - top = last; + top = (int)last; else - top = limit - last + 1; + top = (int)(limit - last + 1); - size_t i = i_nrmax + 1; + int i = (int)i_nrmax + 1; while (i < top && errmax < elist[order[i]]) { order[i - 1] = order[i]; @@ -379,7 +379,7 @@ inline void gsl_compat_qpsrt(gsl_integration_workspace *w) { const double errmin = elist[last]; - size_t k = top - 1; + int k = top - 1; while (k > i - 2 && errmin >= elist[order[k]]) { order[k + 1] = order[k];