// proto-A-gpu-query.cu
// ---------------------------------------------------------------------------
// PROTOTYPE A: GPU AtomSpace "Query" as the set-intersection primitive.
//
// Grounding (HARDWARE-KB.md §G): across every GPU subgraph matcher the hot op
// (~95% of runtime) is parallel SET INTERSECTION of adjacency / incoming-set
// lists. AtomSpace clause-grounding with a shared variable IS that op.
//
// Minimal faithful instance: "common neighbours of a and b" =
// find all $X such that Edge(a,$X) AND Edge(b,$X)
// = intersect( adj[a], adj[b] )
// That is a 2-clause conjunctive query sharing one variable $X. This is the
// seed of worst-case-optimal multiway join (cuMatch); >2 clauses = intersect
// more lists.
//
// Discipline: standalone (no atomspace deps), CPU-verified, runs on RTX 2070.
// Build: nvcc -O2 -arch=sm_75 -o proto-A proto-A-gpu-query.cu
// Run: ./proto-A
// ---------------------------------------------------------------------------
#include <cstdio>
#include <cstdint>
#include <vector>
#include <algorithm>
#include <random>
#include <cuda_runtime.h>
#define CUDA_OK(call) do { cudaError_t e=(call); if(e!=cudaSuccess){ \
printf("CUDA ERROR %s:%d: %s\n",__FILE__,__LINE__,cudaGetErrorString(e)); \
return -1; } } while(0)
// CSR graph: row_off[N+1], col_idx[] with each row's neighbours sorted ascending.
__device__ __forceinline__ int bsearch_dev(const uint32_t* arr, int len, uint32_t key){
int lo=0, hi=len-1;
while(lo<=hi){ int mid=(lo+hi)>>1; uint32_t v=arr[mid];
if(v==key) return mid; if(v<key) lo=mid+1; else hi=mid-1; }
return -1;
}
// One BLOCK per query. Threads stride the SMALLER list, binary-search the larger.
// (Iterate-smaller / search-larger is the standard sorted-list intersection.)
__global__ void intersect_queries(
const uint32_t* __restrict__ row_off,
const uint32_t* __restrict__ col_idx,
const uint32_t* __restrict__ qa,
const uint32_t* __restrict__ qb,
int num_queries,
uint32_t* __restrict__ out_buf, int out_stride,
int* __restrict__ out_count)
{
int q = blockIdx.x;
if(q >= num_queries) return;
uint32_t a = qa[q], b = qb[q];
const uint32_t* A = col_idx + row_off[a]; int lenA = row_off[a+1]-row_off[a];
const uint32_t* B = col_idx + row_off[b]; int lenB = row_off[b+1]-row_off[b];
const uint32_t* S; int lenS; const uint32_t* L; int lenL;
if(lenA <= lenB){ S=A; lenS=lenA; L=B; lenL=lenB; }
else { S=B; lenS=lenB; L=A; lenL=lenA; }
for(int i=threadIdx.x; i<lenS; i+=blockDim.x){
uint32_t key = S[i];
if(bsearch_dev(L, lenL, key) >= 0){
int idx = atomicAdd(&out_count[q], 1); // claim an output slot
if(idx < out_stride) out_buf[q*out_stride + idx] = key;
}
}
}
// ---- CPU reference -------------------------------------------------------
static std::vector<uint32_t> cpu_intersect(
const std::vector<uint32_t>& row_off,
const std::vector<uint32_t>& col_idx, uint32_t a, uint32_t b){
std::vector<uint32_t> r;
std::set_intersection(
col_idx.begin()+row_off[a], col_idx.begin()+row_off[a+1],
col_idx.begin()+row_off[b], col_idx.begin()+row_off[b+1],
std::back_inserter(r));
return r;
}
static int fails = 0;
static void check(const char* name, std::vector<uint32_t> got, std::vector<uint32_t> exp){
std::sort(got.begin(), got.end());
std::sort(exp.begin(), exp.end());
if(got==exp){ printf("PASS %s (|result|=%zu)\n", name, exp.size()); }
else { fails++; printf("FAIL %s\n expected:", name);
for(auto x:exp) printf(" %u",x); printf("\n got: ");
for(auto x:got) printf(" %u",x); printf("\n"); }
}
// Run all queries on GPU, return per-query result vectors.
static int run_gpu(const std::vector<uint32_t>& row_off,
const std::vector<uint32_t>& col_idx,
const std::vector<uint32_t>& qa,
const std::vector<uint32_t>& qb,
int out_stride,
std::vector<std::vector<uint32_t>>& results,
float* ms_out)
{
int nq = qa.size();
uint32_t *d_row,*d_col,*d_qa,*d_qb,*d_out; int* d_cnt;
CUDA_OK(cudaMalloc(&d_row, row_off.size()*4));
CUDA_OK(cudaMalloc(&d_col, col_idx.size()*4));
CUDA_OK(cudaMalloc(&d_qa, nq*4));
CUDA_OK(cudaMalloc(&d_qb, nq*4));
CUDA_OK(cudaMalloc(&d_out, (size_t)nq*out_stride*4));
CUDA_OK(cudaMalloc(&d_cnt, nq*4));
CUDA_OK(cudaMemcpy(d_row,row_off.data(),row_off.size()*4,cudaMemcpyHostToDevice));
CUDA_OK(cudaMemcpy(d_col,col_idx.data(),col_idx.size()*4,cudaMemcpyHostToDevice));
CUDA_OK(cudaMemcpy(d_qa,qa.data(),nq*4,cudaMemcpyHostToDevice));
CUDA_OK(cudaMemcpy(d_qb,qb.data(),nq*4,cudaMemcpyHostToDevice));
CUDA_OK(cudaMemset(d_cnt,0,nq*4));
cudaEvent_t t0,t1; cudaEventCreate(&t0); cudaEventCreate(&t1);
cudaEventRecord(t0);
intersect_queries<<<nq,128>>>(d_row,d_col,d_qa,d_qb,nq,d_out,out_stride,d_cnt);
cudaEventRecord(t1);
CUDA_OK(cudaDeviceSynchronize());
cudaEventElapsedTime(ms_out,t0,t1);
std::vector<int> cnt(nq);
std::vector<uint32_t> out((size_t)nq*out_stride);
CUDA_OK(cudaMemcpy(cnt.data(),d_cnt,nq*4,cudaMemcpyDeviceToHost));
CUDA_OK(cudaMemcpy(out.data(),d_out,(size_t)nq*out_stride*4,cudaMemcpyDeviceToHost));
results.assign(nq,{});
for(int q=0;q<nq;q++){ int c=std::min(cnt[q],out_stride);
for(int i=0;i<c;i++) results[q].push_back(out[(size_t)q*out_stride+i]); }
cudaFree(d_row);cudaFree(d_col);cudaFree(d_qa);cudaFree(d_qb);
cudaFree(d_out);cudaFree(d_cnt);
return 0;
}
int main(){
printf("=== Prototype A: GPU AtomSpace Query (set-intersection primitive) ===\n\n");
// ---- Test 1: tiny hand-built graph with KNOWN answers ----------------
// adj[0]={2,3,4,5} adj[1]={3,5,6,7} -> common(0,1)={3,5}
// adj[2]={3,4,8} adj[3]={4,5,8,9} -> common(2,3)={4,8}
// common(0,2)={3,4} common(0,3)={4,5}
std::vector<std::vector<uint32_t>> adj = {
{2,3,4,5}, {3,5,6,7}, {3,4,8}, {4,5,8,9},
{}, {}, {}, {}, {}, {} };
int N = adj.size();
std::vector<uint32_t> row_off(N+1,0), col_idx;
for(int i=0;i<N;i++){ std::sort(adj[i].begin(),adj[i].end());
row_off[i+1]=row_off[i]+adj[i].size();
for(auto x:adj[i]) col_idx.push_back(x); }
std::vector<uint32_t> qa={0,2,0,0}, qb={1,3,2,3};
std::vector<std::vector<uint32_t>> gpu; float ms;
if(run_gpu(row_off,col_idx,qa,qb,64,gpu,&ms)) return 1;
check("common(0,1)", gpu[0], {3,5});
check("common(2,3)", gpu[1], {4,8});
check("common(0,2)", gpu[2], {3,4});
check("common(0,3)", gpu[3], {4,5});
// ---- Test 2: larger random graph, GPU vs CPU reference ---------------
const int RN = 20000, AVG_DEG = 32, NQ = 4096;
std::mt19937 rng(12345); // fixed seed (no Math.random ban here, but deterministic anyway)
std::uniform_int_distribution<uint32_t> node(0,RN-1);
std::vector<std::vector<uint32_t>> radj(RN);
for(int i=0;i<RN;i++){
int d = AVG_DEG/2 + (rng()% AVG_DEG);
for(int k=0;k<d;k++) radj[i].push_back(node(rng));
std::sort(radj[i].begin(),radj[i].end());
radj[i].erase(std::unique(radj[i].begin(),radj[i].end()),radj[i].end());
}
std::vector<uint32_t> rrow(RN+1,0), rcol;
for(int i=0;i<RN;i++){ rrow[i+1]=rrow[i]+radj[i].size();
for(auto x:radj[i]) rcol.push_back(x); }
std::vector<uint32_t> rqa(NQ), rqb(NQ);
for(int i=0;i<NQ;i++){ rqa[i]=node(rng); rqb[i]=node(rng); }
std::vector<std::vector<uint32_t>> rgpu; float rms;
if(run_gpu(rrow,rcol,rqa,rqb,256,rgpu,&rms)) return 1;
int mism=0; size_t total=0;
for(int q=0;q<NQ;q++){
auto exp = cpu_intersect(rrow,rcol,rqa[q],rqb[q]);
auto got = rgpu[q]; std::sort(got.begin(),got.end());
total += exp.size();
if(got!=exp) mism++;
}
printf("\nTest 2: %d queries on %d-node graph (~%d avg deg)\n", NQ, RN, AVG_DEG);
printf(" GPU kernel time: %.3f ms (%.1f Mqueries/s)\n", rms, NQ/rms/1000.0);
printf(" total matches: %zu mismatches vs CPU: %d\n", total, mism);
if(mism){ fails++; printf("FAIL GPU disagrees with CPU on %d queries\n",mism); }
else printf("PASS GPU == CPU on all %d queries\n", NQ);
printf("\n%s\n", fails? "SOME CHECKS FAILED" : "ALL CHECKS PASSED — Query primitive works on GPU.");
return fails?1:0;
}