diff --git a/src/main.go b/src/main.go
index 46f0d20..968f19d 100644
--- a/src/main.go
+++ b/src/main.go
@@ -10,6 +10,12 @@ import (
 	"time"
 
 	"github.com/james-bowman/sparse"
+	"gonum.org/v1/gonum/mat"
+)
+
+var (
+	gNumMatmuls = 10
+	gNumSolves  = 10
 )
 
 type SparseMatrixTiming struct {
@@ -116,15 +122,6 @@ func getSparseBenchmarkCase(path string) SparseBenchmarkCase {
 }
 
 func timeSparseMatmuls(bcase *SparseBenchmarkCase) {
-	rows := bcase.Timing.Rows
-	cols := bcase.Timing.Cols
-
-	if rows != cols {
-		panic("sparse matmul benchmark assumes square matrix")
-	}
-
-	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 {
@@ -134,7 +131,7 @@ func timeSparseMatmuls(bcase *SparseBenchmarkCase) {
 	// The CSR x CSR matrix multiplication is supposed to have a specific optimised
 	// path for matmuls
 	var out sparse.CSR
-	numberOfMults := 100
+	numberOfMults := gNumMatmuls
 
 	fmt.Printf("NNZ before matmuls: %d\n", mCSR.NNZ())
 
@@ -157,15 +154,63 @@ func timeSparseMatmuls(bcase *SparseBenchmarkCase) {
 	bcase.Timing.MatMulAvgNs = timeAvgNS
 }
 
+func timeSparseSolve(bcase *SparseBenchmarkCase) {
+	matrix := bcase.Matrix
+	rows, _ := matrix.Dims()
+	rhsData := make([]float64, rows) // We have checked square matrix already
+	for i := range rhsData {
+		rhsData[i] = 1.0
+	}
+	b := mat.NewVecDense(rows, rhsData)
+	x := mat.NewVecDense(rows, nil)
+
+	var chol sparse.Cholesky
+
+	// warmup
+	err := chol.SolveVecTo(x, b)
+	if err != nil {
+		panic(err)
+	}
+
+	numberOfSolves := gNumSolves
+	for i := 0; i < 3; i += 1 {
+		timeBegin := time.Now()
+		chol.Factorize(matrix)
+		err = chol.SolveVecTo(x, b)
+		if err != nil {
+			panic(err)
+		}
+		timeElapsed := time.Since(timeBegin)
+	}
+
+	bcase.Timing.SolveRuns = numberOfSolves
+	bcase.Timing.SolveTotalNs = timeElapsed.Nanoseconds()
+	bcase.Timing.SolveAvgNs = timeElapsed.Nanoseconds() / int64(numberOfSolves)
+}
+
 func timeNanoToMS(timeNS int64) float64 {
 	return float64(timeNS) / float64(1e6)
 }
 
 func doTimings(path string) {
 	bcase := getSparseBenchmarkCase(path)
+	rows := bcase.Timing.Rows
+	cols := bcase.Timing.Cols
+
+	if rows != cols {
+		panic("Sparse benchmark assumes square matrix")
+	}
+
+	fmt.Printf("Timing sparse matrix %s with size: %d x %d \n", bcase.Timing.Label, rows, cols)
+	fmt.Printf("Matmul:\n")
 	timeSparseMatmuls(&bcase)
 	avgMatMulTimeMS := timeNanoToMS(bcase.Timing.MatMulAvgNs)
 	fmt.Printf("Avg matmul time for %s: %.4f ms \n", bcase.Timing.Label, avgMatMulTimeMS)
+
+	fmt.Printf("Cholesky solve:\n")
+	timeSparseSolve(&bcase)
+	avgSolveTimeMS := timeNanoToMS(bcase.Timing.SolveAvgNs)
+	fmt.Printf("Avg Cholesky solve time for %s: %.4f ms \n", bcase.Timing.Label, avgSolveTimeMS)
 }
 
 func main() {