diff --git a/src/main.go b/src/main.go
index 968f19d..1dc44ec 100644
--- a/src/main.go
+++ b/src/main.go
@@ -13,10 +13,7 @@ import (
 	"gonum.org/v1/gonum/mat"
 )
 
-var (
-	gNumMatmuls = 10
-	gNumSolves  = 10
-)
+var gNumTestIterations = 30
 
 type SparseMatrixTiming struct {
 	Label string `json:"label"`
@@ -28,9 +25,9 @@ type SparseMatrixTiming struct {
 	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"`
+	SpMVRuns    int   `json:"spmv_runs"`
+	SpMVTotalNs int64 `json:"spmv_total_ns"`
+	SpMVAvgNs   int64 `json:"spmv_avg_ns"`
 }
 
 type SparseBenchmarkCase struct {
@@ -131,7 +128,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 := gNumMatmuls
+	numberOfMults := gNumTestIterations
 
 	fmt.Printf("NNZ before matmuls: %d\n", mCSR.NNZ())
 
@@ -154,38 +151,41 @@ 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)
+func timeSparseMatVec(bcase *SparseBenchmarkCase) {
+	A := bcase.Matrix
+	rows, cols := A.Dims()
 
-	var chol sparse.Cholesky
+	x := mat.NewVecDense(cols, nil)
+	y := mat.NewVecDense(rows, nil)
+
+	for i := 0; i < cols; i++ {
+		x.SetVec(i, 1.0)
+	}
 
 	// warmup
-	err := chol.SolveVecTo(x, b)
-	if err != nil {
-		panic(err)
+	A.MulVecTo(y.RawVector().Data, false, x.RawVector().Data)
+	norm := y.Norm(2)
+
+	if norm < 1e-12 {
+		panic("Norm of resulting SpMV is zero. Something went wrong.\n")
 	}
 
-	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)
+	numberOfRuns := gNumTestIterations
+
+	timeBegin := time.Now()
+
+	for i := 0; i < numberOfRuns; i++ {
+		// fmt.Printf("spmv iteration %d\n", i)
+		A.MulVecTo(y.RawVector().Data, false, x.RawVector().Data)
 	}
 
-	bcase.Timing.SolveRuns = numberOfSolves
-	bcase.Timing.SolveTotalNs = timeElapsed.Nanoseconds()
-	bcase.Timing.SolveAvgNs = timeElapsed.Nanoseconds() / int64(numberOfSolves)
+	timeElapsed := time.Since(timeBegin)
+
+	total := timeElapsed.Nanoseconds()
+	avg := total / int64(numberOfRuns)
+	bcase.Timing.SpMVTotalNs = total
+	bcase.Timing.SpMVAvgNs = avg
+	bcase.Timing.SpMVRuns = numberOfRuns
 }
 
 func timeNanoToMS(timeNS int64) float64 {
@@ -201,16 +201,22 @@ func doTimings(path string) {
 		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)
+	doMatMul := true
+	if doMatMul {
+		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)
+	doSpMV := true
+	if doSpMV {
+		fmt.Printf("SpMV:\n")
+		timeSparseMatVec(&bcase)
+		avgSpMVTimeMS := timeNanoToMS(bcase.Timing.SpMVAvgNs)
+		fmt.Printf("Avg SpMV time for %s: %.4f ms\n", bcase.Timing.Label, avgSpMVTimeMS)
+	}
 }
 
 func main() {