split for gpu optimization

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)

gpu-opt/Makefile

@@ -0,0 +1,6 @@
+main : main.cu util.h
+	nvcc -o main -std=c++11 main.cu
+
+clean :
+	rm main
+

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);
+}
+

gpu-opt/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;
+}
+

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();
+}
+

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};
+}