From 391ce3d5b9d2b536a5f023b76e5ca13bd906e1fa Mon Sep 17 00:00:00 2001 From: AlexandruChiri <63013513+AlexandruChiri@users.noreply.github.com> Date: Thu, 21 Dec 2023 19:50:38 +0200 Subject: [PATCH] Added sparse benchmarks I have added a few sparse benchmarks to NPBench: - Banded Matrix-Matrix-Transposed multiplication (A * B * A^T). NumPy already had an implementation for banded matrix multiplication, but it required the data to be in dense format (2D numpy arrays). - Conjugate Gradient Method for CSR matrices. - BiConjugate Gradient Method for CSR matrices. - BiConjugate Stabilized Gradient Method for CSR matrices. - Minimum Residual Method for CSR matrices. - Generalized Minimum Residual Method for CSR matrices. - CSR Matrix-Matrix multiplication. --- bench_info/banded_mmt.json | 26 +++++++ bench_info/sp_bicg.json | 26 +++++++ bench_info/sp_bicgstab.json | 26 +++++++ bench_info/sp_cg.json | 26 +++++++ bench_info/sp_gmres.json | 26 +++++++ bench_info/sp_minres.json | 26 +++++++ bench_info/spmm.json | 26 +++++++ .../sparse/banded_mmt/banded_mmt.py | 36 ++++++++++ .../sparse/banded_mmt/banded_mmt_numpy.py | 70 ++++++++++++++++++ .../benchmarks/sparse/solvers/bicg/bicg.py | 51 +++++++++++++ .../sparse/solvers/bicg/bicg_numpy.py | 25 +++++++ .../sparse/solvers/bicgstab/bicgstab.py | 69 ++++++++++++++++++ .../sparse/solvers/bicgstab/bicgstab_numpy.py | 25 +++++++ npbench/benchmarks/sparse/solvers/cg/cg.py | 71 ++++++++++++++++++ .../benchmarks/sparse/solvers/cg/cg_numpy.py | 18 +++++ .../benchmarks/sparse/solvers/gmres/gmres.py | 52 ++++++++++++++ .../sparse/solvers/gmres/gmres_numpy.py | 36 ++++++++++ .../sparse/solvers/minres/minres.py | 72 +++++++++++++++++++ .../sparse/solvers/minres/minres_numpy.py | 19 +++++ npbench/benchmarks/sparse/spmm/spmm.py | 14 ++++ npbench/benchmarks/sparse/spmm/spmm_numpy.py | 5 ++ 21 files changed, 745 insertions(+) create mode 100644 bench_info/banded_mmt.json create mode 100644 bench_info/sp_bicg.json create mode 100644 bench_info/sp_bicgstab.json create mode 100644 bench_info/sp_cg.json create mode 100644 bench_info/sp_gmres.json create mode 100644 bench_info/sp_minres.json create mode 100644 bench_info/spmm.json create mode 100644 npbench/benchmarks/sparse/banded_mmt/banded_mmt.py create mode 100644 npbench/benchmarks/sparse/banded_mmt/banded_mmt_numpy.py create mode 100644 npbench/benchmarks/sparse/solvers/bicg/bicg.py create mode 100644 npbench/benchmarks/sparse/solvers/bicg/bicg_numpy.py create mode 100644 npbench/benchmarks/sparse/solvers/bicgstab/bicgstab.py create mode 100644 npbench/benchmarks/sparse/solvers/bicgstab/bicgstab_numpy.py create mode 100644 npbench/benchmarks/sparse/solvers/cg/cg.py create mode 100644 npbench/benchmarks/sparse/solvers/cg/cg_numpy.py create mode 100644 npbench/benchmarks/sparse/solvers/gmres/gmres.py create mode 100644 npbench/benchmarks/sparse/solvers/gmres/gmres_numpy.py create mode 100644 npbench/benchmarks/sparse/solvers/minres/minres.py create mode 100644 npbench/benchmarks/sparse/solvers/minres/minres_numpy.py create mode 100644 npbench/benchmarks/sparse/spmm/spmm.py create mode 100644 npbench/benchmarks/sparse/spmm/spmm_numpy.py diff --git a/bench_info/banded_mmt.json b/bench_info/banded_mmt.json new file mode 100644 index 0000000..7b88ab0 --- /dev/null +++ b/bench_info/banded_mmt.json @@ -0,0 +1,26 @@ +{ + "benchmark": { + "name": "Banded Matrix-Matrix-Transposed Multiplication", + "short_name": "banded_mmt", + "relative_path": "sparse/banded_mmt", + "module_name": "banded_mmt", + "func_name": "banded_mmt", + "kind": "microapp", + "domain": "Other", + "dwarf": "sparse_linear_algebra", + "parameters": { + "S": { "N": 100, "a_lbound": 3, "a_ubound": 4, "b_lbound": 5, "b_ubound": 2 }, + "M": { "N": 1000, "a_lbound": 11, "a_ubound": 9, "b_lbound": 6, "b_ubound": 13 }, + "L": { "N": 10000, "a_lbound": 15, "a_ubound": 12, "b_lbound": 13, "b_ubound": 11 }, + "paper": { "N": 20000, "a_lbound": 43, "a_ubound": 31, "b_lbound": 15, "b_ubound": 50 } + }, + "init": { + "func_name": "initialize", + "input_args": ["N", "a_lbound", "a_ubound", "b_lbound", "b_ubound"], + "output_args": ["A", "B"] + }, + "input_args": ["A", "a_lbound", "a_ubound", "B", "b_lbound", "b_ubound"], + "array_args": ["A", "B"], + "output_args": [] + } +} diff --git a/bench_info/sp_bicg.json b/bench_info/sp_bicg.json new file mode 100644 index 0000000..96d69fd --- /dev/null +++ b/bench_info/sp_bicg.json @@ -0,0 +1,26 @@ +{ + "benchmark": { + "name": "Sparse Biconjugate Gradient", + "short_name": "sp_bicg", + "relative_path": "sparse/solvers/bicg", + "module_name": "bicg", + "func_name": "bicg", + "kind": "microapp", + "domain": "Sparse solver", + "dwarf": "sparse_linear_algebra", + "parameters": { + "S": { "N": 2048, "nnz": 409600 }, + "M": { "N": 4096, "nnz": 1638400 }, + "L": { "N": 8192, "nnz": 6553600 }, + "paper": { "N": 8192, "nnz": 6553600 } + }, + "init": { + "func_name": "initialize", + "input_args": ["N", "nnz"], + "output_args": ["A", "x", "b"] + }, + "input_args": ["A", "x", "b"], + "array_args": ["x", "b"], + "output_args": [] + } +} \ No newline at end of file diff --git a/bench_info/sp_bicgstab.json b/bench_info/sp_bicgstab.json new file mode 100644 index 0000000..e9a8d5f --- /dev/null +++ b/bench_info/sp_bicgstab.json @@ -0,0 +1,26 @@ +{ + "benchmark": { + "name": "Sparse Biconjugate Stabilized Gradient", + "short_name": "sp_bicgstab", + "relative_path": "sparse/solvers/bicgstab", + "module_name": "bicgstab", + "func_name": "bicgstab", + "kind": "microapp", + "domain": "Sparse solver", + "dwarf": "sparse_linear_algebra", + "parameters": { + "S": { "N": 2048, "nnz": 409600 }, + "M": { "N": 4096, "nnz": 1638400 }, + "L": { "N": 8192, "nnz": 6553600 }, + "paper": { "N": 8192, "nnz": 6553600 } + }, + "init": { + "func_name": "initialize", + "input_args": ["N", "nnz"], + "output_args": ["A", "x", "b"] + }, + "input_args": ["A", "x", "b"], + "array_args": ["x", "b"], + "output_args": [] + } +} \ No newline at end of file diff --git a/bench_info/sp_cg.json b/bench_info/sp_cg.json new file mode 100644 index 0000000..5319255 --- /dev/null +++ b/bench_info/sp_cg.json @@ -0,0 +1,26 @@ +{ + "benchmark": { + "name": "Sparse Conjugate Gradient", + "short_name": "sp_cg", + "relative_path": "sparse/solvers/cg", + "module_name": "cg", + "func_name": "cg", + "kind": "microapp", + "domain": "Sparse solver", + "dwarf": "sparse_linear_algebra", + "parameters": { + "S": { "N": 2048, "nnz": 409600 }, + "M": { "N": 4096, "nnz": 1638400 }, + "L": { "N": 8192, "nnz": 6553600 }, + "paper": { "N": 8192, "nnz": 6553600 } + }, + "init": { + "func_name": "initialize", + "input_args": ["N", "nnz"], + "output_args": ["A", "x", "b"] + }, + "input_args": ["A", "x", "b"], + "array_args": ["x", "b"], + "output_args": [] + } +} \ No newline at end of file diff --git a/bench_info/sp_gmres.json b/bench_info/sp_gmres.json new file mode 100644 index 0000000..574c2a6 --- /dev/null +++ b/bench_info/sp_gmres.json @@ -0,0 +1,26 @@ +{ + "benchmark": { + "name": "Sparse General Minimal Residual", + "short_name": "sp_gmres", + "relative_path": "sparse/solvers/gmres", + "module_name": "gmres", + "func_name": "hand_gmres", + "kind": "microapp", + "domain": "Sparse solver", + "dwarf": "sparse_linear_algebra", + "parameters": { + "S": { "N": 2048, "nnz": 409600 }, + "M": { "N": 4096, "nnz": 1638400 }, + "L": { "N": 8192, "nnz": 6553600 }, + "paper": { "N": 8192, "nnz": 6553600 } + }, + "init": { + "func_name": "initialize", + "input_args": ["N", "nnz"], + "output_args": ["A", "x", "b"] + }, + "input_args": ["A", "x", "b"], + "array_args": ["x", "b"], + "output_args": [] + } +} \ No newline at end of file diff --git a/bench_info/sp_minres.json b/bench_info/sp_minres.json new file mode 100644 index 0000000..5cdc91f --- /dev/null +++ b/bench_info/sp_minres.json @@ -0,0 +1,26 @@ +{ + "benchmark": { + "name": "Sparse Minimal Residual", + "short_name": "sp_minres", + "relative_path": "sparse/solvers/minres", + "module_name": "minres", + "func_name": "hand_minres", + "kind": "microapp", + "domain": "Sparse solver", + "dwarf": "sparse_linear_algebra", + "parameters": { + "S": { "N": 2048, "nnz": 409600 }, + "M": { "N": 4096, "nnz": 1638400 }, + "L": { "N": 8192, "nnz": 6553600 }, + "paper": { "N": 8192, "nnz": 6553600 } + }, + "init": { + "func_name": "initialize", + "input_args": ["N", "nnz"], + "output_args": ["A", "x", "b"] + }, + "input_args": ["A", "x", "b"], + "array_args": ["x", "b"], + "output_args": [] + } +} \ No newline at end of file diff --git a/bench_info/spmm.json b/bench_info/spmm.json new file mode 100644 index 0000000..b35b2a9 --- /dev/null +++ b/bench_info/spmm.json @@ -0,0 +1,26 @@ +{ + "benchmark": { + "name": "Sparse Matrix-Matrix Multiplication", + "short_name": "spmm", + "relative_path": "sparse/spmm", + "module_name": "spmm", + "func_name": "spmm", + "kind": "microapp", + "domain": "LinAlg", + "dwarf": "sparse_linear_algebra", + "parameters": { + "S": { "NI": 4096, "NJ": 4096, "NK": 4096, "nnz_A": 163840, "nnz_B": 163840 }, + "M": { "NI": 8192, "NJ": 8192, "NK": 8192, "nnz_A": 655360, "nnz_B": 655360 }, + "L": { "NI": 16384, "NJ": 16384, "NK": 16384, "nnz_A": 2621440, "nnz_B": 2621440 }, + "paper": { "NI": 4096, "NJ": 4096, "NK": 4096, "nnz_A": 1638400, "nnz_B": 1638400 } + }, + "init": { + "func_name": "initialize", + "input_args": ["NI", "NJ", "NK", "nnz_A", "nnz_B"], + "output_args": ["alpha", "beta", "C", "A", "B"] + }, + "input_args": ["alpha", "beta", "C", "A", "B"], + "array_args": ["C"], + "output_args": [] + } +} diff --git a/npbench/benchmarks/sparse/banded_mmt/banded_mmt.py b/npbench/benchmarks/sparse/banded_mmt/banded_mmt.py new file mode 100644 index 0000000..20b6fa0 --- /dev/null +++ b/npbench/benchmarks/sparse/banded_mmt/banded_mmt.py @@ -0,0 +1,36 @@ +# Copyright 2023 University Politehnica of Bucharest and the NPBench authors. All rights reserved. + +import numpy as np + +# Function which stores and returns a banded square matrix in +# the compressed form with random elements +def generate_banded(lbound : int, ubound : int, size : int, dtype : type = np.float64): + # Allocates the matrix and initialises its elements with 0 + ret = np.zeros([size, min(lbound + ubound + 1, size)], dtype) + for i in range(0, size): + # Calculates the position of the first non-zero element on the current line + start = max(i - lbound, 0) + # Calculates the position of the first zero element after all the + # non-zero elements within the given bounds + stop = min(size, i + ubound + 1) + # Stores the non-zero elements from the current line + ret[i][0 : stop - start] = np.random.rand(stop - start).astype(dtype) + return ret + +# Function which stores and returns a banded square matrix in +# the compressed form with random elements +def generate_banded_scipy(lbound : int, ubound : int, size : int, dtype : type = np.float64): + # A = generate_banded(lbound, ubound, size, dtype=dtype) + diag_indexes = np.arange(-lbound, ubound + 1) + diags = np.empty(lbound + ubound + 1, dtype=object) + for i in range(diag_indexes.size): + diags[i] = np.random.rand(size - abs(diag_indexes[i])).astype(dtype) + import scipy.sparse as sp + return sp.diags(diags, diag_indexes, shape=(size, size)) + + +def initialize(N : int, a_lbound : int, a_ubound : int, b_lbound : int, b_ubound : int, dtype : type = np.float64): + np.random.seed(42) + A = generate_banded(a_lbound, a_ubound, N, dtype = dtype) + B = generate_banded(b_lbound, b_ubound, N, dtype = dtype) + return A, B diff --git a/npbench/benchmarks/sparse/banded_mmt/banded_mmt_numpy.py b/npbench/benchmarks/sparse/banded_mmt/banded_mmt_numpy.py new file mode 100644 index 0000000..eb8ecd1 --- /dev/null +++ b/npbench/benchmarks/sparse/banded_mmt/banded_mmt_numpy.py @@ -0,0 +1,70 @@ +# Bounded Matrix_1 * Matrix_2 * Transposed_1 +import numpy as np + +# Returns the transposed of a banded square matrix +def transposed(A : np.ndarray, lbound : int, ubound : int): + size = A.shape[0] + width = A.shape[1] + ret = np.zeros([size, width]) + # Stores the indexes of the first non-zero element for + # each line in the transposed matrix + ret_start = np.array(list(map(lambda i: max(i - ubound, 0), range(0, size)))) + # for i in range(0, size): + # ret_start[i] = max(i - ubound, 0) + for i in range(0, size): + start = max(i - lbound, 0) + stop = min(size, i + ubound + 1) + for j in range(0, stop - start): + # Get the column of the current element in the coresponding dense matrix + dense_j = j + start + ret[dense_j][i - ret_start[dense_j]] = A[i][j] + return ret + + +# Returns the product of a banded square matrix with the transposed of +# the another banded square matrix +# Also returns the bounds of the resulted matrix +def banded_dgemt(A : np.ndarray, a_lbound : int, a_ubound : int, B : np.ndarray, b_lbound : int, b_ubound : int): + if not A.shape[0] == B.shape[0]: + print(f"Cannot multiply square matrixes with different sizes {A.shape[0]} {B.shape[0]}") + return None + size = A.shape[0] + # A bound cannot be less than -size or bigger than size - 1 + lbound = max(-size, min(a_lbound + b_ubound, size - 1)) + ubound = max(-size, min(a_ubound + b_lbound, size - 1)) + ret = np.zeros([size, min(size, 1 + lbound + ubound)]) + for i in range(0, size): + start = max(i - lbound, 0) + stop = min(size, i + ubound + 1) + a_start = max(0, i - a_lbound) + a_cnt = 1 + min(size - i - 1, a_ubound) + min(i, a_lbound) + a_stop = min(size, i + a_ubound + 1) + for j in range(start, stop): + b_start = max(0, j - b_lbound) + b_cnt = 1 + min(size - j - 1, b_ubound) + min(j, b_lbound) + b_stop = min(size, i + b_ubound + 1) + acc = A.dtype.type(0) + offset_a = 0 + offset_b = 0 + if a_start >= b_start: + offset_a = a_start - b_start + else: + offset_b = b_start - a_start + interval = min(a_cnt - offset_b, b_cnt - offset_a) + ret[i][j - start] = A[i][offset_b : offset_b + interval] @ B[j][offset_a : offset_a + interval] + return ret, lbound, ubound + +# Returns the product of a banded square matrix another banded square matrix +# Also returns the bounds of the resulted matrix +def banded_dgemm(A : np.ndarray, a_lbound : int, a_ubound : int, B : np.ndarray, b_lbound : int, b_ubound : int): + if not A.shape[0] == B.shape[0]: + print(f"Cannot multiply square matrixes with different sizes {A.shape[0]} {B.shape[0]}") + return None + return banded_dgemt(A, a_lbound, a_ubound, transposed(B, b_lbound, b_ubound), b_ubound, b_lbound) + +# Returns the result of the A @ B @ A^T and its bounds +def banded_mmt(A : np.ndarray, a_lbound : int, a_ubound : int, + B : np.ndarray, b_lbound : int, b_ubound : int): + ret, lbound, ubound = banded_dgemm(A, a_lbound, a_ubound, B, b_lbound, b_ubound) + ret, lbound, ubound = banded_dgemt(ret, lbound, ubound, A, a_lbound, a_ubound) + return ret, lbound, ubound diff --git a/npbench/benchmarks/sparse/solvers/bicg/bicg.py b/npbench/benchmarks/sparse/solvers/bicg/bicg.py new file mode 100644 index 0000000..795c9db --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/bicg/bicg.py @@ -0,0 +1,51 @@ +# Copyright 2023 University Politehnica of Bucharest and the NPBench authors. All rights reserved. + +import numpy as np + +# Generate required elements for a matrix with at least one non-zero term +# in the determinant +def required_elems(n): + rows = np.random.choice(n, n, replace=False) + cols = np.random.choice(n, n, replace=False) + vals = np.array(np.random.rand(n) * 10 - 5, dtype=np.float64) + return list(rows), list(cols), list(vals) + +# Generate more elements to make sure there are nnz or nnz + 1 +def ensure_nnz(n, nnz, rows, cols, vals): + nnz -= len(rows) + if nnz <= 0: + return + coords = set(zip(rows, cols)) + # Generate at least as many random values as they are required + new_vals = np.array(np.random.rand(nnz) * 10 - 5, dtype=np.float64) + i = 0 + while nnz > 0: + x, y = np.random.choice(n, 2) + stop = y + # The following loop changes x and y untill they are a new pair of coordinates + while (x, y) in coords: + y += 1 + if y == n: + y = 0 + if y == stop: + x += 1 + if x == n: + x = 0 + generated_pair = (x, y) + coords.add(generated_pair) + rows.append(generated_pair[0]) + cols.append(generated_pair[1]) + vals.append(new_vals[i]) + nnz -= 1 + i += 1 + +def initialize(n : int, nnz : int, dtype=np.float64): + np.random.seed(42) + rows, cols, vals = required_elems(n) + ensure_nnz(n, nnz, rows, cols, vals) + from scipy.sparse import coo_matrix + A = coo_matrix((vals, (rows, cols)), shape=(n, n)).asformat('csr') + b = A @ np.random.rand(n) + # Generate starting solution + x = np.random.rand(n) + return A, x, b diff --git a/npbench/benchmarks/sparse/solvers/bicg/bicg_numpy.py b/npbench/benchmarks/sparse/solvers/bicg/bicg_numpy.py new file mode 100644 index 0000000..2bfc0f3 --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/bicg/bicg_numpy.py @@ -0,0 +1,25 @@ +import numpy as np + +# Solves A @ x = b where A is a Compressed Sparse Row matrix using the Biconjugate Gradient method +def bicg(A, b, x, max_iter=100, tol=np.float64(1e-6)): + n = A.shape[0] + r = b - A @ x + r_tilde = np.copy(r) + p = np.copy(r) + p_tilde = np.copy(r_tilde) + x = np.copy(x) + rho = r_tilde.T @ r + for _ in range(max_iter): + Ap = A @ p + alpha = rho / (p_tilde.T @ Ap) + x = x + alpha * p + r = r - alpha * Ap + r_tilde = r_tilde - alpha * (A.T @ p_tilde) + rho_new = r_tilde.T @ r + beta = rho_new / rho + p = r + beta * p + p_tilde = r_tilde + beta * p_tilde + if np.linalg.norm(r) < tol: + break + rho = rho_new + return x diff --git a/npbench/benchmarks/sparse/solvers/bicgstab/bicgstab.py b/npbench/benchmarks/sparse/solvers/bicgstab/bicgstab.py new file mode 100644 index 0000000..a764b1b --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/bicgstab/bicgstab.py @@ -0,0 +1,69 @@ +# Copyright 2023 University Politehnica of Bucharest and the NPBench authors. All rights reserved. + +import numpy as np + +# Generate required elements for a matrix with at least one non-zero term +# in the determinant +def required_elems(n): + rows = np.random.choice(n, n, replace=False) + cols = np.random.choice(n, n, replace=False) + vals = np.array(np.random.rand(n) * 10 - 5, dtype=np.float64) + return to_symetric(rows, cols, vals) + +# Generate more elements to make sure there are nnz or nnz + 1 +def ensure_nnz(n, nnz, rows, cols, vals): + nnz -= len(rows) + if nnz <= 0: + return + coords = set(zip(rows, cols)) + # Generate at least as many random values as they are required + new_vals = np.array(np.random.rand(nnz) * 10 - 5, dtype=np.float64) + i = 0 + while nnz > 0: + x, y = np.random.choice(n, 2) + stop = y + # The following loop changes x and y untill they are a new pair of coordinates + while (x, y) in coords: + y += 1 + if y == n: + y = 0 + if y == stop: + x += 1 + if x == n: + x = 0 + generated_pair = (x, y) + coords.add(generated_pair) + rows.append(generated_pair[0]) + cols.append(generated_pair[1]) + vals.append(new_vals[i]) + nnz -= 1 + if not generated_pair[0] == generated_pair[1]: + coords.add((generated_pair[1], generated_pair[0])) + rows.append(generated_pair[1]) + cols.append(generated_pair[0]) + vals.append(new_vals[i]) + nnz -= 1 + i += 1 + +def to_symetric(rows, cols, vals): + new_cols, new_rows, new_vals = [], [], [] + for i in range(0, vals.size): + new_rows.append(rows[i]) + new_cols.append(cols[i]) + new_vals.append(vals[i]) + if not rows[i] == cols[i]: + new_rows.append(cols[i]) + new_cols.append(rows[i]) + new_vals.append(vals[i]) + return new_rows, new_cols, new_vals + +def initialize(n : int, nnz : int, dtype=np.float64): + np.random.seed(42) + rows, cols, vals = required_elems(n) + ensure_nnz(n, nnz, rows, cols, vals) + from scipy.sparse import coo_matrix + A = coo_matrix((vals, (rows, cols)), shape=(n, n)).asformat('csr') + b = A @ np.random.rand(n) + # Generate starting solution + x = np.random.rand(n) + return A, x, b diff --git a/npbench/benchmarks/sparse/solvers/bicgstab/bicgstab_numpy.py b/npbench/benchmarks/sparse/solvers/bicgstab/bicgstab_numpy.py new file mode 100644 index 0000000..bd1c1c1 --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/bicgstab/bicgstab_numpy.py @@ -0,0 +1,25 @@ +import numpy as np + +# Solves A @ x = b where A is a Compressed Sparse Row matrix using the Biconjugate Gradient Stabilized method +def bicgstab(A, b, x, max_iter=100, tol=np.float64(1e-6)): + n = A.shape[0] + r = b - A @ x + rho_prev = alpha = omega = 1.0 + p = v = np.zeros_like(b) + r_tilde = np.copy(r) + x = np.copy(x) + for i in range(max_iter): + rho = np.dot(r_tilde.T, r) + beta = (rho / rho_prev) * (alpha / omega) + p = r + beta * (p - omega * v) + v = A @ p + alpha = rho / np.dot(r_tilde.T, v) + s = r - alpha * v + t = A @ s + omega = np.dot(t.T, s) / np.dot(t.T, t) + x = x + alpha * p + omega * s + r = s - omega * t + if np.linalg.norm(r) < tol: + break + rho_prev = rho + return x diff --git a/npbench/benchmarks/sparse/solvers/cg/cg.py b/npbench/benchmarks/sparse/solvers/cg/cg.py new file mode 100644 index 0000000..1a22a84 --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/cg/cg.py @@ -0,0 +1,71 @@ +# Copyright 2023 University Politehnica of Bucharest and the NPBench authors. All rights reserved. + +import numpy as np + +# Generate required elements for a matrix with at least one non-zero term +# in the determinant +def required_elems(n): + rows = np.random.choice(n, n, replace=False) + cols = np.random.choice(n, n, replace=False) + vals = np.array(np.random.rand(n) * 10 - 5, dtype=np.float64) + return to_symetric(rows, cols, vals) + +# Generate more elements to make sure there are nnz or nnz + 1 +def ensure_nnz(n, nnz, rows, cols, vals): + nnz -= len(rows) + if nnz <= 0: + return + coords = set(zip(rows, cols)) + # Generate at least as many random values as they are required + new_vals = np.array(np.random.rand(nnz) * 10 - 5, dtype=np.float64) + i = 0 + while nnz > 0: + x, y = np.random.choice(n, 2) + stop = y + # The following loop changes x and y untill they are a new pair of coordinates + while (x, y) in coords: + y += 1 + if y == n: + y = 0 + if y == stop: + x += 1 + if x == n: + x = 0 + generated_pair = (x, y) + coords.add(generated_pair) + rows.append(generated_pair[0]) + cols.append(generated_pair[1]) + vals.append(new_vals[i]) + nnz -= 1 + # Adds the symetric of the newly added element in the set + if not generated_pair[0] == generated_pair[1]: + coords.add((generated_pair[1], generated_pair[0])) + rows.append(generated_pair[1]) + cols.append(generated_pair[0]) + vals.append(new_vals[i]) + nnz -= 1 + i += 1 + +# Adds the symetrics of the existent elements in the matrix +def to_symetric(rows, cols, vals): + new_cols, new_rows, new_vals = [], [], [] + for i in range(0, vals.size): + new_rows.append(rows[i]) + new_cols.append(cols[i]) + new_vals.append(vals[i]) + if not rows[i] == cols[i]: + new_rows.append(cols[i]) + new_cols.append(rows[i]) + new_vals.append(vals[i]) + return new_rows, new_cols, new_vals + +def initialize(n : int, nnz : int, dtype=np.float64): + np.random.seed(42) + rows, cols, vals = required_elems(n) + ensure_nnz(n, nnz, rows, cols, vals) + from scipy.sparse import coo_matrix + A = coo_matrix((vals, (rows, cols)), shape=(n, n)).asformat('csr') + b = A @ np.random.rand(n) + # Generate starting solution + x = np.random.rand(n) + return A, x, b diff --git a/npbench/benchmarks/sparse/solvers/cg/cg_numpy.py b/npbench/benchmarks/sparse/solvers/cg/cg_numpy.py new file mode 100644 index 0000000..e4fddd9 --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/cg/cg_numpy.py @@ -0,0 +1,18 @@ +import numpy as np + +# Solves A @ x = b where A is a Compressed Sparse Row matrix using the Conjugate Gradient method +def cg(A, x, b, max_iter=100, tol=np.float64(1e-6)): + r = b - A @ x + p = r + rsold = r @ r + for i in range(max_iter): + Ap = A @ p + alpha = rsold / (p @ Ap) + x = x + alpha * p + r = r - alpha * Ap + rsnew = r @ r + if np.sqrt(rsnew) < tol: + break + p = r + (rsnew / rsold) * p + rsold = rsnew + return x diff --git a/npbench/benchmarks/sparse/solvers/gmres/gmres.py b/npbench/benchmarks/sparse/solvers/gmres/gmres.py new file mode 100644 index 0000000..a951ad9 --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/gmres/gmres.py @@ -0,0 +1,52 @@ +# Copyright 2023 University Politehnica of Bucharest and the NPBench authors. All rights reserved. + +import numpy as np + +# Generate required elements for a matrix with at least one non-zero term +# in the determinant +def required_elems(n): + rows = np.random.choice(n, n, replace=False) + cols = np.random.choice(n, n, replace=False) + vals = np.array(np.random.rand(n) * 10 - 5, dtype=np.float64) + return list(rows), list(cols), list(vals) + +# Generate more elements to make sure there are nnz or nnz + 1 +def ensure_nnz(n, nnz, rows, cols, vals): + nnz -= len(rows) + # No need to add more elements if there are already enough + if nnz <= 0: + return + coords = set(zip(rows, cols)) + # Generate at least as many random values as they are required + new_vals = np.array(np.random.rand(nnz) * 10 - 5, dtype=np.float64) + i = 0 + while nnz > 0: + x, y = np.random.choice(n, 2) + stop = y + # The following loop changes x and y untill they are a new pair of coordinates + while (x, y) in coords: + y += 1 + if y == n: + y = 0 + if y == stop: + x += 1 + if x == n: + x = 0 + generated_pair = (x, y) + coords.add(generated_pair) + rows.append(generated_pair[0]) + cols.append(generated_pair[1]) + vals.append(new_vals[i]) + nnz -= 1 + i += 1 + +def initialize(n : int, nnz : int, dtype=np.float64): + np.random.seed(42) + rows, cols, vals = required_elems(n) + ensure_nnz(n, nnz, rows, cols, vals) + from scipy.sparse import coo_matrix + A = coo_matrix((vals, (rows, cols)), shape=(n, n)).asformat('csr') + b = A @ np.random.rand(n) + # Generate starting solution + x = np.random.rand(n) + return A, x, b diff --git a/npbench/benchmarks/sparse/solvers/gmres/gmres_numpy.py b/npbench/benchmarks/sparse/solvers/gmres/gmres_numpy.py new file mode 100644 index 0000000..89553c0 --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/gmres/gmres_numpy.py @@ -0,0 +1,36 @@ +import numpy as np + +# Solves A @ x = b where A is a Compressed Sparse Row matrix using the Generalized Minimum Residual method +def hand_gmres(A, x, b, max_iter=100, tol=1e-6): + n = b.shape[0] + # Setting the dimensions of the Krylov subspace + m = min(max_iter, n) + + Q = np.empty((n, m + 1)) + H = np.zeros((m + 1, m)) + + r = b - A @ x + beta = np.linalg.norm(r) + Q[:, 0] = r / beta + + for k in range(m): + y = A @ Q[:, k] + for j in range(k + 1): + H[j, k] = Q[:, j] @ y + y -= H[j, k] * Q[:, j] + H[k + 1, k] = np.linalg.norm(y) + + if abs(H[k + 1, k]) < tol: + m = k + 1 + break + + Q[:, k + 1] = y / H[k + 1, k] + + e1 = np.zeros(m + 1) + e1[0] = 1.0 + + y = np.linalg.lstsq(H[:m, :], beta * e1[:m], rcond=None)[0] + + x += Q[:, :m] @ y + + return x diff --git a/npbench/benchmarks/sparse/solvers/minres/minres.py b/npbench/benchmarks/sparse/solvers/minres/minres.py new file mode 100644 index 0000000..e682f8b --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/minres/minres.py @@ -0,0 +1,72 @@ +# Copyright 2023 University Politehnica of Bucharest and the NPBench authors. All rights reserved. + +import numpy as np + +# Generate required elements for a matrix with at least one non-zero term +# in the determinant +def required_elems(n): + rows = np.random.choice(n, n, replace=False) + cols = np.random.choice(n, n, replace=False) + vals = np.array(np.random.rand(n) * 10 - 5, dtype=np.float64) + return to_symetric(rows, cols, vals) + +# Generate more elements to make sure there are nnz or nnz + 1 +def ensure_nnz(n, nnz, rows, cols, vals): + nnz -= len(rows) + # No need to add more elements if there are already enough + if nnz <= 0: + return + coords = set(zip(rows, cols)) + # Generate at least as many random values as they are required + new_vals = np.array(np.random.rand(nnz) * 10 - 5, dtype=np.float64) + i = 0 + while nnz > 0: + x, y = np.random.choice(n, 2) + stop = y + # The following loop changes x and y untill they are a new pair of coordinates + while (x, y) in coords: + y += 1 + if y == n: + y = 0 + if y == stop: + x += 1 + if x == n: + x = 0 + generated_pair = (x, y) + coords.add(generated_pair) + rows.append(generated_pair[0]) + cols.append(generated_pair[1]) + vals.append(new_vals[i]) + nnz -= 1 + # Adding the same element on the symetric position + if not generated_pair[0] == generated_pair[1]: + coords.add((generated_pair[1], generated_pair[0])) + rows.append(generated_pair[1]) + cols.append(generated_pair[0]) + vals.append(new_vals[i]) + nnz -= 1 + i += 1 + +# Adds the symetrics of the existent elements in the matrix +def to_symetric(rows, cols, vals): + new_cols, new_rows, new_vals = [], [], [] + for i in range(0, vals.size): + new_rows.append(rows[i]) + new_cols.append(cols[i]) + new_vals.append(vals[i]) + if not rows[i] == cols[i]: + new_rows.append(cols[i]) + new_cols.append(rows[i]) + new_vals.append(vals[i]) + return new_rows, new_cols, new_vals + +def initialize(n : int, nnz : int, dtype = np.float64): + np.random.seed(42) + rows, cols, vals = required_elems(n) + ensure_nnz(n, nnz, rows, cols, vals) + from scipy.sparse import coo_matrix + A = coo_matrix((vals, (rows, cols)), shape=(n, n)).asformat('csr') + b = A @ np.random.rand(n) + # Generate starting solution + x = np.random.rand(n) + return A, x, b diff --git a/npbench/benchmarks/sparse/solvers/minres/minres_numpy.py b/npbench/benchmarks/sparse/solvers/minres/minres_numpy.py new file mode 100644 index 0000000..60ac91c --- /dev/null +++ b/npbench/benchmarks/sparse/solvers/minres/minres_numpy.py @@ -0,0 +1,19 @@ +import numpy as np + +# Solves A @ x = b where A is a Compressed Sparse Row matrix using the Minimum Residual method +def hand_minres(A, b, x, max_iter=100, tol=1e-6): + # Residual vector + r = b - A @ x + # Initial search direction + p = r + for _ in range(max_iter): + Ap = A @ p + alpha = (r @ r) / (p @ Ap) + x = x + alpha * p + r_new = r - alpha * Ap + if np.linalg.norm(r_new) < tol: + break + beta = (r_new @ r_new) / (r @ r) + p = r_new + beta * p + r = r_new + return x diff --git a/npbench/benchmarks/sparse/spmm/spmm.py b/npbench/benchmarks/sparse/spmm/spmm.py new file mode 100644 index 0000000..5b63f40 --- /dev/null +++ b/npbench/benchmarks/sparse/spmm/spmm.py @@ -0,0 +1,14 @@ +# Copyright 2023 University Politehnica of Bucharest and the NPBench authors. All rights reserved. + +import numpy as np + +def initialize(NI, NJ, NK, nnz_A, nnz_B, datatype=np.float64): + import scipy.sparse as sp + rng = np.random.default_rng(42) + alpha = datatype(0.8) + beta = datatype(0.3) + np.random.seed(42) + C = np.fromfunction(lambda i, j: np.random.rand(i.shape[0], j.shape[0]), (NI, NJ), dtype=datatype) + A = sp.random(NI, NK, density=nnz_A / (NI * NK), format='csr', dtype=datatype, random_state=rng) + B = sp.random(NK, NJ, density=nnz_B / (NK * NJ), format='csr', dtype=datatype, random_state=rng) + return alpha, beta, C, A, B diff --git a/npbench/benchmarks/sparse/spmm/spmm_numpy.py b/npbench/benchmarks/sparse/spmm/spmm_numpy.py new file mode 100644 index 0000000..9cb1b00 --- /dev/null +++ b/npbench/benchmarks/sparse/spmm/spmm_numpy.py @@ -0,0 +1,5 @@ +import numpy as np + +# Stores the result of the product of 2 Compressed Sparse Row matrices A and B in C as a dense matrice +def spmm(alpha, beta, C, A, B): + C[:] = alpha * A @ B + beta * C