From f81c97bb1e0f2b167c650356360f37c5a1cf5201 Mon Sep 17 00:00:00 2001
From: antonl <anton.ljungdahl@gmail.com>
Date: Sat, 14 Mar 2026 21:45:42 +0100
Subject: [PATCH] some more work

---
 src/main.go | 126 +++++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 100 insertions(+), 26 deletions(-)

diff --git a/src/main.go b/src/main.go
index 370f87b..46f0d20 100644
--- a/src/main.go
+++ b/src/main.go
@@ -4,6 +4,7 @@ import (
 	"bufio"
 	"fmt"
 	"os"
+	"path/filepath"
 	"strconv"
 	"strings"
 	"time"
@@ -11,9 +12,29 @@ import (
 	"github.com/james-bowman/sparse"
 )
 
+type SparseMatrixTiming struct {
+	Label string `json:"label"`
+	Rows  int    `json:"rows"`
+	Cols  int    `json:"cols"`
+	NNZ   int    `json:"nnz"`
+
+	MatMulRuns    int   `json:"matmul_runs"`
+	MatMulTotalNs int64 `json:"matmul_total_ns"`
+	MatMulAvgNs   int64 `json:"matmul_avg_ns"`
+
+	SolveRuns    int   `json:"solve_runs"`
+	SolveTotalNs int64 `json:"solve_total_ns"`
+	SolveAvgNs   int64 `json:"solve_avg_ns"`
+}
+
+type SparseBenchmarkCase struct {
+	Timing SparseMatrixTiming
+	Matrix *sparse.CSR
+}
+
 // The "SuiteSparse" collection of matrices come in a format called
 // "market matrix" and so we parse and load them to a sparse.COO format here.
-func matrix_load_market(path string) *sparse.COO {
+func matrixLoadMarket(path string) *sparse.CSR {
 	f, err := os.Open(path)
 	if err != nil {
 		panic(err)
@@ -60,47 +81,100 @@ func matrix_load_market(path string) *sparse.COO {
 
 	coo := sparse.NewCOO(rows, cols, r, c, v)
 
-	return coo
+	return coo.ToCSR()
 }
 
-func main() {
-	// small_test_matrix_path := "suitesparse_test_matrices/1138_bus.mtx"
-	big_8mill_nnz_matrix_path := "suitesparse_test_matrices/thermal2.mtx"
-	m := matrix_load_market(big_8mill_nnz_matrix_path)
-	rows, cols := m.Dims()
-	fmt.Printf("matrix size: %d x %d \n", rows, cols)
+// small helper just to strip path and file ending
+func matrixGetLabel(path string) string {
+	base := filepath.Base(path)          // f. ex "1138_bus.mtx"
+	ext := filepath.Ext(base)            // ".mtx"
+	out := strings.TrimSuffix(base, ext) // just "1138_bus"
+	return out
+}
+
+func getSparseBenchmarkCase(path string) SparseBenchmarkCase {
+	var out SparseBenchmarkCase
+
+	loadMarketStart := time.Now()
+	label := matrixGetLabel(path)
+	m := matrixLoadMarket(path)
+	loadMarketEnd := time.Since(loadMarketStart)
+	loadedMarketMillisecs := float64(loadMarketEnd.Milliseconds())
+	fmt.Printf("Loaded market matrix %s in %.4f ms \n", label, loadedMarketMillisecs)
+
+	out.Matrix = m
+	out.Timing.Label = label
+
+	rows, cols := out.Matrix.Dims()
+	nnz := m.NNZ()
+
+	out.Timing.Rows = rows
+	out.Timing.Cols = cols
+	out.Timing.NNZ = nnz
+
+	return out
+}
+
+func timeSparseMatmuls(bcase *SparseBenchmarkCase) {
+	rows := bcase.Timing.Rows
+	cols := bcase.Timing.Cols
 
 	if rows != cols {
 		panic("sparse matmul benchmark assumes square matrix")
 	}
 
-	// Convert to Compressed-Sparse-Row-storage
-	m_csr := m.ToCSR()
-	// Warm up caches to not include "first call" effects
+	fmt.Printf("Timing sparse matrix %s with size: %d x %d \n", bcase.Timing.Label, rows, cols)
+
+	mCSR := bcase.Matrix
 	var warm sparse.CSR
 	for i := 0; i < 3; i += 1 {
-		warm.Mul(m_csr, m_csr)
+		warm.Mul(mCSR, mCSR)
 	}
 
 	// The CSR x CSR matrix multiplication is supposed to have a specific optimised
 	// path for matmuls
 	var out sparse.CSR
-	number_of_mults := 100
+	numberOfMults := 100
 
-	time_begin := time.Now()
-	for i := 0; i < number_of_mults; i += 1 {
-		out.Mul(m_csr, m_csr)
+	fmt.Printf("NNZ before matmuls: %d\n", mCSR.NNZ())
+
+	timeBegin := time.Now()
+	for i := 0; i < numberOfMults; i += 1 {
+		out.Mul(mCSR, mCSR)
+	}
+	timeElapsed := time.Since(timeBegin)
+
+	fmt.Printf("NNZ after matmuls: %d\n", out.NNZ())
+
+	if warm.NNZ() != out.NNZ() {
+		panic("Sparsity pattern changed unexpectedly after matmul!")
 	}
 
-	elapsed := time.Since(time_begin)
-	elapsed_ms := elapsed.Milliseconds()
+	bcase.Timing.MatMulRuns = numberOfMults
+	bcase.Timing.MatMulTotalNs = timeElapsed.Nanoseconds()
 
-	time_per_matmul := float32(elapsed_ms) / float32(number_of_mults)
-
-	out_rows, out_cols := out.Dims()
-	if out_rows != rows || out_cols != cols {
-		panic("Output matrix dimensions not same as input. We are supposed to matmul square matrices")
-	}
-
-	fmt.Printf("Elapsed time per mat mul: %.4f ms\n", time_per_matmul)
+	timeAvgNS := timeElapsed.Nanoseconds() / int64(numberOfMults)
+	bcase.Timing.MatMulAvgNs = timeAvgNS
+}
+
+func timeNanoToMS(timeNS int64) float64 {
+	return float64(timeNS) / float64(1e6)
+}
+
+func doTimings(path string) {
+	bcase := getSparseBenchmarkCase(path)
+	timeSparseMatmuls(&bcase)
+	avgMatMulTimeMS := timeNanoToMS(bcase.Timing.MatMulAvgNs)
+	fmt.Printf("Avg matmul time for %s: %.4f ms \n", bcase.Timing.Label, avgMatMulTimeMS)
+}
+
+func main() {
+	mainBegin := time.Now()
+
+	thermal2Path := "suitesparse_test_matrices/thermal2.mtx"
+	doTimings(thermal2Path)
+
+	mainElapsed := time.Since(mainBegin)
+
+	fmt.Printf("Program finished in %.4f seconds \n", float64(mainElapsed.Milliseconds())/1e3)
 }