diff --git a/.gitignore b/.gitignore
index 48e421e..86f9481 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@
 CMakeFiles
 cmake-build*
 *.swp
+main
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0f7d6f5..ab49ecb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,9 @@
 cmake_minimum_required(VERSION 2.8)
 project(hpc)
 
-set(CMAKE_CXX_FLAGS "-std=gnu++11 -fopenmp -O3")
+add_subdirectory(common)
 
-add_executable(hpc main.cpp)
+add_subdirectory(cpu-slo)
+add_subdirectory(cpu-opt)
+
+add_subdirectory(gpu-slo)
diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt
new file mode 100644
index 0000000..7cdff4e
--- /dev/null
+++ b/common/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_library(common INTERFACE)
+target_include_directories(common INTERFACE include)
diff --git a/common/include/groups.hpp b/common/include/groups.hpp
new file mode 100644
index 0000000..e8cb22f
--- /dev/null
+++ b/common/include/groups.hpp
@@ -0,0 +1,110 @@
+#pragma once
+
+#include <types.hpp>
+
+struct Mult {
+    int from, to, multiplicity;
+};
+
+Table mults(const std::vector<Mult>& ms) {
+    Table res;
+    for (const auto &m : ms) {
+        int N = res.size();
+        res.emplace_back(m.multiplicity * 2, m.to);
+        for (int i = 0; i < m.multiplicity * 2; i += 2) {
+            res[N][i] = m.from;
+        }
+    }
+    return res;
+}
+
+std::vector<Mult> ezmults(int ngens, const std::vector<Mult> &ms) {
+    bool table[ngens][ngens];
+
+    for (int i = 0; i < ngens; i++) {
+        for (int j = 0; j < ngens; j++) {
+            table[i][j] = false;
+        }
+    }
+
+    for (const auto &m : ms) {
+        table[m.from][m.to] = true;
+        table[m.to][m.from] = true;
+    }
+
+    std::vector<Mult> res(ms);
+
+    for (int i = 0; i < ngens; i++) {
+        for (int j = i + 1; j < ngens; j++) {
+            if (!table[i][j]) {
+                res.push_back({i, j, 2});
+            }
+        }
+    }
+
+    return res;
+}
+
+/*
+ * Order 4*res*res
+ */
+std::pair<Table, int> torus(int res) {
+
+    return std::make_pair(mults(ezmults(4, {
+        {0, 1, res},
+        {2, 3, res},
+    })), 4);
+}
+
+/*
+ * Order 14,400
+ */
+std::pair<Table, int> H4() {
+    return std::make_pair(mults(ezmults(4, {
+        {0, 1, 5},
+        {1, 2, 3},
+        {2, 3, 3},
+    })), 4);
+}
+
+/*
+ * Order 51,840
+ */
+std::pair<Table, int> E6() {
+    return std::make_pair(mults(ezmults(6, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+    })), 6);
+}
+
+/*
+ * Order 2,903,040
+ */
+std::pair<Table, int> E7() {
+    return std::make_pair(mults(ezmults(7, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+    })), 7);
+}
+
+/*
+ * Order 696,729,600
+ */
+std::pair<Table, int> E8() {
+    return std::make_pair(mults(ezmults(8, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+        {6, 7, 3},
+    })), 8);
+}
diff --git a/common/include/types.hpp b/common/include/types.hpp
new file mode 100644
index 0000000..18b23be
--- /dev/null
+++ b/common/include/types.hpp
@@ -0,0 +1,4 @@
+#pragma once
+
+using Gens=std::vector<int>;
+using Table=std::vector<Gens>;
diff --git a/cpu-opt/CMakeLists.txt b/cpu-opt/CMakeLists.txt
new file mode 100644
index 0000000..4d7e286
--- /dev/null
+++ b/cpu-opt/CMakeLists.txt
@@ -0,0 +1,3 @@
+set(CMAKE_CXX_FLAGS "-std=gnu++11 -fopenmp -O3")
+
+add_executable(cpu-opt main.cpp solver.hpp)
diff --git a/cpu-opt/main.cpp b/cpu-opt/main.cpp
new file mode 100644
index 0000000..a55e4ca
--- /dev/null
+++ b/cpu-opt/main.cpp
@@ -0,0 +1,440 @@
+#include <cstdlib>
+#include <immintrin.h>
+
+#include <iostream>
+#include <vector>
+#include <iomanip>
+#include <chrono>
+
+#include <omp.h>
+
+using Ind=int;
+
+using Gen=int;
+using Gens=std::vector<Gen>;
+using Table=std::vector<Gens>;
+
+using Cos=int;
+
+/*
+ * COXETER GROUP DEFINITIONS
+ */
+
+const size_t ALIGN_SIZE=64;
+const size_t GENS_PER_LINE=ALIGN_SIZE/sizeof(Gen);
+
+struct Coxeter {
+    Gen *gen[2];
+    int *size; // multiplicity * 2
+    int ngens;
+    int nrels;
+    Coxeter(int ngens): ngens(ngens){
+        nrels = (ngens*(ngens-1))>>1;
+        int allocsize = ((nrels-1)/(GENS_PER_LINE) + 1)*ALIGN_SIZE;
+        int err = 0;
+        err |= posix_memalign((void**)&(gen[0]),ALIGN_SIZE,allocsize);
+        err |= posix_memalign((void**)&(gen[1]),ALIGN_SIZE,allocsize);
+        err |= posix_memalign((void**)&(size),ALIGN_SIZE,allocsize);
+        if (err != 0) {
+            std::cerr << "Error allocating Coxeter!" << std::endl;
+            exit(EXIT_FAILURE);
+        }
+    }
+    void clean() {
+        free(gen[0]);
+        free(gen[1]);
+        free(size);
+    }
+};
+
+struct Mult {
+    int from, to, multiplicity;
+    Mult() {}
+    Mult(int from, int to, int multiplicity): from(from), to(to), multiplicity(multiplicity) {}
+};
+
+Coxeter make_coxeter(int ngens, const std::vector<Mult> &ms) {
+    int mults[ngens][ngens];
+
+    for (int i = 0; i < ngens; i++) {
+        for (int j = 0; j < ngens; j++) {
+            mults[i][j] = 2;
+            mults[j][i] = 2;
+        }
+    }
+
+    for (const auto &m : ms) {
+        mults[m.from][m.to] = m.multiplicity;
+        mults[m.to][m.from] = m.multiplicity;
+    }
+
+    Coxeter c(ngens);
+    int k=0;
+    for (int i = 0; i < ngens; i++) {
+        for (int j = i + 1; j < ngens; j++) {
+            int size = mults[i][j]<<1;
+            c.gen[0][k] = i;
+            c.gen[1][k] = j;
+            c.size[k] = size;
+            k++;
+        }
+    }
+
+    return c;
+}
+
+/*
+ * Order 4*res*res
+ */
+Coxeter torus(int res) {
+    return make_coxeter(4, {
+        {0, 1, res},
+        {2, 3, res},
+    });
+}
+
+Coxeter hypercube(int dim) {
+    std::vector<Mult> mults;
+    mults.emplace_back(0,1,4);
+    for (int i = 2; i < dim; i++) {
+        mults.emplace_back(i-1, i, 3);
+    }
+    return make_coxeter(dim, mults);
+}
+
+/*
+ * Order 14,400
+ */
+Coxeter H4() {
+    return make_coxeter(4, {
+        {0, 1, 5},
+        {1, 2, 3},
+        {2, 3, 3},
+    });
+}
+
+/*
+ * Order 51,840
+ */
+Coxeter E6() {
+    return make_coxeter(6, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+    });
+}
+
+/*
+ * Order 2,903,040
+ */
+Coxeter E7() {
+    return make_coxeter(7, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+    });
+}
+
+/*
+ * Order 696,729,600
+ */
+Coxeter E8() {
+    return make_coxeter(8, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+        {6, 7, 3},
+    });
+}
+
+/*
+ * LEARNING / RelTable DIFINITIONS
+ */
+
+struct RelTable {
+    std::vector<int> coset_poss;
+    std::vector<Cos> init_cosets;
+    std::vector<Cos> start_cosets;
+    std::vector<Cos> end_cosets;
+    std::vector<Ind> start_inds;
+    std::vector<Ind> end_inds;
+    int num_rows;
+    Gen gen[2];
+    Ind end_ind;
+    RelTable(Gen gen0, Gen gen1, Ind end_ind): end_ind(end_ind), num_rows(0) {
+        gen[0] = gen0;
+        gen[1] = gen1;
+    }
+    void add_row(Cos new_coset) {
+        coset_poss.push_back(num_rows);
+        init_cosets.push_back(new_coset);
+        start_cosets.push_back(new_coset);
+        end_cosets.push_back(new_coset);
+        start_inds.push_back(0);
+        end_inds.push_back(end_ind);
+        num_rows++;
+    }
+    void rem_row(int idx) {
+        num_rows--;
+
+        coset_poss[init_cosets[num_rows]] = idx;
+        coset_poss[init_cosets[idx]] = -1;
+
+        init_cosets[idx] = init_cosets[num_rows];
+        init_cosets.pop_back();
+
+        start_cosets[idx] = start_cosets[num_rows];
+        start_cosets.pop_back();
+
+        end_cosets[idx] = end_cosets[num_rows];
+        end_cosets.pop_back();
+
+        start_inds[idx] = start_inds[num_rows];
+        start_inds.pop_back();
+
+        end_inds[idx] = end_inds[num_rows];
+        end_inds.pop_back();
+    }
+};
+
+struct CosetTable {
+    std::vector<Cos> table;
+    int num_cosets;
+    int ngens;
+    CosetTable (int ngens): ngens(ngens), num_cosets(0) {}
+    void add_row() {
+        num_cosets++;
+        table.resize(table.size()+ngens, -1);
+    }
+    inline Cos &operator[](int idx) {
+        return table[idx];
+    }
+};
+
+void pp(const Gens &g, int w) {
+    for (const auto &e : g) {
+        std::cerr << std::setw(w) << e << " ";
+    }
+    std::cerr << std::endl;
+}
+
+void pp(const Table &t) {
+    std::cerr << "| table:" << std::endl;
+    int w = 3;
+
+    for (size_t i = 0; i < t.size(); ++i) {
+        std::cerr << std::setw(w) << i << " | ";
+        pp(t[i], w);
+    }
+}
+
+void add_row(const Coxeter &cox,
+    CosetTable &cosets, std::vector<RelTable> &reltables) {
+
+    int C = cosets.num_cosets;
+
+    cosets.add_row();
+
+    for (RelTable &rt : reltables)
+        rt.add_row(C);
+}
+
+int add_coset(const Coxeter &cox,
+    CosetTable &cosets, std::vector<RelTable> &reltables,
+    int coset_scan_hint) {
+
+    const int C = cosets.num_cosets;
+    const int ngens = cox.ngens;
+    
+    int idx = coset_scan_hint*ngens;
+    for (int c = coset_scan_hint; c < C; ++c) {
+        for (int g = 0; g < ngens; ++g, ++idx) {
+            if (cosets[idx] == -1) {
+                cosets[idx] = C;
+                add_row(cox, cosets, reltables);
+                cosets[C*ngens+g] = c;
+                return c;
+            }
+        }
+    }
+
+    return -1;
+}
+
+
+/**
+ * learn until it can't
+ */
+inline void learn(const Coxeter &cox, CosetTable &cosets,
+    std::vector<RelTable> &reltables) {
+
+    const int nrels = cox.nrels;
+    const int ngens = cox.ngens;
+
+    Gen gens[2];
+    Ind s_i, e_i;
+    Cos s_c, e_c, i_c;
+    int lookup, idx;
+    Gen g;
+#pragma omp parallel for schedule(dynamic, 1) private(gens, s_i, e_i, s_c, e_c, i_c, lookup, idx, g)
+    for (unsigned int r = 0; r < nrels; ++r) { 
+        auto &table = reltables[r];
+        auto &start_inds = table.start_inds;
+        auto &end_inds = table.end_inds;
+        auto &start_cosets = table.start_cosets;
+        auto &end_cosets = table.end_cosets;
+        auto &init_cosets = table.init_cosets;
+        auto &coset_poss = table.coset_poss;
+        gens[0] = table.gen[0];
+        gens[1] = table.gen[1];
+
+        for (unsigned int c = 0; c < table.num_rows; c++) {
+            s_i = start_inds[c];
+            e_i = end_inds[c];
+            s_c = start_cosets[c];
+            e_c = end_cosets[c];
+            i_c = init_cosets[c];
+
+            while (s_i < e_i) {
+                lookup = cosets[gens[s_i&1] + s_c*ngens];
+                if (lookup < 0) break;
+
+                s_i++;
+                s_c = lookup;
+
+                if (s_c > i_c) {
+                    idx = coset_poss[s_c];
+                    if (idx >= 0)
+                        table.rem_row(idx); 
+                }
+            }
+
+            start_inds[c] = s_i;
+            start_cosets[c] = s_c;
+
+            while (s_i < e_i) {
+                lookup = cosets[gens[e_i&1] + e_c*ngens];
+                if (lookup < 0) break;
+
+                e_i--;
+                e_c = lookup;
+
+                if (e_c > i_c) {
+                    idx = coset_poss[e_c];
+                    if (idx >= 0)
+                        table.rem_row(idx);
+                }
+            }
+
+            end_inds[c] = e_i;
+            end_cosets[c] = e_c;
+
+            if (s_i == e_i) {
+                g = gens[s_i&1];
+                cosets[s_c*ngens + g] = e_c;
+                cosets[e_c*ngens + g] = s_c;
+
+                table.rem_row(c);
+                c--;
+            }
+        }
+    }
+}
+
+CosetTable solve_tc(const Coxeter &cox, const Gens &subgens) {
+    CosetTable cosets(cox.ngens);
+    std::vector<RelTable> reltables;
+
+    for (int i=0; i<cox.nrels; i++) {
+        reltables.emplace_back(cox.gen[0][i], cox.gen[1][i], cox.size[i]-1);
+    }
+
+    // set up initial coset
+    add_row(cox, cosets, reltables);
+    for (const auto &gen : subgens) {
+        cosets[gen] = 0;
+    }
+
+    int coset_scan_hint = 0;
+    char a;
+    while (coset_scan_hint >= 0) {
+        learn(cox, cosets, reltables);
+        coset_scan_hint = add_coset(cox, cosets, reltables, coset_scan_hint);
+    }
+
+    return cosets;
+}
+
+/*
+ * returns coxeter group based on the arguments
+ * prints out type and arguments, without an endline
+ */
+Coxeter proc_args(int argc, const char* argv[]) {
+    if (argc < 2) {
+        std::cerr << "missing type argument." << std::endl;
+        exit(EXIT_FAILURE);
+    }
+
+    int type = std::strtol(argv[1], nullptr, 10);
+    std::cout << type << ',';
+
+    int arg;
+    switch (type) {
+    case 0:
+        if (argc < 3) {
+            std::cerr << "Must provide a size for torus!" << std::endl;
+            exit(EXIT_FAILURE);
+        }
+        arg = std::strtol(argv[2], nullptr, 10);
+        std::cout << arg << ',';
+        return torus(arg);
+    case 1:
+        std::cout << -1 << ',';
+        return H4();
+    case 2:
+        std::cout << -1 << ',';
+        return E6();
+    case 3:
+        std::cout << -1 << ',';
+        return E7();
+    case 4:
+        std::cout << -1 << ',';
+        return E8();
+    case 5:
+        if (argc < 3) {
+            std::cerr << "Must provide a dimension for hypercube!" << std::endl;
+            exit(EXIT_FAILURE);
+        }
+        arg = std::strtol(argv[2], nullptr, 10);
+        std::cout << arg << ',';
+        return hypercube(arg);
+    }
+
+    std::cerr << "Not a valid type!" << std::endl;
+    exit(EXIT_FAILURE);
+}
+
+int main(int argc, const char *argv[]) {
+    Coxeter cox = proc_args(argc, argv);
+
+    auto s = std::chrono::system_clock::now();
+    auto cosets = solve_tc(cox, {});
+    auto e = std::chrono::system_clock::now();
+
+    std::chrono::duration<float> diff = e - s;
+    size_t order = cosets.num_cosets;
+
+    // type,arg,ngens,time,order
+    std::cout << cox.ngens << ',' << diff.count() << ',' << order << std::endl;
+
+    cox.clean();
+    return 0;
+}
diff --git a/cpu-opt/solver.hpp b/cpu-opt/solver.hpp
new file mode 100644
index 0000000..79469d3
--- /dev/null
+++ b/cpu-opt/solver.hpp
@@ -0,0 +1,195 @@
+#pragma once
+
+#include <map>
+#include <vector>
+#include <algorithm>
+
+struct Mults;
+struct Table;
+struct IRow;
+
+struct Mults {
+    const int dim;
+    std::map<std::tuple<int, int>, int> mults;
+
+    explicit Mults(const int dim) : dim(dim) {
+        mults = std::map<std::tuple<int, int>, int>();
+    }
+
+    void set(int a, int b, int mult) {
+        if (a > b) std::swap(a, b);
+        mults[std::make_tuple(a, b)] = mult;
+    }
+
+    [[nodiscard]] int get(int a, int b) const {
+        if (a == b) return 1;
+        if (a > b) std::swap(a, b);
+
+        const auto &tup = std::make_tuple(a, b);
+        const auto &res = mults.find(tup);
+
+        if (res == mults.end()) return 2;
+        return res->second;
+    }
+
+    [[nodiscard]] std::vector<int> irelation(int a, int b) const {
+        std::vector<int> res{};
+        int mult = get(a, b);
+
+        for (int i = 0; i < mult; ++i) {
+            res.push_back(a);
+            res.push_back(b);
+        }
+
+        return res;
+    }
+
+    [[nodiscard]] Mults sub(const std::vector<int> &igens) const {
+        Mults res(igens.size());
+        for (int a = 0; a < (int) igens.size(); ++a) {
+            for (int b = a + 1; b < (int) igens.size(); ++b) {
+                res.set(a, b, get(igens[a], igens[b]));
+            }
+        }
+        return res;
+    }
+
+    [[nodiscard]] Table *isolve(const std::vector<int> &isubgens) const;
+};
+
+struct Table {
+    const Mults mults;
+    std::vector<std::vector<int>> fwd;
+
+    explicit Table(Mults mults) : mults(mults) {
+        add_row();
+    }
+
+    void add_row() {
+        fwd.emplace_back(mults.dim, -1);
+    }
+
+    int add_coset() {
+        for (int from = 0; from < (int) size(); ++from) {
+            for (int igen = 0; igen < (int) mults.dim; igen++) {
+                if (iget(from, igen) < 0) {
+                    int to = (int) size();
+                    add_row();
+                    iset(from, igen, to);
+                    return to;
+                }
+            }
+        }
+        return 0;
+    }
+
+    [[nodiscard]] unsigned size() const {
+        return fwd.size();
+    }
+
+    void iset(int from, int igen, int to) {
+        fwd[from][igen] = to;
+        fwd[to][igen] = from;
+    }
+
+    [[nodiscard]] int iget(int from, int igen) const {
+        return fwd[from][igen];
+    }
+
+    [[nodiscard]] int irget(int igen, int to) const {
+        return fwd[to][igen];
+    }
+
+};
+
+struct IRow {
+    std::vector<int>::const_iterator l;
+    std::vector<int>::const_iterator r;
+
+    int from;
+    int to;
+
+    IRow(const std::vector<int> &rel, int cos)
+        : l(rel.begin()), r(rel.end() - 1), from(cos), to(cos) {
+    }
+
+    [[nodiscard]] bool learn(Table *table) {
+        if (r - l == 0) {
+            return false;
+        }
+
+        while (r - l > 0) {
+            int next = table->iget(from, *l);
+            if (next < 0) break;
+            l++;
+            from = next;
+        }
+
+        while (r - l > 0) {
+            int next = table->irget(*r, to);
+            if (next < 0) break;
+            r--;
+            to = next;
+        }
+
+        if (r - l == 0) {
+            table->iset(from, *l, to);
+            return true;
+        }
+
+        return false;
+    }
+};
+
+Mults schlafli(const std::vector<int> &symbol) {
+    unsigned int dim = symbol.size();
+    Mults mults(dim + 1);
+    for (int i = 0; i < dim; ++i) {
+        mults.set(i, i + 1, symbol[i]);
+    }
+    return mults;
+}
+
+Table *Mults::isolve(const std::vector<int> &isubgens) const {
+    auto *table = new Table(*this);
+    for (int igen : isubgens)
+        table->iset(0, igen, 0);
+
+    std::vector<std::vector<int>> irels{};
+    for (unsigned a = 0; a < dim; ++a) {
+        for (unsigned b = a + 1; b < dim; ++b) {
+            irels.push_back(irelation(a, b));
+        }
+    }
+
+    std::vector<IRow> irows;
+    irows.reserve(irels.size());
+    for (const auto &irel : irels)
+        irows.emplace_back(irel, 0);
+
+    while (!irows.empty()) {
+        while (true) {
+            bool learned = false;
+            for (int i = (int) irows.size() - 1; i >= 0; i--) {
+                if (irows[i].learn(table)) {
+                    learned = true;
+                    irows.erase(irows.begin() + i);
+                }
+            }
+            if (!learned)
+                break;
+        }
+
+        int i = (int) table->size();
+        if (table->add_coset() > 0) {
+            for (const auto &irel : irels)
+                irows.emplace_back(irel, i);
+        } else {
+            break;
+        }
+    }
+
+    return table;
+}
+
+
diff --git a/cpu-slo/CMakeLists.txt b/cpu-slo/CMakeLists.txt
new file mode 100644
index 0000000..2b63699
--- /dev/null
+++ b/cpu-slo/CMakeLists.txt
@@ -0,0 +1,4 @@
+set(CMAKE_CXX_FLAGS "-std=gnu++11 -fopenmp -O3")
+
+add_executable(cpu-slo main.cpp)
+target_link_libraries(cpu-slo PRIVATE common)
diff --git a/Makefile b/cpu-slo/Makefile
similarity index 100%
rename from Makefile
rename to cpu-slo/Makefile
diff --git a/main.cpp b/cpu-slo/main.cpp
similarity index 62%
rename from main.cpp
rename to cpu-slo/main.cpp
index 52ada27..7b798a5 100644
--- a/main.cpp
+++ b/cpu-slo/main.cpp
@@ -8,6 +8,123 @@
 using Gens=std::vector<int>;
 using Table=std::vector<Gens>;
 
+struct Mult {
+    int from, to, multiplicity;
+    Mult() {}
+    Mult(int from, int to, int multiplicity): from(from), to(to), multiplicity(multiplicity) {}
+};
+
+Table mults(const std::vector<Mult>& ms) {
+    Table res;
+    for (const auto &m : ms) {
+        int N = res.size();
+        res.emplace_back(m.multiplicity * 2, m.to);
+        for (int i = 0; i < m.multiplicity * 2; i += 2) {
+            res[N][i] = m.from;
+        }
+    }
+    return res;
+}
+
+std::vector<Mult> ezmults(int ngens, const std::vector<Mult> &ms) {
+    bool table[ngens][ngens];
+
+    for (int i = 0; i < ngens; i++) {
+        for (int j = 0; j < ngens; j++) {
+            table[i][j] = false;
+        }
+    }
+
+    for (const auto &m : ms) {
+        table[m.from][m.to] = true;
+        table[m.to][m.from] = true;
+    }
+
+    std::vector<Mult> res(ms);
+
+    for (int i = 0; i < ngens; i++) {
+        for (int j = i + 1; j < ngens; j++) {
+            if (!table[i][j]) {
+                res.push_back({i, j, 2});
+            }
+        }
+    }
+
+    return res;
+}
+
+/*
+ * Order 4*res*res
+ */
+std::pair<Table, int> torus(int res) {
+    return std::make_pair(mults(ezmults(4, {
+        {0, 1, res},
+        {2, 3, res},
+    })), 4);
+}
+
+std::pair<Table, int> hypercube(int dim) {
+    std::vector<Mult> hc_mults;
+    hc_mults.emplace_back(0,1,4);
+    for (int i = 2; i < dim; i++) {
+        hc_mults.emplace_back(i-1, i, 3);
+    }
+    return std::make_pair(mults(ezmults(dim,hc_mults)), dim); 
+}
+
+/*
+ * Order 14,400
+ */
+std::pair<Table, int> H4() {
+    return std::make_pair(mults(ezmults(4, {
+        {0, 1, 5},
+        {1, 2, 3},
+        {2, 3, 3},
+    })), 4);
+}
+
+/*
+ * Order 51,840
+ */
+std::pair<Table, int> E6() {
+    return std::make_pair(mults(ezmults(6, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+    })), 6);
+}
+
+/*
+ * Order 2,903,040
+ */
+std::pair<Table, int> E7() {
+    return std::make_pair(mults(ezmults(7, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+    })), 7);
+}
+
+/*
+ * Order 696,729,600
+ */
+std::pair<Table, int> E8() {
+    return std::make_pair(mults(ezmults(8, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+        {6, 7, 3},
+    })), 8);
+}
+
 void pp(const Gens &g, int w) {
     for (const auto &e : g) {
         std::cout << std::setw(w) << e << " ";
@@ -149,144 +266,47 @@ Table solve_tc(int ngens, const Gens &subgens, const std::vector<Gens> &rels) {
     return cosets;
 }
 
-struct Mult {
-    int from, to, multiplicity;
-};
-
-Table mults(std::vector<Mult> ms) {
-    Table res;
-    for (const auto &m : ms) {
-	     int N = res.size();
-        res.emplace_back(m.multiplicity * 2, m.to);
-        for (int i = 0; i < m.multiplicity * 2; i += 2) {
-            res[N][i] = m.from;
-        }
-    }
-    return res;
-}
-
-Table ezmults(int ngens, std::vector<Mult> ms) {
-   bool table[ngens][ngens];
-
-   for (int i=0; i<ngens; i++) {
-      for (int j=0; j<ngens; j++) {
-         table[i][j] = false;
-      }
-   }
-
-   for (const auto &m : ms) {
-      table[m.from][m.to] = true;
-      table[m.to][m.from] = true;
-   }
-
-   Table res = mults(ms);
-   for (int i=0; i<ngens; i++) {
-      for (int j=i+1; j<ngens; j++) {
-         if (!table[i][j]) {
-            res.push_back({i,j,i,j});
-         }
-      }
-   }
-   return res;
-}
-
-/*
- * Order 4*res*res
- */
-std::pair<Table,int> torus(int res) {
-    return std::make_pair(ezmults(4,{
-        {0,1,res},
-        {2,3,res},
-    }),4);
-}
-
-/*
- * Order 14,400
- */
-std::pair<Table,int> H4() {
-   return std::make_pair(ezmults(4,{
-      {0,1,5},
-      {1,2,3},
-      {2,3,3},
-   }),4);
-}
-
-/*
- * Order 51,840
- */
-std::pair<Table,int> E6() {
-   return std::make_pair(ezmults(6,{
-      {0,1,3},
-      {1,2,3},
-      {2,3,3},
-      {2,4,3},
-      {4,5,3},
-   }),6);
-}
-
-/*
- * Order 2,903,040
- */
-std::pair<Table,int> E7() {
-   return std::make_pair(ezmults(7,{
-      {0,1,3},
-      {1,2,3},
-      {2,3,3},
-      {2,4,3},
-      {4,5,3},
-      {5,6,3},
-   }),7);
-}
-
-/*
- * Order 696,729,600
- */
-std::pair<Table,int> E8() {
-   return std::make_pair(ezmults(8,{
-      {0,1,3},
-      {1,2,3},
-      {2,3,3},
-      {2,4,3},
-      {4,5,3},
-      {5,6,3},
-      {6,7,3},
-   }),8);
-}
-
-int main(int argc, char *argv[]) {
+int main(int argc, const char *argv[]) {
     if (argc < 2) {
-        std::cerr << "REEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE" << std::endl;
-        std::cerr << "gimme more ~~tendies~~ arguments" << std::endl;
+        std::cerr << "missing type argument." << std::endl;
         return 1;
     }
 
     int type = std::strtol(argv[1], nullptr, 10);
-    int arg = 0;
-    std::pair<Table,int> res;
+    int arg = -1;
+    std::pair<Table, int> res;
     switch (type) {
-      case 0:
-         if (argc < 3) {
+    case 0:
+        if (argc < 3) {
             std::cerr << "Must provide a size for torus!" << std::endl;
             return 2;
-         }
-         arg = std::strtol(argv[2], nullptr, 10);
-         res = torus(arg);
-         break;
-      case 1:
-         res = H4();
-         break;
-      case 2:
-         res = E6();
-         break;
-      case 3:
-         res = E7();
-         break;
-      case 4:
-         res = E8();
-         break;
-      default:
-         std::cerr << "Not a valid type!" << std::endl;
-         return 3;
+        }
+        arg = std::strtol(argv[2], nullptr, 10);
+        res = torus(arg);
+        break;
+    case 1:
+        res = H4();
+        break;
+    case 2:
+        res = E6();
+        break;
+    case 3:
+        res = E7();
+        break;
+    case 4:
+        res = E8();
+        break;
+    case 5:
+        if (argc < 3) {
+            std::cerr << "Must provide a dimension for hypercube!" << std::endl;
+            return 3;
+        }
+        arg = std::strtol(argv[2], nullptr, 10);
+        res = hypercube(arg);
+        break;
+    default:
+        std::cerr << "Not a valid type!" << std::endl;
+        return 3;
     }
 
     const Table &rels = res.first;
@@ -300,7 +320,7 @@ int main(int argc, char *argv[]) {
     size_t order = cosets.size();
 
     // ngens,type,arg,time,order
-    std::cout << ngens << ',' << type << ',' << arg  << ',' << diff.count() << ',' << order << std::endl;
+    std::cout << type << ',' << arg << ',' << ngens << ',' << diff.count() << ',' << order << std::endl;
 
     return 0;
 }
diff --git a/gpu-opt/CMakeLists.txt b/gpu-opt/CMakeLists.txt
new file mode 100644
index 0000000..6ae6437
--- /dev/null
+++ b/gpu-opt/CMakeLists.txt
@@ -0,0 +1,6 @@
+cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
+project(main LANGUAGES CXX CUDA)
+
+add_executable(main main.cu)
+target_compile_features(main PUBLIC cxx_std_11)
+set_target_properties(main PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
diff --git a/gpu-opt/Makefile b/gpu-opt/Makefile
new file mode 100644
index 0000000..0e5d935
--- /dev/null
+++ b/gpu-opt/Makefile
@@ -0,0 +1,9 @@
+main : main.cu util.h
+	nvcc -o main -std=c++11 -O3 main.cu
+
+clean :
+	rm main
+
+queue : main
+	qsub -q mamba -l walltime=24:00:00 -l nodes=1:ppn=1:gpus=1 -d `pwd` run_gpu_opt.sh
+
diff --git a/gpu-opt/groups.h b/gpu-opt/groups.h
new file mode 100644
index 0000000..8f7768a
--- /dev/null
+++ b/gpu-opt/groups.h
@@ -0,0 +1,127 @@
+#pragma once
+
+#include "util.h"
+
+#include <iostream>
+
+/*
+ * Order 4*res*res
+ */
+Coxeter torus(int res) {
+    return make_coxeter(4, {
+        {0, 1, res},
+        {2, 3, res},
+    });
+}
+
+Coxeter hypercube(int dim) {
+    std::vector<Rel> rels;
+    rels.push_back({0, 1, 4});
+    for (int d = 2; d < dim; d++) {
+        rels.push_back({d-1, d, 3});
+    }
+    return make_coxeter(dim, rels);
+}
+
+/*
+ * Order 14,400
+ */
+Coxeter H4() {
+    return make_coxeter(4, {
+        {0, 1, 5},
+        {1, 2, 3},
+        {2, 3, 3},
+    });
+}
+
+/*
+ * Order 51,840
+ */
+Coxeter E6() {
+    return make_coxeter(6, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+    });
+}
+
+/*
+ * Order 2,903,040
+ */
+Coxeter E7() {
+    return make_coxeter(7, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+    });
+}
+
+/*
+ * Order 696,729,600
+ */
+Coxeter E8() {
+    return make_coxeter(8, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+        {6, 7, 3},
+    });
+}
+
+/*
+ * returns coxeter group based on the arguments
+ * prints out type and arguments, without an endline
+ */
+Coxeter proc_args(int argc, const char* argv[]) {
+    if (argc < 2) {
+        std::cerr << "missing type argument." << std::endl;
+        exit(EXIT_FAILURE);
+    }
+
+    int type = std::strtol(argv[1], nullptr, 10);
+    std::cout << type << ',';
+
+    int arg;
+    switch (type) {
+    case 0:
+        if (argc < 3) {
+            std::cerr << "Must provide a size for torus!" << std::endl;
+            exit(EXIT_FAILURE);
+        }
+        arg = std::strtol(argv[2], nullptr, 10);
+        std::cout << arg << ',';
+        return torus(arg);
+    case 1:
+        std::cout << -1 << ',';
+        return H4();
+    case 2:
+        std::cout << -1 << ',';
+        return E6();
+    case 3:
+        std::cout << -1 << ',';
+        return E7();
+    case 4:
+        std::cout << -1 << ',';
+        return E8();
+    case 5:
+        if (argc < 3) {
+            std::cerr << "Must provide a dimension for hypercube!" << std::endl;
+            exit(EXIT_FAILURE);
+        }
+        arg = std::strtol(argv[2], nullptr, 10);
+        std::cout << arg << ',';
+        return hypercube(arg);
+    }
+
+    std::cerr << "Not a valid type!" << std::endl;
+    exit(EXIT_FAILURE);
+}
+
diff --git a/gpu-opt/main.cu b/gpu-opt/main.cu
new file mode 100644
index 0000000..df671d3
--- /dev/null
+++ b/gpu-opt/main.cu
@@ -0,0 +1,245 @@
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <thrust/logical.h>
+
+#include <vector>
+#include <iostream>
+#include <chrono>
+
+#include "util.h"
+#include "groups.h"
+
+__constant__ Rel c_rels[128];
+__constant__ int c_nrels[1];
+__constant__ int c_ngens[1];
+
+struct Row {
+    int rel;
+
+    int l, r;
+
+    int from, to;
+
+    __host__ __device__
+    Row() : rel(0), l(0), r(0), from(0), to(0) {}
+
+    __device__
+    Row(int rel, int cos, int size) {
+        l = 0;
+        r = size - 1;
+        
+        from = to = cos;
+        
+        this->rel = rel;
+    }
+};
+
+std::ostream &operator<<(std::ostream &o, const Row &r) {
+    return o << "Row[" << r.rel << "]{" << r.l << ":" << r.from << "-" << r.to << ":" << r.r << "}";
+}
+
+// this performs a pass on one relation table row, applying learned data to the coset table.
+struct Solver {
+    int *cosets;
+    
+    Solver(thrust::device_vector<int> &cosets)
+        : cosets(thrust::raw_pointer_cast(cosets.data())) {
+    }
+    
+    __device__
+    void operator()(Row &drow) {
+        Row row = drow;
+
+        if (row.r - row.l <= 0) {
+            return;
+        }
+        
+        while (row.r - row.l > 0) {
+            int gen = c_rels[row.rel].gens[row.l & 1];
+            int next = cosets[row.from * c_ngens[0] + gen];
+            if (next < 0) break;
+            row.l++;
+            row.from = next;
+        }
+
+        while (row.r - row.l > 0) {
+            int gen = c_rels[row.rel].gens[row.r & 1];
+            int next = cosets[row.to * c_ngens[0] + gen];
+            if (next < 0) break;
+            row.r--;
+            row.to = next;
+        }
+
+        drow = row;
+            
+        if (row.r - row.l <= 0) { 
+            int gen = c_rels[row.rel].gens[row.l & 1];
+            cosets[row.from * c_ngens[0] + gen] = row.to;
+            cosets[row.to * c_ngens[0] + gen] = row.from;
+            return;
+        }
+    }
+};
+
+// this sets the inital row in the coset table based on the subgroup generators
+struct CosetInitializer {
+    int *cosets;
+
+    CosetInitializer(thrust::device_vector<int> &cosets)
+        : cosets(thrust::raw_pointer_cast(cosets.data())) {
+    }
+
+    __device__
+    void operator()(int gen) {
+        cosets[gen] = 0;
+    }
+};
+
+// this creates rows for cosets by index of each relation table
+struct RowGen {
+    int coset;
+
+    RowGen(int coset) 
+        : coset(coset) {
+    }
+
+    __device__
+    Row operator()(int rel) {
+        return Row(rel, coset, c_rels[rel].mul * 2);
+    }
+};
+
+// determines if rows are incomplete; used to remove completed rows
+struct RowIncomplete {
+    __device__
+    bool operator()(Row r) {
+        return r.r - r.l > 1;
+    }
+};
+
+// add a row to the coset table filled with -1
+void add_row(
+        int ngens,
+        thrust::device_vector<int> &cosets) {
+    cosets.resize(cosets.size() + ngens, -1);
+};
+
+// add a new coset to the coset table, picking up where the last call left off.
+bool add_coset(
+        int ngens,
+        int *coset,
+        int *hint,
+        thrust::device_vector<int> &dcosets) {
+    int offset = *hint;
+    thrust::host_vector<int> cosets(dcosets.begin() + offset, dcosets.end());
+    *coset = dcosets.size() / ngens;
+
+    while (cosets[*hint - offset] >= 0) {
+        *hint = *hint + 1;
+        if (*hint - offset >= cosets.size()) 
+            return true;
+    }
+    int from = *hint / ngens;
+    int gen = *hint % ngens;
+    
+    add_row(ngens, dcosets);
+    
+    dcosets[*hint] = *coset;
+    dcosets[*coset * ngens + gen] = from;
+
+    return false;
+}
+
+// add a row for each relation table for some coset
+void gen_rows(
+        int coset,
+        int nrels,
+        thrust::device_vector<Row> &rows) {
+    rows.resize(rows.size() + nrels);
+
+    thrust::counting_iterator<int> counter(0);
+    thrust::transform(
+            thrust::device,
+            counter, counter + nrels,
+            rows.end() - nrels,
+            RowGen(coset));
+}
+
+// do everything. data is implicitly passed to the device via device_vector.
+thrust::device_vector<int> solve(
+        int ngens,
+        int nrels,
+        thrust::device_vector<int> subs) {
+    
+    thrust::device_vector<int> cosets;
+    thrust::device_vector<Row> rows;
+
+    // create the inital row and populate it from subs
+    add_row(ngens, cosets);
+    thrust::for_each(
+            thrust::device,
+            subs.begin(), subs.end(), 
+            CosetInitializer(cosets));
+
+    // generate initial relation table rows for coset 0
+    gen_rows(0, nrels, rows);
+
+    // these keep track of what progress has been made
+    int coset = 0;
+    int hint = 0;
+
+    // will break out later
+    while (true) {
+        // create a solver and apply it until nothing is being learned
+        Solver solve(cosets);
+        thrust::for_each(
+                thrust::device,
+                rows.begin(), rows.end(),
+                solve);
+
+        // fails if hint passes the end of the table. in that case, break.
+        bool done = add_coset(
+                ngens,
+                &coset, &hint,
+                cosets);
+        if (done) break;
+
+        // generate relation table rows for new coset
+        gen_rows(coset, nrels, rows);
+
+        // move completed rows to the end of the list and remove.
+        auto cut = thrust::partition(
+                thrust::device, 
+                rows.begin(), rows.end(), 
+                RowIncomplete());
+        rows.erase(cut, rows.end());
+    }
+
+    return cosets;
+}
+
+
+int main(int argc, const char* argv[]) {
+    Coxeter cox;
+    cox = proc_args(argc, argv);
+    std::vector<int> subs = {};
+    int nrels = cox.rels.size();
+    int ngens = cox.ngens;
+
+    cudaMemcpyToSymbol(c_ngens, &ngens, sizeof(int));
+    cudaMemcpyToSymbol(c_nrels, &nrels, sizeof(int));
+    cudaMemcpyToSymbol(c_rels, cox.rels.data(), cox.rels.size() * sizeof(Rel));
+
+    auto s = std::chrono::system_clock::now();
+    thrust::host_vector<int> cosets = solve(cox.ngens, nrels, subs);
+    auto e = std::chrono::system_clock::now();
+
+    std::chrono::duration<float> diff = e - s;
+    int order = cosets.size() / cox.ngens;
+
+    // type, arg, ngens, time, order
+    std::cout << cox.ngens << ',' << diff.count() << ',' << order << std::endl;
+
+    return 0;
+}
+
diff --git a/gpu-opt/run_gpu_opt.sh b/gpu-opt/run_gpu_opt.sh
new file mode 100644
index 0000000..9376285
--- /dev/null
+++ b/gpu-opt/run_gpu_opt.sh
@@ -0,0 +1,39 @@
+MAIN='./main'
+
+echo "type,arg,nrels,secs,cosets"
+
+TYPE=0
+for ARG in $(seq 25 25 250); do
+   $MAIN $TYPE $ARG
+done
+
+TYPE=5
+for ARG in $(seq 2 5); do
+   $MAIN $TYPE $ARG
+done
+
+TYPE=1
+$MAIN $TYPE
+
+TYPE=5
+$MAIN $TYPE 6
+
+TYPE=2
+$MAIN $TYPE
+
+TYPE=5
+$MAIN $TYPE 7
+
+TYPE=3
+$MAIN $TYPE
+
+TYPE=5
+$MAIN $TYPE 8
+$MAIN $TYPE 9
+
+TYPE=4
+$MAIN $TYPE
+
+TYPE=5
+$MAIN $TYPE 10
+
diff --git a/gpu-opt/tests.cu b/gpu-opt/tests.cu
new file mode 100644
index 0000000..cad1f96
--- /dev/null
+++ b/gpu-opt/tests.cu
@@ -0,0 +1,104 @@
+//#include <cstdio>
+//#include <cstdlib>
+//
+//#include <thrust/host_vector.h>
+//#include <thrust/device_vector.h>
+//#include <thrust/sequence.h>
+//
+//#define N 50
+//
+//__global__
+//void vector_add(float* out, float* a, float* b, int n) {
+//    for(int i = 0; i < n; i++){
+//        out[i] = a[i] + b[i];
+//    }
+//}
+//
+//int main(){
+//    thrust::host_vector<float> a(N);
+//    thrust::sequence(a.begin(), a.end());
+//
+//    thrust::host_vector<float> b(N);
+//    thrust::sequence(b.begin(), b.end());
+//    thrust::reverse(b.begin(), b.end());
+//
+//    for (int i = 0; i < N; ++i) {
+//        printf("%.1f ", a[i]);
+//    } printf("\n");
+//
+//    for (int i = 0; i < N; ++i) {
+//        printf("%.1f ", b[i]);
+//    } printf("\n");
+//
+//    thrust::device_vector<float> aD = a;
+//    thrust::device_vector<float> bD = b;
+//    thrust::device_vector<float> outD(N);
+//
+//    vector_add<<<1, 1>>>(
+//        thrust::raw_pointer_cast(&outD[0]),
+//        thrust::raw_pointer_cast(&aD[0]),
+//        thrust::raw_pointer_cast(&bD[0]),
+//        N);
+//
+//    thrust::host_vector<float> out = outD;
+//
+//    for (int i = 0; i < N; ++i) {
+//        printf("%.1f ", out[i]);
+//    } printf("\n");
+//
+//    return 0;
+//}
+
+#include <cstdio>
+#include <cstdlib>
+#include <chrono>
+
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <thrust/sequence.h>
+
+void add_proc(int *c, int *a, int *b) {
+    *c = *a + *b;
+}
+
+void test_proc(){
+    int a = 0;
+    int b = 1;
+
+    auto start = std::chrono::system_clock::now();
+
+    for (int i = 0; i < 1000000; ++i) {
+        add_proc(&a, &a, &b);
+    }
+
+    auto end = std::chrono::system_clock::now();
+
+    std::chrono::duration<float, std::micro> diff = end - start;
+
+    printf("proc: %d: 1B in %.3f micro\n", a, diff.count());
+}
+
+__global__
+void add_gpu(int *c, int *a, int *b) {
+    *c = *a + *b;
+}
+
+void test_gpu(){
+    thrust::device_vector<int> vals(2, 0);
+    vals[0] = 0;
+    vals[1] = 1;
+    printf(" gpu: %d: 1B in %.3f micro\n", vals[0], 0.0f);
+
+    int *a = thrust::raw_pointer_cast(&vals[0]);
+    int *b = thrust::raw_pointer_cast(&vals[1]);
+
+    add_gpu<<<1, 1>>>(a, a, a);
+
+    printf(" gpu: %d: 1B in %.3f micro\n", vals[0], 0.0f);
+}
+
+int main(int argc, char *argv[]) {
+    test_proc();
+    test_gpu();
+}
+
diff --git a/gpu-opt/util.h b/gpu-opt/util.h
new file mode 100644
index 0000000..c26a4a5
--- /dev/null
+++ b/gpu-opt/util.h
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <vector>
+#include <iostream>
+
+template<class T>
+std::ostream &operator<<(std::ostream &o, const thrust::host_vector<T> &vec) {
+    if (vec.size() == 0 || vec.size() > 15) 
+        return o << "host_vector{size=" << vec.size() << "}";
+
+    o << "[";
+    
+    for (int i = 0; i < vec.size() - 1; i++) o << vec[i] << ", ";
+    
+    if (vec.size() > 0) o << vec[vec.size() - 1];
+    
+    o << "]";
+
+    return o;
+}
+
+template<class T>
+std::ostream &operator<<(std::ostream &o, const thrust::device_vector<T> &vec) {
+    return o << "device_vector{size=" << vec.size() << "}";
+}
+
+struct Rel {
+    int gens[2];
+    int mul;
+};
+
+struct Coxeter {
+    int ngens;
+    std::vector<Rel> rels;
+};
+
+Coxeter make_coxeter(int ngens, const std::vector<Rel> &rels) {
+    int mults[ngens][ngens];
+
+    for (int i = 0; i < ngens; i++) {
+        for (int j = 0; j < ngens; j++) {
+            mults[i][j] = 2;
+            mults[j][i] = 2;
+        }
+    }
+
+    for (const auto &r : rels) {
+        mults[r.gens[0]][r.gens[1]] = r.mul;
+        mults[r.gens[1]][r.gens[0]] = r.mul;
+    }
+
+    std::vector<Rel> res;
+    
+    for (int i = 0; i < ngens; i++) {
+        for (int j = i + 1; j < ngens; j++) {
+            res.push_back({i, j, mults[i][j]});
+        }
+    }
+
+    return {ngens, res};
+}
diff --git a/gpu-slo/CMakeLists.txt b/gpu-slo/CMakeLists.txt
new file mode 100644
index 0000000..6ae6437
--- /dev/null
+++ b/gpu-slo/CMakeLists.txt
@@ -0,0 +1,6 @@
+cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
+project(main LANGUAGES CXX CUDA)
+
+add_executable(main main.cu)
+target_compile_features(main PUBLIC cxx_std_11)
+set_target_properties(main PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
diff --git a/gpu-slo/Makefile b/gpu-slo/Makefile
new file mode 100644
index 0000000..8d087cc
--- /dev/null
+++ b/gpu-slo/Makefile
@@ -0,0 +1,9 @@
+main : main.cu util.h
+	nvcc -o main -std=c++11 -O3 main.cu
+
+clean :
+	rm main
+
+queue : main
+	qsub -q mamba -l walltime=24:00:00 -l nodes=1:ppn=1:gpus=1 -d `pwd` run_gpu.sh
+
diff --git a/gpu-slo/groups.h b/gpu-slo/groups.h
new file mode 100644
index 0000000..8f7768a
--- /dev/null
+++ b/gpu-slo/groups.h
@@ -0,0 +1,127 @@
+#pragma once
+
+#include "util.h"
+
+#include <iostream>
+
+/*
+ * Order 4*res*res
+ */
+Coxeter torus(int res) {
+    return make_coxeter(4, {
+        {0, 1, res},
+        {2, 3, res},
+    });
+}
+
+Coxeter hypercube(int dim) {
+    std::vector<Rel> rels;
+    rels.push_back({0, 1, 4});
+    for (int d = 2; d < dim; d++) {
+        rels.push_back({d-1, d, 3});
+    }
+    return make_coxeter(dim, rels);
+}
+
+/*
+ * Order 14,400
+ */
+Coxeter H4() {
+    return make_coxeter(4, {
+        {0, 1, 5},
+        {1, 2, 3},
+        {2, 3, 3},
+    });
+}
+
+/*
+ * Order 51,840
+ */
+Coxeter E6() {
+    return make_coxeter(6, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+    });
+}
+
+/*
+ * Order 2,903,040
+ */
+Coxeter E7() {
+    return make_coxeter(7, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+    });
+}
+
+/*
+ * Order 696,729,600
+ */
+Coxeter E8() {
+    return make_coxeter(8, {
+        {0, 1, 3},
+        {1, 2, 3},
+        {2, 3, 3},
+        {2, 4, 3},
+        {4, 5, 3},
+        {5, 6, 3},
+        {6, 7, 3},
+    });
+}
+
+/*
+ * returns coxeter group based on the arguments
+ * prints out type and arguments, without an endline
+ */
+Coxeter proc_args(int argc, const char* argv[]) {
+    if (argc < 2) {
+        std::cerr << "missing type argument." << std::endl;
+        exit(EXIT_FAILURE);
+    }
+
+    int type = std::strtol(argv[1], nullptr, 10);
+    std::cout << type << ',';
+
+    int arg;
+    switch (type) {
+    case 0:
+        if (argc < 3) {
+            std::cerr << "Must provide a size for torus!" << std::endl;
+            exit(EXIT_FAILURE);
+        }
+        arg = std::strtol(argv[2], nullptr, 10);
+        std::cout << arg << ',';
+        return torus(arg);
+    case 1:
+        std::cout << -1 << ',';
+        return H4();
+    case 2:
+        std::cout << -1 << ',';
+        return E6();
+    case 3:
+        std::cout << -1 << ',';
+        return E7();
+    case 4:
+        std::cout << -1 << ',';
+        return E8();
+    case 5:
+        if (argc < 3) {
+            std::cerr << "Must provide a dimension for hypercube!" << std::endl;
+            exit(EXIT_FAILURE);
+        }
+        arg = std::strtol(argv[2], nullptr, 10);
+        std::cout << arg << ',';
+        return hypercube(arg);
+    }
+
+    std::cerr << "Not a valid type!" << std::endl;
+    exit(EXIT_FAILURE);
+}
+
diff --git a/gpu-slo/main.cu b/gpu-slo/main.cu
new file mode 100644
index 0000000..4040ad2
--- /dev/null
+++ b/gpu-slo/main.cu
@@ -0,0 +1,281 @@
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <thrust/logical.h>
+
+#include <vector>
+#include <iostream>
+#include <chrono>
+
+#include "util.h"
+#include "groups.h"
+
+struct Row {
+    int rel;
+
+    int l, r;
+
+    int from, to;
+
+    bool learning;
+
+    __host__ __device__
+    Row() : rel(0), l(0), r(0), from(0), to(0), learning(true) {}
+
+    __device__
+    Row(int rel, int cos, int size) {
+        l = 0;
+        r = size - 1;
+        
+        from = to = cos;
+        
+        this->rel = rel;
+
+        learning = true;
+    }
+};
+
+std::ostream &operator<<(std::ostream &o, const Row &r) {
+    return o << "Row[" << r.rel << "]{" << r.l << ":" << r.from << "-" << r.to << ":" << r.r << "}(" << r.learning << ")";
+}
+
+// this performs a pass on one relation table row, applying learned data to the coset table.
+struct Solver {
+    int ngens;
+    int *cosets;
+    Rel *rels;
+    
+    Solver(int ngens,
+           thrust::device_vector<int> &cosets,
+           thrust::device_vector<Rel> &rels)
+        : ngens(ngens),
+           cosets(thrust::raw_pointer_cast(cosets.data())),
+           rels(thrust::raw_pointer_cast(rels.data())) {
+        }
+    
+    __device__
+    void operator()(Row &r) {
+        if (r.r - r.l <= 0) {
+            r.learning = false;
+            return;
+        }
+        
+        while (r.r - r.l > 0) {
+            int gen = rels[r.rel].gens[r.l & 1];
+            int next = cosets[r.from * ngens + gen];
+            if (next < 0) break;
+            r.l++;
+            r.from = next;
+        }
+
+        while (r.r - r.l > 0) {
+            int gen = rels[r.rel].gens[r.r & 1];
+            int next = cosets[r.to * ngens + gen];
+            if (next < 0) break;
+            r.r--;
+            r.to = next;
+        }
+            
+        if (r.r - r.l <= 0) { 
+            int gen = rels[r.rel].gens[r.l & 1];
+            cosets[r.from * ngens + gen] = r.to;
+            cosets[r.to * ngens + gen] = r.from;
+
+            r.learning = true;
+            return;
+        }
+
+        r.learning = false;
+    }
+};
+
+// this sets the inital row in the coset table based on the subgroup generators
+struct CosetInitializer {
+    int *cosets;
+
+    CosetInitializer(thrust::device_vector<int> &cosets)
+        : cosets(thrust::raw_pointer_cast(cosets.data())) {
+    }
+
+    __device__
+    void operator()(int gen) {
+        cosets[gen] = 0;
+    }
+};
+
+// this creates rows for cosets by index of each relation table
+struct RowGen {
+    Rel *rels;
+
+    int coset;
+
+    RowGen(int coset, thrust::device_vector<Rel> &rels) 
+        : coset(coset),
+          rels(thrust::raw_pointer_cast(rels.data())) {}
+
+    __device__
+    Row operator()(int rel) {
+        return Row(rel, coset, rels[rel].mul * 2);
+    }
+};
+
+// determines if rows are incomplete; used to remove completed rows
+struct RowIncomplete {
+    __device__
+    bool operator()(Row r) {
+        return r.r - r.l > 1;
+    }
+};
+
+// re-set rows to be learning for a next pass
+struct Relearn {
+    __device__
+    void operator()(Row &r) {
+        r.learning = true;
+    }
+};
+
+// determine if rows are learning. used for exit condition
+struct Learning {
+    __device__
+    bool operator()(Row r) {
+        return r.learning;
+    }
+};
+
+// add a row to the coset table filled with -1
+void add_row(
+        int ngens,
+        thrust::device_vector<int> &cosets) {
+    cosets.resize(cosets.size() + ngens, -1);
+};
+
+// add a new coset to the coset table, picking up where the last call left off.
+// todo: this part is _real_ slow.
+bool add_coset(
+        int ngens,
+        int *coset,
+        int *hint,
+        thrust::device_vector<int> &cosets) {
+    *coset = cosets.size() / ngens;
+
+    // todo: this part especially.
+    while (cosets[*hint] >= 0) {
+        *hint = *hint + 1;
+        if (*hint >= cosets.size()) 
+            return true;
+    }
+    int from = *hint / ngens;
+    int gen = *hint % ngens;
+    
+    add_row(ngens, cosets);
+    
+    cosets[*hint] = *coset;
+    cosets[*coset * ngens + gen] = from;
+
+    return false;
+}
+
+// add a row for each relation table for some coset
+void gen_rows(
+        int coset,
+        thrust::device_vector<Rel> &rels,
+        thrust::device_vector<Row> &rows) {
+    rows.resize(rows.size() + rels.size());
+
+    thrust::counting_iterator<int> counter(0);
+    thrust::transform(
+            thrust::device,
+            counter, counter + rels.size(),
+            rows.end() - rels.size(),
+            RowGen(coset, rels));
+}
+
+// do everything. data is implicitly passed to the device via device_vector.
+thrust::device_vector<int> solve(
+        int ngens,
+        thrust::device_vector<int> subs,
+        thrust::device_vector<Rel> rels) {
+    
+    thrust::device_vector<int> cosets;
+    thrust::device_vector<Row> rows;
+
+    // create the inital row and populate it from subs
+    add_row(ngens, cosets);
+    thrust::for_each(
+            thrust::device,
+            subs.begin(), subs.end(), 
+            CosetInitializer(cosets));
+
+    // generate initial relation table rows for coset 0
+    gen_rows(0, rels, rows);
+
+    // these keep track of what progress has been made
+    int coset = 0;
+    int hint = 0;
+
+    // will break out later
+    while (true) {
+        // reset learning=true for all rows.
+        thrust::for_each(
+                thrust::device, 
+                rows.begin(), 
+                rows.end(),
+                Relearn());
+
+        // create a solver and apply it until nothing is being learned
+        Solver solve(ngens, cosets, rels);
+        while (true) {
+            thrust::for_each(
+                    thrust::device,
+                    rows.begin(), rows.end(),
+                    solve);
+
+            // if not any row is learning, then break.
+            bool r = thrust::any_of(
+                    thrust::device,
+                    rows.begin(), rows.end(),
+                    Learning());
+            if (!r) break;
+        }
+
+
+        // fails if hint passes the end of the table. in that case, break.
+        bool done = add_coset(
+                ngens,
+                &coset, &hint,
+                cosets);
+        if (done) break;
+
+        // generate relation table rows for new coset
+        gen_rows(coset, rels, rows);
+
+        // move completed rows to the end of the list and remove.
+        auto cut = thrust::partition(
+                thrust::device, 
+                rows.begin(), rows.end(), 
+                RowIncomplete());
+        rows.erase(cut, rows.end());
+    }
+
+    return cosets;
+}
+
+
+int main(int argc, const char* argv[]) {
+    Coxeter cox;
+    cox = proc_args(argc, argv);
+    std::vector<int> subs = {};
+
+    auto s = std::chrono::system_clock::now();
+    thrust::host_vector<int> cosets = solve(cox.ngens, subs, cox.rels);
+    auto e = std::chrono::system_clock::now();
+
+    std::chrono::duration<float> diff = e - s;
+    int order = cosets.size() / cox.ngens;
+
+    // type, arg, ngens, time, order
+    std::cout << cox.ngens << ',' << diff.count() << ',' << order << std::endl;
+
+    return 0;
+}
+
diff --git a/gpu-slo/run_gpu.sh b/gpu-slo/run_gpu.sh
new file mode 100644
index 0000000..9376285
--- /dev/null
+++ b/gpu-slo/run_gpu.sh
@@ -0,0 +1,39 @@
+MAIN='./main'
+
+echo "type,arg,nrels,secs,cosets"
+
+TYPE=0
+for ARG in $(seq 25 25 250); do
+   $MAIN $TYPE $ARG
+done
+
+TYPE=5
+for ARG in $(seq 2 5); do
+   $MAIN $TYPE $ARG
+done
+
+TYPE=1
+$MAIN $TYPE
+
+TYPE=5
+$MAIN $TYPE 6
+
+TYPE=2
+$MAIN $TYPE
+
+TYPE=5
+$MAIN $TYPE 7
+
+TYPE=3
+$MAIN $TYPE
+
+TYPE=5
+$MAIN $TYPE 8
+$MAIN $TYPE 9
+
+TYPE=4
+$MAIN $TYPE
+
+TYPE=5
+$MAIN $TYPE 10
+
diff --git a/gpu-slo/tests.cu b/gpu-slo/tests.cu
new file mode 100644
index 0000000..cad1f96
--- /dev/null
+++ b/gpu-slo/tests.cu
@@ -0,0 +1,104 @@
+//#include <cstdio>
+//#include <cstdlib>
+//
+//#include <thrust/host_vector.h>
+//#include <thrust/device_vector.h>
+//#include <thrust/sequence.h>
+//
+//#define N 50
+//
+//__global__
+//void vector_add(float* out, float* a, float* b, int n) {
+//    for(int i = 0; i < n; i++){
+//        out[i] = a[i] + b[i];
+//    }
+//}
+//
+//int main(){
+//    thrust::host_vector<float> a(N);
+//    thrust::sequence(a.begin(), a.end());
+//
+//    thrust::host_vector<float> b(N);
+//    thrust::sequence(b.begin(), b.end());
+//    thrust::reverse(b.begin(), b.end());
+//
+//    for (int i = 0; i < N; ++i) {
+//        printf("%.1f ", a[i]);
+//    } printf("\n");
+//
+//    for (int i = 0; i < N; ++i) {
+//        printf("%.1f ", b[i]);
+//    } printf("\n");
+//
+//    thrust::device_vector<float> aD = a;
+//    thrust::device_vector<float> bD = b;
+//    thrust::device_vector<float> outD(N);
+//
+//    vector_add<<<1, 1>>>(
+//        thrust::raw_pointer_cast(&outD[0]),
+//        thrust::raw_pointer_cast(&aD[0]),
+//        thrust::raw_pointer_cast(&bD[0]),
+//        N);
+//
+//    thrust::host_vector<float> out = outD;
+//
+//    for (int i = 0; i < N; ++i) {
+//        printf("%.1f ", out[i]);
+//    } printf("\n");
+//
+//    return 0;
+//}
+
+#include <cstdio>
+#include <cstdlib>
+#include <chrono>
+
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <thrust/sequence.h>
+
+void add_proc(int *c, int *a, int *b) {
+    *c = *a + *b;
+}
+
+void test_proc(){
+    int a = 0;
+    int b = 1;
+
+    auto start = std::chrono::system_clock::now();
+
+    for (int i = 0; i < 1000000; ++i) {
+        add_proc(&a, &a, &b);
+    }
+
+    auto end = std::chrono::system_clock::now();
+
+    std::chrono::duration<float, std::micro> diff = end - start;
+
+    printf("proc: %d: 1B in %.3f micro\n", a, diff.count());
+}
+
+__global__
+void add_gpu(int *c, int *a, int *b) {
+    *c = *a + *b;
+}
+
+void test_gpu(){
+    thrust::device_vector<int> vals(2, 0);
+    vals[0] = 0;
+    vals[1] = 1;
+    printf(" gpu: %d: 1B in %.3f micro\n", vals[0], 0.0f);
+
+    int *a = thrust::raw_pointer_cast(&vals[0]);
+    int *b = thrust::raw_pointer_cast(&vals[1]);
+
+    add_gpu<<<1, 1>>>(a, a, a);
+
+    printf(" gpu: %d: 1B in %.3f micro\n", vals[0], 0.0f);
+}
+
+int main(int argc, char *argv[]) {
+    test_proc();
+    test_gpu();
+}
+
diff --git a/gpu-slo/util.h b/gpu-slo/util.h
new file mode 100644
index 0000000..c26a4a5
--- /dev/null
+++ b/gpu-slo/util.h
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <thrust/device_vector.h>
+#include <thrust/host_vector.h>
+#include <vector>
+#include <iostream>
+
+template<class T>
+std::ostream &operator<<(std::ostream &o, const thrust::host_vector<T> &vec) {
+    if (vec.size() == 0 || vec.size() > 15) 
+        return o << "host_vector{size=" << vec.size() << "}";
+
+    o << "[";
+    
+    for (int i = 0; i < vec.size() - 1; i++) o << vec[i] << ", ";
+    
+    if (vec.size() > 0) o << vec[vec.size() - 1];
+    
+    o << "]";
+
+    return o;
+}
+
+template<class T>
+std::ostream &operator<<(std::ostream &o, const thrust::device_vector<T> &vec) {
+    return o << "device_vector{size=" << vec.size() << "}";
+}
+
+struct Rel {
+    int gens[2];
+    int mul;
+};
+
+struct Coxeter {
+    int ngens;
+    std::vector<Rel> rels;
+};
+
+Coxeter make_coxeter(int ngens, const std::vector<Rel> &rels) {
+    int mults[ngens][ngens];
+
+    for (int i = 0; i < ngens; i++) {
+        for (int j = 0; j < ngens; j++) {
+            mults[i][j] = 2;
+            mults[j][i] = 2;
+        }
+    }
+
+    for (const auto &r : rels) {
+        mults[r.gens[0]][r.gens[1]] = r.mul;
+        mults[r.gens[1]][r.gens[0]] = r.mul;
+    }
+
+    std::vector<Rel> res;
+    
+    for (int i = 0; i < ngens; i++) {
+        for (int j = i + 1; j < ngens; j++) {
+            res.push_back({i, j, mults[i][j]});
+        }
+    }
+
+    return {ngens, res};
+}
diff --git a/hpclog/log-gpu-opt-o3 b/hpclog/log-gpu-opt-o3
new file mode 100644
index 0000000..5eb86e4
--- /dev/null
+++ b/hpclog/log-gpu-opt-o3
@@ -0,0 +1,20 @@
+type,arg,nrels,secs,cosets
+0,25,4,0.38605,2500
+0,50,4,1.6768,10000
+0,75,4,7.64476,22500
+0,100,4,19.5545,40000
+0,125,4,44.443,62500
+0,150,4,82.011,90000
+0,175,4,124.86,122500
+0,200,4,174.554,160000
+0,225,4,230.284,202500
+0,250,4,288.602,250000
+5,2,2,0.00299317,8
+5,3,3,0.00585187,48
+5,4,4,0.0421133,384
+5,5,5,0.439017,3840
+1,-1,4,1.63292,14400
+5,6,6,14.8635,46080
+2,-1,6,17.6105,51840
+5,7,7,734.339,645120
+3,-1,7,6729.78,2903040
diff --git a/hpclog/log-gpu-slo b/hpclog/log-gpu-slo
new file mode 100644
index 0000000..d15e575
--- /dev/null
+++ b/hpclog/log-gpu-slo
@@ -0,0 +1,20 @@
+type,arg,nrels,secs,cosets
+0,25,4,1.11799,2500
+0,50,4,4.33001,10000
+0,75,4,10.4824,22500
+0,100,4,23.661,40000
+0,125,4,56.4877,62500
+0,150,4,103.069,90000
+0,175,4,163.208,122500
+0,200,4,246.577,160000
+0,225,4,354.257,202500
+0,250,4,494.363,250000
+5,2,2,0.0894585,8
+5,3,3,0.0971892,48
+5,4,4,0.189169,384
+5,5,5,1.35973,3840
+1,-1,4,4.76805,14400
+5,6,6,22.5943,46080
+2,-1,6,25.1293,51840
+5,7,7,1350.11,645120
+3,-1,7,14657.7,2903040
diff --git a/hpclog/log-gpu-slo-o3 b/hpclog/log-gpu-slo-o3
new file mode 100644
index 0000000..394b80d
--- /dev/null
+++ b/hpclog/log-gpu-slo-o3
@@ -0,0 +1,20 @@
+type,arg,nrels,secs,cosets
+0,25,4,1.03598,2500
+0,50,4,3.89485,10000
+0,75,4,9.97716,22500
+0,100,4,22.3504,40000
+0,125,4,60.0055,62500
+0,150,4,107.5,90000
+0,175,4,169.508,122500
+0,200,4,256.386,160000
+0,225,4,367.362,202500
+0,250,4,508.957,250000
+5,2,2,0.0747127,8
+5,3,3,0.0833176,48
+5,4,4,0.166347,384
+5,5,5,1.24269,3840
+1,-1,4,4.30961,14400
+5,6,6,20.8581,46080
+2,-1,6,23.1169,51840
+5,7,7,1382.62,645120
+3,-1,7,14806.1,2903040