Added sparse benchmarks

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": []
+	}
+}

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": []
+	}
+}

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": []
+    }
+}

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": []
+	}
+}

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": []
+    }
+}

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": []
+    }
+}

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": []
+	}
+}

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

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

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

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