diff --git a/go.sum b/go.sum
new file mode 100644
index 0000000..d45948d
--- /dev/null
+++ b/go.sum
@@ -0,0 +1,2 @@
+gonum.org/v1/gonum v0.17.0 h1:VbpOemQlsSMrYmn7T2OUvQ4dqxQXU+ouZFQsZOx50z4=
+gonum.org/v1/gonum v0.17.0/go.mod h1:El3tOrEuMpv2UdMrbNlKEh9vd86bmQ6vqIcDwxEOc1E=
diff --git a/src/main.go b/src/main.go
new file mode 100644
index 0000000..4f1a6ce
--- /dev/null
+++ b/src/main.go
@@ -0,0 +1,90 @@
+package main
+
+import (
+	"bufio"
+	"fmt"
+	"io"
+	"os"
+	"strconv"
+	"strings"
+
+	"github.com/james-bowman/sparse"
+	"gonum.org/v1/gonum/mat"
+)
+
+func matrix_read_market(path string) (*sparse.CSR, error) {
+	f, err := os.Open(path)
+	if err != nil {
+		return nil, fmt.Errorf("open %q: %w", path, err)
+	}
+
+	defer f.Close()
+
+	sc := bufio.NewScanner(f)
+
+	const maxCapacity = 1024 * 1024 // Long lines
+	buf := make([]byte, 64*1024)
+	sc.Buffer(buf, maxCapacity)
+
+	if !sc.Scan() {
+		if err := sc.Err(); err != nil {
+			return nil, fmt.Errorf("read header: %w", err)
+		}
+		return nil, fmt.Errorf("empty file!")
+	}
+
+	header := strings.Fields(strings.TrimSpace(sc.Text()))
+	if len(header) != 5 || header[0] != "%%MatrixMarket" {
+		return nil, fmt.Errorf("invalid matrix martket header: %q", sc.Text())
+	}
+
+	object := strings.ToLower(header[1])
+	format := strings.ToLower(header[2])
+	field := strings.ToLower(header[3])
+	symmetry := strings.ToLower(header[4])
+
+	if object != "matrix" {
+		return nil, fmt.Errorf("unsupported object: %q", object)
+	}
+	if format != "coordinate" {
+		return nil, fmt.Errorf("unsupported format %q (need coordinate)", format)
+	}
+
+	switch field {
+	case "real", "integer", "pattern":
+	default:
+		return nil, fmt.Errorf("unsupported field %q", field)
+	}
+
+	switch symmetry {
+	case "general", "symmetric":
+	default:
+		return nil, fmt.Errorf("unsupported symmetry %q", symmetry)
+	}
+
+	// Skip comments and then read the size line
+	var size_line string
+	for sc.Scan() {
+		line := strings.TrimSpace(sc.Text())
+		if line == "" || strings.HasPrefix(line, "%") {
+			continue // skip comment
+		}
+		sizeLine
+}
+
+func main() {
+	A := mat.NewDense(2, 2, []float64{
+		1, 2,
+		3, 4,
+	})
+
+	B := mat.NewDense(2, 2, []float64{
+		5, 6,
+		7, 8,
+	})
+
+	var C mat.Dense
+	C.Mul(A, B)
+
+	fmt.Printf("%v\n", mat.Formatted(&C))
+}