diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c8cfa4f0..b19a26caf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -360,6 +360,7 @@ add_library(spla SHARED src/binding/c_scalar.cpp src/binding/c_type.cpp src/binding/c_vector.cpp + src/binding/c_descriptor.cpp # C++ optional part ${SRC_OPENCL}) diff --git a/examples/sssp.cpp b/examples/sssp.cpp index d5ddf04fb..7a5b40289 100644 --- a/examples/sssp.cpp +++ b/examples/sssp.cpp @@ -27,10 +27,8 @@ #include "common.hpp" #include "options.hpp" - #include - int main(int argc, const char* const* argv) { std::shared_ptr options = make_options("sssp", "sssp (single-source shortest path) algorithm with spla library"); cxxopts::ParseResult args; @@ -73,12 +71,12 @@ int main(int argc, const char* const* argv) { desc->set_traversal_mode(static_cast(args[OPT_PUSH_PULL].as() - 1)); desc->set_front_factor(args[OPT_FRONT_FACTOR].as()); - const spla::uint w = 1.0f; - const auto& Ai = loader.get_Ai(); - const auto& Aj = loader.get_Aj(); + const auto& Ai = loader.get_Ai(); + const auto& Aj = loader.get_Aj(); + const auto& Aw = loader.get_Aw(); for (std::size_t k = 0; k < loader.get_n_values(); ++k) { - A->set_float(Ai[k], Aj[k], w); + A->set_float(Ai[k], Aj[k], Aw[k]); } const int n_iters = args[OPT_NITERS].as(); @@ -114,7 +112,7 @@ int main(int argc, const char* const* argv) { for (std::size_t k = 0; k < loader.get_n_values(); ++k) { ref_Ai[Ai[k]].push_back(Aj[k]); - ref_Ax[Ai[k]].push_back(w); + ref_Ax[Ai[k]].push_back(Aw[k]); } timer_ref.lap_begin(); @@ -141,4 +139,4 @@ int main(int argc, const char* const* argv) { std::cout << std::endl; return 0; -} \ No newline at end of file +} diff --git a/examples/tc.cpp b/examples/tc.cpp index 90b07e3cd..9ab8ce47a 100644 --- a/examples/tc.cpp +++ b/examples/tc.cpp @@ -80,7 +80,7 @@ int main(int argc, const char* const* argv) { const auto& Aj = loader.get_Aj(); for (std::size_t k = 0; k < loader.get_n_values(); ++k) { - if (Ai[k] > Aj[k]) { + if (Ai[k] < Aj[k]) { A->set_int(Ai[k], Aj[k], 1); } } @@ -118,7 +118,7 @@ int main(int argc, const char* const* argv) { // for tc it is important to have no duplicates for (std::size_t k = 0; k < loader.get_n_values(); ++k) { - if (Ai[k] > Aj[k]) { + if (Ai[k] < Aj[k]) { ref_A_set[Ai[k]].insert(Aj[k]); } } diff --git a/include/spla.h b/include/spla.h index 907865fab..336ba97f8 100644 --- a/include/spla.h +++ b/include/spla.h @@ -334,6 +334,7 @@ SPLA_API spla_Status spla_Vector_get_float(spla_Vector v, spla_uint row_id, floa SPLA_API spla_Status spla_Vector_build(spla_Vector v, spla_MemView keys, spla_MemView values); SPLA_API spla_Status spla_Vector_read(spla_Vector v, spla_MemView* keys, spla_MemView* values); SPLA_API spla_Status spla_Vector_clear(spla_Vector v); +SPLA_API spla_Status spla_Vector_fill_with(spla_Vector v, spla_Scalar value); ////////////////////////////////////////////////////////////////////////////////////// @@ -389,6 +390,13 @@ SPLA_API spla_Status spla_Exec_v_count_mf(spla_Scalar r, spla_Vector v, spla_Des ////////////////////////////////////////////////////////////////////////////////////// +/* Descriptor object methods */ + +SPLA_API spla_Status spla_Descriptor_make(spla_Descriptor* descriptor); +SPLA_API spla_Status spla_Descriptor_set_struct_only(spla_Descriptor descriptor, spla_bool value); +SPLA_API spla_Status spla_Descriptor_set_early_exit(spla_Descriptor descriptor, spla_bool value); + +////////////////////////////////////////////////////////////////////////////////////// #if defined(__cplusplus) } #endif diff --git a/include/spla/io.hpp b/include/spla/io.hpp index 90b0d5c54..c8d5049ab 100644 --- a/include/spla/io.hpp +++ b/include/spla/io.hpp @@ -61,7 +61,7 @@ namespace spla { */ SPLA_API bool load(std::filesystem::path file_path, bool offset_indices = true, - bool make_undirected = true, + bool make_undirected = false, bool remove_loops = true); /** @@ -80,8 +80,9 @@ namespace spla { [[nodiscard]] SPLA_API const std::vector& get_Ai() const; [[nodiscard]] SPLA_API const std::vector& get_Aj() const; - [[nodiscard]] SPLA_API uint get_n_rows() const; - [[nodiscard]] SPLA_API uint get_n_cols() const; + [[nodiscard]] SPLA_API const std::vector& get_Aw() const; + [[nodiscard]] SPLA_API uint get_n_rows() const; + [[nodiscard]] SPLA_API uint get_n_cols() const; [[nodiscard]] SPLA_API std::size_t get_n_values() const; private: @@ -89,6 +90,7 @@ namespace spla { std::filesystem::path m_file_path; std::vector m_Ai; std::vector m_Aj; + std::vector m_Aw; bool m_base_is_zero = false; uint m_n_rows = 0; uint m_n_cols = 0; diff --git a/python/algorithms/README.md b/python/algorithms/README.md new file mode 100644 index 000000000..3c74a9863 --- /dev/null +++ b/python/algorithms/README.md @@ -0,0 +1,111 @@ +# Graph Algorithms with SPLA + +This project provides a Python implementation of graph algorithms using sparse +linear algebra primitives from SPLA on GPU + +Tests (`test_compare.py`) checks that both implementations produce identical results. + +Run tests: +```bash +python -m unittest test_compare.py -v +``` +--- + +## Supported Algorithms + +| Algorithm | Flag | +|-----------|------| +| BFS (Breadth-First Search) | `--algo bfs` | +| SSSP (Single-Source Shortest Paths) | `--algo sssp` | +| PageRank | `--algo pr` | +| Triangle Counting | `--algo tc` | + +--- + +SSSP should be used with undirected graphs only + +## Installation + +```bash +cd python +python -m venv venv +source venv/bin/activate +pip install -e . +cd algorithms +``` + +--- + +## Usage +```text +usage: main_spla.py [-h] --algo {bfs,sssp,pr,tc} [-m MATRIX] [-v VECTORS] [-o OUTPUT] [-s START] [-a ALPHA] [-e EPS] + +options: + -h, --help show this help message and exit + --algo {bfs,sssp,pr,tc} + Algorithm to run: + bfs - Breadth-First Search + sssp - Single-Source Shortest Paths + pr - PageRank + tc - Triangle Counting + -m MATRIX, --matrix MATRIX + Path to graph in Matrix Market format (.mtx) + -v VECTORS, --vectors VECTORS + Path to graph in vectors format (.txt) + -o OUTPUT, --output OUTPUT + Output file path (default: result.txt) + -s START, --start START + Start vertex (used in bfs, sssp; default: 0) + -a ALPHA, --alpha ALPHA + Damping factor (used in pr; default: 0.85) + -e EPS, --eps EPS + Convergence tolerance (used in pr; default: 1e-4) +``` +--- + +## Examples + +```bash +# BFS +python main_spla.py --algo bfs -m graph.mtx -s 0 -o bfs_result.txt + +# SSSP +python main_spla.py --algo sssp -v graph.txt -s 0 -o sssp_result.txt + +# PageRank +python main_spla.py --algo pr -v graph.txt -a 0.85 -e 1e-6 -o pr_result.txt + +# Triangle Counting +python main_spla.py --algo tc -m graph.mtx -o tc_result.txt +``` + +--- + +## Input format + +Two options: + +- `-v graph.txt` — Vectors format +- `-m graph.mtx` — Matrix Market format (from https://sparse.tamu.edu/) + + +### graph.txt + +Example +```text +0 1 2 # I: source vertex indices +1 2 0 # J: target vertex indices +5 3 2 # V: edge weights +``` + +### graph.mtx (Matrix Market) +Example +```text +%%MatrixMarket matrix coordinate real general +% rows cols nnz +4 4 4 #rows, cols, non-zero values +1 2 1.0 #source vertex, target vertex, weight +2 3 2.0 +3 4 3.0 +4 1 4.0 +``` \ No newline at end of file diff --git a/python/algorithms/bfs_classic.py b/python/algorithms/bfs_classic.py new file mode 100644 index 000000000..62967c99f --- /dev/null +++ b/python/algorithms/bfs_classic.py @@ -0,0 +1,16 @@ +from collections import deque + +def bfs_naive(start, graph, n): + depth = [None] * n + q = deque() + depth[start] = 0 + q.append(start) + + while q: + u = q.popleft() + for v in graph[u]: + if depth[v] is None: + depth[v] = depth[u] + 1 + q.append(v) + + return depth \ No newline at end of file diff --git a/python/algorithms/bfs_spla.py b/python/algorithms/bfs_spla.py new file mode 100644 index 000000000..a900104f7 --- /dev/null +++ b/python/algorithms/bfs_spla.py @@ -0,0 +1,41 @@ +from pyspla import INT, Matrix, Scalar, Vector, Descriptor + +def bfs_spla(s: int, A: Matrix, + push_only: bool = False, + pull_only: bool = False, + push_pull: bool = True, + front_factor: float = 0.05): + n = A.n_rows + if not (push_only or pull_only or push_pull): + push_only = True + + v = Vector(n, INT) + front_prev = Vector(n, INT) + front_new = Vector(n, INT) + front_prev.set(s, 1) + + depth = Scalar(INT, 0) + count = 0 + front_size = 1 + + desc = Descriptor() + desc.set_early_exit(True) + desc.set_struct_only(True) + + while front_size > 0: + depth += 1 + count += front_size + v.assign(front_prev, depth, op_assign=INT.SECOND, op_select=INT.NQZERO) + front_density = front_size / n + is_push_better = front_density <= front_factor + + if push_only or (push_pull and is_push_better): + front_prev.vxm(v, A,op_mult=INT.LAND, op_add=INT.LOR, op_select=INT.EQZERO, out=front_new, desc=desc) + else: + A.mxv(v, front_prev, op_mult=INT.LAND, op_add=INT.LOR, op_select=INT.EQZERO, out=front_new, desc=desc) + + front_size = front_new.count_mf().get() + + front_prev, front_new = front_new, front_prev + + return v, count, depth.get() \ No newline at end of file diff --git a/python/algorithms/graph.txt b/python/algorithms/graph.txt new file mode 100644 index 000000000..166ad99ca --- /dev/null +++ b/python/algorithms/graph.txt @@ -0,0 +1,3 @@ +0 1 2 +1 2 0 +5 3 2 diff --git a/python/algorithms/graph_spla.py b/python/algorithms/graph_spla.py new file mode 100644 index 000000000..6c1c8f3f3 --- /dev/null +++ b/python/algorithms/graph_spla.py @@ -0,0 +1,56 @@ +from scipy.io import mmread +from pyspla import Matrix, INT, FLOAT + + +def sort_by_ij(I, J, V): + return map(list, zip(*sorted(zip(I, J, V)))) + + +def read_spla(filepath, dtype=FLOAT): + coo = mmread(filepath).tocoo() + n = int(max(coo.shape)) + rows = coo.row.astype(int).tolist() + cols = coo.col.astype(int).tolist() + vals = coo.data.astype(float).tolist() + if dtype == INT: + vals = [int(v) for v in vals] + rows, cols, vals = sort_by_ij(rows, cols, vals) + return Matrix.from_lists(rows, cols, vals, shape=(n, n), dtype=dtype) + + +def read_vectors_spla(filepath, dtype=FLOAT): + with open(filepath) as f: + I = list(map(int, f.readline().split())) + J = list(map(int, f.readline().split())) + V = list(map(float, f.readline().split())) + if dtype == INT: + V = [int(v) for v in V] + n = max(max(I), max(J)) + 1 + I, J, V = sort_by_ij(I, J, V) + return Matrix.from_lists(I, J, V, shape=(n, n), dtype=dtype) + + +def read_mtx_pr_spla(file_path, alpha=0.85): + coo = mmread(file_path).tocoo() + n = int(max(coo.shape)) + rows = coo.row.astype(int).tolist() + cols = coo.col.astype(int).tolist() + out_deg = [0] * n + for i in rows: + out_deg[i] += 1 + weights = [alpha / out_deg[i] if out_deg[i] > 0 else 0.0 for i in rows] + rows, cols, weights = sort_by_ij(rows, cols, weights) + return Matrix.from_lists(rows, cols, weights, shape=(n, n), dtype=FLOAT) + + +def read_vectors_pr_spla(file_path, alpha=0.85): + with open(file_path) as f: + I = list(map(int, f.readline().split())) + J = list(map(int, f.readline().split())) + n = max(max(I), max(J)) + 1 + out_deg = [0] * n + for i in I: + out_deg[i] += 1 + V = [alpha / out_deg[i] if out_deg[i] > 0 else 0.0 for i in I] + I, J, V = sort_by_ij(I, J, V) + return Matrix.from_lists(I, J, V, shape=(n, n), dtype=FLOAT) \ No newline at end of file diff --git a/python/algorithms/main_spla.py b/python/algorithms/main_spla.py new file mode 100644 index 000000000..06834317d --- /dev/null +++ b/python/algorithms/main_spla.py @@ -0,0 +1,122 @@ +import argparse +import sys +import math +from pathlib import Path +from bfs_spla import bfs_spla +from sssp_spla import sssp_spla +from pr_spla import pr_spla +from tc_spla import tc_spla +from graph_spla import ( + read_spla, + read_vectors_spla, + read_mtx_pr_spla, + read_vectors_pr_spla, +) +from pyspla import INT, FLOAT + +INF = math.inf + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument( + "--algo", + choices=["bfs", "sssp", "pr", "tc"], + required=True, + help="Algorithm to run: " + "bfs - Breadth-First Search, " + "sssp - Single-Source Shortest Paths, " + "pr - PageRank, " + "tc - Triangle Counting", + ) + parser.add_argument( + "-m", "--matrix", + type=Path, + help="Path to graph in Matrix Market format (.mtx)", + ) + parser.add_argument( + "-v", "--vectors", + type=Path, + help="Path to graph in vectors format (.txt)", + ) + parser.add_argument( + "-o", "--output", + type=Path, + default="result.txt", + help="Output file path", + ) + parser.add_argument( + "-s", "--start", + type=int, + default=0, + help="Start vertex(default: 0)", + ) + parser.add_argument( + "-a", "--alpha", + type=float, + default=0.85, + help="Damping factor(default: 0.85)", + ) + parser.add_argument( + "-e", "--eps", + type=float, + default=1e-4, + help="Convergence tolerance(default: 1e-4)", + ) + args = parser.parse_args() + + if not args.matrix and not args.vectors: + sys.exit("Error: no graph file provided") + + if args.algo == "bfs": + if args.matrix: + A = read_spla(str(args.matrix), dtype=INT) + v, count, depth = bfs_spla(args.start, A) + else: + A = read_vectors_spla(str(args.vectors), dtype=INT) + v, count, depth = bfs_spla(args.start, A, push_only=True, pull_only=False, push_pull=False) + with open(args.output, "w") as f_out: + f_out.write(f"Reached vertices: {count}\n") + f_out.write(f"Max depth: {depth - 1}\n") + idx, vals = v.to_lists() + for k in range(len(idx)): + f_out.write(f"{idx[k]} {vals[k] - 1}\n") + print(f"Result saved to {args.output}") + + elif args.algo == "sssp": + if args.matrix: + A = read_spla(str(args.matrix), dtype=FLOAT) + v = sssp_spla(args.start, A) + else: + A = read_vectors_spla(str(args.vectors), dtype=FLOAT) + v = sssp_spla(args.start, A, push_only=True, pull_only=False, push_pull=False) + with open(args.output, "w") as f_out: + idx, vals = v.to_lists() + for k in range(len(idx)): + if vals[k] < INF: + f_out.write(f"{idx[k]} {vals[k]}\n") + print(f"Result saved to {args.output}") + + elif args.algo == "pr": + if args.matrix: + A = read_mtx_pr_spla(str(args.matrix), alpha=args.alpha) + else: + A = read_vectors_pr_spla(str(args.vectors), alpha=args.alpha) + p = pr_spla(A, args.alpha, args.eps) + with open(args.output, "w") as f_out: + idx, vals = p.to_lists() + for k in range(len(idx)): + f_out.write(f"{idx[k]} {vals[k]:.6f}\n") + print(f"Result saved to {args.output}") + + elif args.algo == "tc": + if args.matrix: + A = read_spla(str(args.matrix), dtype=INT) + else: + A = read_vectors_spla(str(args.vectors), dtype=INT) + triangles = tc_spla(A) + with open(args.output, "w") as f_out: + f_out.write(f"Triangles: {triangles}\n") + print(f"Result saved to {args.output}") + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/python/algorithms/pr_classic.py b/python/algorithms/pr_classic.py new file mode 100644 index 000000000..66823f499 --- /dev/null +++ b/python/algorithms/pr_classic.py @@ -0,0 +1,21 @@ +import math + +def pagerank_naive(Ai, Ax, alpha,eps): + N=len(Ai) + p_prev=[1.0/N]*N + error=eps+1.0 + while error>eps: + p=[0.0]*N + for i in range(N): + s=0.0 + for k in range(len(Ai[i])): + s+=Ax[i][k]*p_prev[Ai[i][k]] + p[i]=s+(1.0-alpha)/N + error=0.0 + for i in range(N): + diff=p[i]-p_prev[i] + error+=diff*diff + error=math.sqrt(error) + p_prev,p=p,p_prev + + return p_prev \ No newline at end of file diff --git a/python/algorithms/pr_spla.py b/python/algorithms/pr_spla.py new file mode 100644 index 000000000..94aa978f3 --- /dev/null +++ b/python/algorithms/pr_spla.py @@ -0,0 +1,26 @@ +import math +from pyspla import FLOAT, Matrix, Scalar, Vector + + +def pr_spla(A: Matrix, alpha: float, eps: float): + N = A.n_rows + mask = Vector(N, FLOAT) + addition = Vector(N, FLOAT).fill_with(Scalar(FLOAT, (1.0 - alpha) / N)) + p_prev = Vector(N, FLOAT).fill_with(Scalar(FLOAT, 1.0 / N)) + p_tmp = Vector(N, FLOAT) + p = Vector(N, FLOAT) + errors = Vector(N, FLOAT) + zero = Scalar(FLOAT, 0.0) + error2 = Scalar(FLOAT) + error = eps + 0.1 + while error > eps: + A.mxv(mask, p_prev, + op_mult=FLOAT.MULT, op_add=FLOAT.PLUS, op_select=FLOAT.ALWAYS, + init=zero, out=p_tmp) + p_tmp.eadd(FLOAT.PLUS, addition, out=p) + p.eadd(FLOAT.MINUS_POW2, p_prev, out=errors) + errors.reduce(FLOAT.PLUS, out=error2, init=zero) + error = math.sqrt(error2.get()) + p, p_prev = p_prev, p + + return p_prev \ No newline at end of file diff --git a/python/algorithms/result_cpu.pdf b/python/algorithms/result_cpu.pdf new file mode 100644 index 000000000..949260fb1 Binary files /dev/null and b/python/algorithms/result_cpu.pdf differ diff --git a/python/algorithms/result_gpu.pdf b/python/algorithms/result_gpu.pdf new file mode 100644 index 000000000..307cec892 Binary files /dev/null and b/python/algorithms/result_gpu.pdf differ diff --git a/python/algorithms/sssp_classic.py b/python/algorithms/sssp_classic.py new file mode 100644 index 000000000..e95b13329 --- /dev/null +++ b/python/algorithms/sssp_classic.py @@ -0,0 +1,13 @@ +import math +INF = math.inf +def sssp_naive(start, graph, n): + dist = [INF] * n + dist[start] = 0 + + for _ in range(n - 1): + for u in range(n): + if dist[u] != INF: + for v, w in graph[u]: + dist[v] = min(dist[v], dist[u] + w) + + return dist diff --git a/python/algorithms/sssp_spla.py b/python/algorithms/sssp_spla.py new file mode 100644 index 000000000..192ac726d --- /dev/null +++ b/python/algorithms/sssp_spla.py @@ -0,0 +1,42 @@ +import math +from pyspla import FLOAT, Matrix, Scalar, Vector +INF = math.inf + + +def make_vector(n): + v = Vector(n, FLOAT) + v.set_fill_value(Scalar(FLOAT, INF)) + return v + + +def sssp_spla(start, A, push_only=False, pull_only=False, push_pull=True, front_factor=0.05): + n = A.n_rows + if not (push_only or pull_only or push_pull): + push_only = True + + mask = Vector(n, FLOAT) + init_inf = Scalar(FLOAT, INF) + dist = make_vector(n) + dist.set(start, 0.0) + feedback = make_vector(n) + feedback.set(start, 0.0) + frontier = make_vector(n) + fs = feedback.count_mf().get() + + while fs > 0: + front_density = fs / n + is_push_better = front_density <= front_factor + + if push_only or (push_pull and is_push_better): + feedback.vxm(mask, A, + op_mult=FLOAT.PLUS, op_add=FLOAT.MIN, + op_select=FLOAT.ALWAYS, init=init_inf, out=frontier) + else: + A.mxv(mask, feedback, + op_mult=FLOAT.PLUS, op_add=FLOAT.MIN, + op_select=FLOAT.ALWAYS, init=init_inf, out=frontier) + + dist.eadd_fdb(frontier, FLOAT.MIN, feedback) + fs = feedback.count_mf().get() + + return dist \ No newline at end of file diff --git a/python/algorithms/tc_classic.py b/python/algorithms/tc_classic.py new file mode 100644 index 000000000..cf091695f --- /dev/null +++ b/python/algorithms/tc_classic.py @@ -0,0 +1,10 @@ +def tc_naive(graph, n): + triangles = 0 + + for i in range(n): + for j in graph[i]: + for k in graph[j]: + if k in graph[i]: + triangles += 1 + + return triangles // 6 diff --git a/python/algorithms/tc_spla.py b/python/algorithms/tc_spla.py new file mode 100644 index 000000000..55fbb042b --- /dev/null +++ b/python/algorithms/tc_spla.py @@ -0,0 +1,6 @@ +from pyspla import INT, Matrix, Scalar + +def tc_spla(A: Matrix): + zero = Scalar(INT, 0) + B = A.mxmT(A, A, op_mult=INT.MULT, op_add=INT.PLUS, op_select=INT.GTZERO, init=zero) + return B.reduce(op_reduce=INT.PLUS, init=zero).get() \ No newline at end of file diff --git a/python/algorithms/test_compare.py b/python/algorithms/test_compare.py new file mode 100644 index 000000000..9f8d6ff81 --- /dev/null +++ b/python/algorithms/test_compare.py @@ -0,0 +1,171 @@ +import warnings + +warnings.simplefilter("ignore", DeprecationWarning) + +import heapq +import math +import sys +import unittest +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).parent)) + +from pyspla import FLOAT, INT +from graph_spla import read_spla, read_mtx_pr_spla +from bfs_spla import bfs_spla +from bfs_classic import bfs_naive as bfs_classic +from sssp_spla import sssp_spla +from pr_spla import pr_spla +from pr_classic import pagerank_naive as pr_classic +from tc_spla import tc_spla +from tc_classic import tc_naive as tc_classic + + +DATA = Path(__file__).parent.parent.parent / "bench" / "dataset" +UNDIRECTED_GRAPH = DATA / "belgium_osm.mtx" +DIRECTED_GRAPH = DATA / "amazon-2008.mtx" +UPPER = DATA / "upper" / "belgium_osm.mtx" + + +def split(A): + i, j, v = A.to_lists() + return A.n_rows, list(zip(i, j)), list(v) + + +def to_list(v, n, fill): + i, x = v.to_lists() + out = [fill] * n + for k, val in zip(i, x): + out[k] = val + return out + + +def adj_plain(n, pairs): + adj = [[] for _ in range(n)] + for u, v in pairs: + adj[u].append(v) + return adj + + +def adj_weighted(n, pairs, weights): + adj = [[] for _ in range(n)] + for (u, v), w in zip(pairs, weights): + adj[u].append((v, w)) + return adj + + +def short_distance(n, adj, start): + dist = [math.inf] * n + dist[start] = 0.0 + + queue = [(0.0, start)] + while queue: + d, u = heapq.heappop(queue) + if d > dist[u]: + continue + for v, w in adj[u]: + new = d + w + if new < dist[v]: + dist[v] = new + heapq.heappush(queue, (new, v)) + + return dist + + +class TestBfs(unittest.TestCase): + def _check(self, mtx, **kwargs): + A = read_spla(str(mtx), INT) + n, pairs, _ = split(A) + + ref = bfs_classic(0, adj_plain(n, pairs), n) + ref = [(x + 1) if x is not None else 0 for x in ref] + + v, _, _ = bfs_spla(0, A, **kwargs) + got = to_list(v, n, 0) + + self.assertEqual(ref, got) + + def test_symmetric_push(self): + self._check(UNDIRECTED_GRAPH, push_only=True, pull_only=False, push_pull=False) + + def test_symmetric_pull(self): + self._check(UNDIRECTED_GRAPH, push_only=False, pull_only=True, push_pull=False) + + def test_symmetric_push_pull(self): + self._check(UNDIRECTED_GRAPH, push_only=False, pull_only=False, push_pull=True) + + def test_directed_push(self): + self._check(DIRECTED_GRAPH, push_only=True, pull_only=False, push_pull=False) + + +class TestSssp(unittest.TestCase): + def _check(self, mtx, **kwargs): + A = read_spla(str(mtx), FLOAT) + n, pairs, weights = split(A) + + ref = short_distance(n, adj_weighted(n, pairs, weights), 0) + + d = sssp_spla(0, A, **kwargs) + got = to_list(d, n, math.inf) + got[0] = 0.0 + + for i, (a, b) in enumerate(zip(ref, got)): + if math.isinf(a) and math.isinf(b): + continue + self.assertLessEqual(abs(a - b), 1e-4, f"v{i}: ref={a} got={b}") + + def test_symmetric_push(self): + self._check(UNDIRECTED_GRAPH, push_only=True, pull_only=False, push_pull=False) + + def test_symmetric_pull(self): + self._check(UNDIRECTED_GRAPH, push_only=False, pull_only=True, push_pull=False) + + def test_directed_push(self): + self._check(DIRECTED_GRAPH, push_only=True, pull_only=False, push_pull=False) + + +class TestPr(unittest.TestCase): + def _check(self, mtx): + A = read_mtx_pr_spla(str(mtx)) + n, pairs, weights = split(A) + + adj = [[] for _ in range(n)] + wgt = [[] for _ in range(n)] + for (u, v), w in zip(pairs, weights): + adj[u].append(v) + wgt[u].append(w) + + alpha, eps = 0.85, 1e-6 + ref = pr_classic(adj, wgt, alpha, eps) + + p = pr_spla(A, alpha, eps) + got = to_list(p, n, 0.0) + + max_diff = max(abs(a - b) for a, b in zip(ref, got)) + self.assertLess(max_diff, 1e-4, f"max diff {max_diff:.2e}") + + def test_symmetric(self): + self._check(UNDIRECTED_GRAPH) + + def test_directed(self): + self._check(DIRECTED_GRAPH) + + +class TestTc(unittest.TestCase): + def test_upper(self): + A = read_spla(str(UPPER), INT) + n, pairs, _ = split(A) + + adj = [[] for _ in range(n)] + for u, v in pairs: + adj[u].append(v) + adj[v].append(u) + + ref = tc_classic(adj, n) + got = tc_spla(A) + + self.assertEqual(ref, got) + + +if __name__ == "__main__": + unittest.main(verbosity=2) \ No newline at end of file diff --git a/python/pyspla/bridge.py b/python/pyspla/bridge.py index 808f6ff17..a9cdb5c0f 100644 --- a/python/pyspla/bridge.py +++ b/python/pyspla/bridge.py @@ -214,6 +214,12 @@ def load_library(lib_path): _spla.spla_Type_UINT.argtypes = [] _spla.spla_Type_FLOAT.restype = _object_t _spla.spla_Type_FLOAT.argtypes = [] + _spla.spla_Descriptor_make.restype = _status_t + _spla.spla_Descriptor_make.argtypes = [_p_object_t] + _spla.spla_Descriptor_set_struct_only.restype = _status_t + _spla.spla_Descriptor_set_struct_only.argtypes = [_object_t, _int] + _spla.spla_Descriptor_set_early_exit.restype = _status_t + _spla.spla_Descriptor_set_early_exit.argtypes = [_object_t, _int] _spla.spla_OpUnary_IDENTITY_INT.restype = _object_t _spla.spla_OpUnary_IDENTITY_UINT.restype = _object_t @@ -494,6 +500,7 @@ def load_library(lib_path): _spla.spla_Vector_build.restype = _status_t _spla.spla_Vector_read.restype = _status_t _spla.spla_Vector_clear.restype = _status_t + _spla.spla_Vector_fill_with.restype = _status_t _spla.spla_Vector_make.argtypes = [_p_object_t, _uint, _object_t] _spla.spla_Vector_set_format.argtypes = [_object_t, ctypes.c_int] @@ -508,6 +515,7 @@ def load_library(lib_path): _spla.spla_Vector_build.argtypes = [_object_t, _object_t, _object_t] _spla.spla_Vector_read.argtypes = [_object_t, _p_object_t, _p_object_t] _spla.spla_Vector_clear.argtypes = [_object_t] + _spla.spla_Vector_fill_with.argtypes = [_object_t, _object_t] _spla.spla_Matrix_make.restype = _status_t _spla.spla_Matrix_set_format.restype = _status_t diff --git a/python/pyspla/descriptor.py b/python/pyspla/descriptor.py index 2e8d014d6..1c5c0fc1a 100644 --- a/python/pyspla/descriptor.py +++ b/python/pyspla/descriptor.py @@ -48,4 +48,16 @@ def __init__(self, label=None, hnd=None): MemView native void* handle to a C counterpart. """ + if hnd is None: + hnd = ctypes.c_void_p(0) + check(backend().spla_Descriptor_make(ctypes.byref(hnd))) + super().__init__(label, hnd) + + def set_struct_only(self, value: bool): + check(backend().spla_Descriptor_set_struct_only(self.hnd, ctypes.c_int(int(value)))) + return self + + def set_early_exit(self, value: bool): + check(backend().spla_Descriptor_set_early_exit(self.hnd, ctypes.c_int(int(value)))) + return self diff --git a/python/pyspla/vector.py b/python/pyspla/vector.py index 61f379ae7..cf796e850 100644 --- a/python/pyspla/vector.py +++ b/python/pyspla/vector.py @@ -246,7 +246,11 @@ def read(self): keys_view_hnd = ctypes.c_void_p(0) values_view_hnd = ctypes.c_void_p(0) - check(backend().spla_Vector_read(self.hnd, ctypes.byref(keys_view_hnd), ctypes.byref(values_view_hnd))) + check( + backend().spla_Vector_read( + self.hnd, ctypes.byref(keys_view_hnd), ctypes.byref(values_view_hnd) + ) + ) return MemView(hnd=keys_view_hnd, owner=self), MemView(hnd=values_view_hnd, owner=self) def clear(self): @@ -275,8 +279,16 @@ def to_lists(self): buffer_I = (UINT._c_type * count)() buffer_V = (self._dtype._c_type * count)() - check(backend().spla_MemView_read(I.hnd, ctypes.c_size_t(0), ctypes.sizeof(buffer_I), buffer_I)) - check(backend().spla_MemView_read(V.hnd, ctypes.c_size_t(0), ctypes.sizeof(buffer_V), buffer_V)) + check( + backend().spla_MemView_read( + I.hnd, ctypes.c_size_t(0), ctypes.sizeof(buffer_I), buffer_I + ) + ) + check( + backend().spla_MemView_read( + V.hnd, ctypes.c_size_t(0), ctypes.sizeof(buffer_V), buffer_V + ) + ) return list(buffer_I), list(buffer_V) @@ -331,7 +343,9 @@ def to_string(self, format_string="{:>%s}", width=2, precision=2, skip_value=0): value = self.get(row) value = value if value != skip_value else "." result += format_string.format(row) + "|" - result += format_string.format(self.dtype.format_value(value, width, precision)).rstrip() + result += format_string.format( + self.dtype.format_value(value, width, precision) + ).rstrip() if row < self.n_rows - 1: result += "\n" @@ -543,9 +557,20 @@ def vxm(self, mask, M, op_mult, op_add, op_select, out=None, init=None, desc=Non assert mask.n_rows == M.n_cols assert M.n_rows == self.n_rows - check(backend().spla_Exec_vxm_masked(out.hnd, mask.hnd, self.hnd, M.hnd, - op_mult.hnd, op_add.hnd, op_select.hnd, - init.hnd, self._get_desc(desc), self._get_task(None))) + check( + backend().spla_Exec_vxm_masked( + out.hnd, + mask.hnd, + self.hnd, + M.hnd, + op_mult.hnd, + op_add.hnd, + op_select.hnd, + init.hnd, + self._get_desc(desc), + self._get_task(None), + ) + ) return out @@ -588,8 +613,11 @@ def eadd(self, op_add, v, out=None, desc=None): assert out.dtype == out.dtype assert op_add - check(backend().spla_Exec_v_eadd(out.hnd, self.hnd, v.hnd, op_add.hnd, - self._get_desc(desc), self._get_task(None))) + check( + backend().spla_Exec_v_eadd( + out.hnd, self.hnd, v.hnd, op_add.hnd, self._get_desc(desc), self._get_task(None) + ) + ) return out @@ -632,8 +660,11 @@ def emult(self, op_mult, v, out=None, desc=None): assert out.dtype == out.dtype assert op_mult - check(backend().spla_Exec_v_emult(out.hnd, self.hnd, v.hnd, op_mult.hnd, - self._get_desc(desc), self._get_task(None))) + check( + backend().spla_Exec_v_emult( + out.hnd, self.hnd, v.hnd, op_mult.hnd, self._get_desc(desc), self._get_task(None) + ) + ) return out @@ -677,8 +708,17 @@ def assign(self, mask, value, op_assign, op_select, desc=None): assert op_assign assert op_select - check(backend().spla_Exec_v_assign_masked(self.hnd, mask.hnd, value.hnd, op_assign.hnd, op_select.hnd, - self._get_desc(desc), self._get_task(None))) + check( + backend().spla_Exec_v_assign_masked( + self.hnd, + mask.hnd, + value.hnd, + op_assign.hnd, + op_select.hnd, + self._get_desc(desc), + self._get_task(None), + ) + ) return self @@ -715,8 +755,11 @@ def map(self, op_map, out=None, desc=None): assert out.n_rows == self.n_rows assert out.dtype == self.dtype - check(backend().spla_Exec_v_map(out.hnd, self.hnd, op_map.hnd, - self._get_desc(desc), self._get_task(None))) + check( + backend().spla_Exec_v_map( + out.hnd, self.hnd, op_map.hnd, self._get_desc(desc), self._get_task(None) + ) + ) return out @@ -759,8 +802,45 @@ def reduce(self, op_reduce, out=None, init=None, desc=None): assert out.dtype == self.dtype assert init.dtype == self.dtype - check(backend().spla_Exec_v_reduce(out.hnd, init.hnd, self.hnd, op_reduce.hnd, - self._get_desc(desc), self._get_task(None))) + check( + backend().spla_Exec_v_reduce( + out.hnd, + init.hnd, + self.hnd, + op_reduce.hnd, + self._get_desc(desc), + self._get_task(None), + ) + ) + + return out + + def eadd_fdb(self, v, op_add, fdb, desc=None): + assert v + assert fdb + assert op_add + assert v.n_rows == self.n_rows + assert fdb.n_rows == self.n_rows + assert v.dtype == self.dtype + assert fdb.dtype == self.dtype + + check( + backend().spla_Exec_v_eadd_fdb( + self.hnd, v.hnd, fdb.hnd, op_add.hnd, self._get_desc(desc), self._get_task(None) + ) + ) + + return self + + def count_mf(self, out=None, desc=None): + if out is None: + out = Scalar(INT) + + check( + backend().spla_Exec_v_count_mf( + out.hnd, self.hnd, self._get_desc(desc), self._get_task(None) + ) + ) return out @@ -776,3 +856,11 @@ def _get_desc(self, desc: Descriptor): def _get_task(self, task): return ctypes.POINTER(ctypes.c_void_p)() + + def set_fill_value(self, scalar): + check(backend().spla_Vector_set_fill_value(self.hnd, scalar.hnd)) + return self + + def fill_with(self, scalar): + check(backend().spla_Vector_fill_with(self.hnd, scalar.hnd)) + return self diff --git a/src/binding/c_vector.cpp b/src/binding/c_vector.cpp index 25619d1bc..a263e2edf 100644 --- a/src/binding/c_vector.cpp +++ b/src/binding/c_vector.cpp @@ -62,6 +62,9 @@ spla_Status spla_Vector_get_float(spla_Vector v, spla_uint row_id, float* value) spla_Status spla_Vector_build(spla_Vector v, spla_MemView keys, spla_MemView values) { return to_c_status(as_ptr(v)->build(as_ref(keys), as_ref(values))); } +spla_Status spla_Vector_fill_with(spla_Vector v, spla_Scalar value) { + return to_c_status(as_ptr(v)->fill_with(as_ref(value))); +} spla_Status spla_Vector_read(spla_Vector v, spla_MemView* keys, spla_MemView* values) { spla::ref_ptr out_keys; spla::ref_ptr out_values; @@ -75,4 +78,4 @@ spla_Status spla_Vector_read(spla_Vector v, spla_MemView* keys, spla_MemView* va } spla_Status spla_Vector_clear(spla_Vector v) { return to_c_status(as_ptr(v)->clear()); -} \ No newline at end of file +} diff --git a/src/io.cpp b/src/io.cpp index 9a962087d..3989bf969 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -63,12 +63,20 @@ namespace spla { std::size_t n_lines = 0; std::size_t n_sort = 0; + std::string first_line; std::string line; while (std::getline(file, line)) { + if (first_line.empty() && line.rfind("%%MatrixMarket", 0) == 0) { + first_line = line; + } if (line[0] != '%') break; n_lines++; } + bool file_is_symmetric = first_line.find("symmetric") != std::string::npos; + bool file_is_pattern = first_line.find("pattern") != std::string::npos; + bool effective_undirected = make_undirected || file_is_symmetric; + std::size_t nnz; std::stringstream header(line); header >> m_n_rows >> m_n_cols >> nnz; @@ -79,7 +87,11 @@ namespace spla { std::cout << " Data: " << nnz << " directed edges" << std::endl; if (remove_loops) std::cout << " Opt: remove self-loops" << std::endl; if (offset_indices) std::cout << " Opt: offset indices by -1" << std::endl; - if (make_undirected) std::cout << " Opt: double edges" << std::endl; + if (effective_undirected) { + std::cout << " Opt: double edges"; + if (file_is_symmetric) std::cout << " (auto from header)"; + std::cout << std::endl; + } std::cout << " Reading data: "; // optimized reading by sliding window @@ -89,16 +101,16 @@ namespace spla { char buffer[BUFFER_CAPACITY + 1]; // read data - std::size_t to_count = 0; - std::size_t to_read = nnz; - std::size_t to_preallocate = to_read * (make_undirected ? 2 : 1); - std::vector Ai; - std::vector Aj; - + std::size_t to_count = 0; + std::size_t to_read = nnz; + std::size_t to_preallocate = to_read * (effective_undirected ? 2 : 1); + std::vector Ai; + std::vector Aj; + std::vector Aw; // preallocate to avoid copy Ai.reserve(to_preallocate); Aj.reserve(to_preallocate); - + Aw.reserve(to_preallocate); float job_done = 0.0f; float job_total = 35.0f; @@ -149,9 +161,15 @@ namespace spla { } } - char* end = nullptr; - auto i = uint(std::strtoll(buffer + buffer_offset, &end, 10)); - auto j = uint(std::strtoll(end, &end, 10)); + char* end = nullptr; + auto i = uint(std::strtoll(buffer + buffer_offset, &end, 10)); + auto j = uint(std::strtoll(end, &end, 10)); + float w = 1.0f; + if (!file_is_pattern && end < buffer + line_end) { + char* val_end = nullptr; + float parsed = std::strtof(end, &val_end); + if (val_end != end && val_end <= buffer + line_end) w = parsed; + } buffer_offset = line_end + 1; assert(i > 0 && j > 0); @@ -163,54 +181,65 @@ namespace spla { i -= 1; j -= 1; } - if (make_undirected) { + if (effective_undirected) { Ai.push_back(j); Aj.push_back(i); + Aw.push_back(w); } Ai.push_back(i); Aj.push_back(j); + Aw.push_back(w); } t.lap_end();// parsing - std::vector sorted; + struct Entry { + uint i; + uint j; + float w; + }; + + std::vector sorted; { sorted.reserve(Ai.size()); n_sort = Ai.size(); for (std::size_t k = 0; k < Ai.size(); k++) { - std::uint64_t entry = 0; - entry |= std::uint64_t(Ai[k]) << 32u; - entry |= std::uint64_t(Aj[k]) << 0u; - sorted.push_back(entry); + sorted.push_back({Ai[k], Aj[k], Aw[k]}); } Ai.clear(); Aj.clear(); + Aw.clear(); - std::sort(sorted.begin(), sorted.end()); + std::sort(sorted.begin(), sorted.end(), + [](const Entry& a, const Entry& b) { + return std::tie(a.i, a.j) < std::tie(b.i, b.j); + }); } t.lap_end();// sorting - std::vector reduced_Ai; - std::vector reduced_Aj; + std::vector reduced_Ai; + std::vector reduced_Aj; + std::vector reduced_Aw; { reduced_Ai.reserve(sorted.size()); reduced_Aj.reserve(sorted.size()); - - std::uint64_t entry_prev = 0xffffffffffffffff; - for (std::uint64_t entry : sorted) { - if (entry_prev != entry) { - uint i = uint((entry >> 32u) & 0xffffffff); - uint j = uint((entry >> 0u) & 0xffffffff); - reduced_Ai.push_back(i); - reduced_Aj.push_back(j); + reduced_Aw.reserve(sorted.size()); + + Entry prev{0xffffffffu, 0xffffffffu, 0.0f}; + for (const auto& e : sorted) { + if (e.i != prev.i || e.j != prev.j) { + reduced_Ai.push_back(e.i); + reduced_Aj.push_back(e.j); + reduced_Aw.push_back(e.w); } - entry_prev = entry; + prev = e; } m_n_values = reduced_Ai.size(); m_Ai = std::move(reduced_Ai); m_Aj = std::move(reduced_Aj); + m_Aw = std::move(reduced_Aw); } t.lap_end();// reducing @@ -367,7 +396,9 @@ namespace spla { const std::vector& MtxLoader::get_Aj() const { return m_Aj; } - + const std::vector& MtxLoader::get_Aw() const { + return m_Aw; + } uint MtxLoader::get_n_rows() const { return m_n_rows; } @@ -378,4 +409,4 @@ namespace spla { return m_n_values; } -}// namespace spla \ No newline at end of file +}// namespace spla diff --git a/src/opencl/cl_v_count_mf.hpp b/src/opencl/cl_v_count_mf.hpp index b57d818e0..544ccfcab 100644 --- a/src/opencl/cl_v_count_mf.hpp +++ b/src/opencl/cl_v_count_mf.hpp @@ -73,9 +73,10 @@ namespace spla { private: Status execute_sp(const DispatchContext& ctx) { - auto t = ctx.task.template cast_safe(); - ref_ptr> v = t->v.template cast_safe>(); - CLCooVec* dec_v = v->template get>(); + auto t = ctx.task.template cast_safe(); + ref_ptr> v = t->v.template cast_safe>(); + v->validate_rw(FormatVector::AccCoo); + CLCooVec* dec_v = v->template get>(); t->r->set_uint(dec_v->values); @@ -83,9 +84,10 @@ namespace spla { } Status execute_dn(const DispatchContext& ctx) { - auto t = ctx.task.template cast_safe(); - ref_ptr> v = t->v.template cast_safe>(); - CLDenseVec* dec_v = v->template get>(); + auto t = ctx.task.template cast_safe(); + ref_ptr> v = t->v.template cast_safe>(); + v->validate_rw(FormatVector::AccDense); + CLDenseVec* dec_v = v->template get>(); std::shared_ptr program; if (!ensure_kernel(program)) return Status::CompilationError;