Repository: alugowski/fast_matrix_market (language: python)

Issue:
Add GraphBLAS bindings



Discussion / hints:

Fix patch:
diff --git a/README.GraphBLAS.md b/README.GraphBLAS.md
new file mode 100644
index 0000000..bfb894a
--- /dev/null
+++ b/README.GraphBLAS.md
@@ -0,0 +1,62 @@
+# GraphBLAS Matrix Market
+
+`fast_matrix_market` provides read and write methods for GraphBLAS GrB_Matrix and GrB_Vector objects.
+
+Tested with [SuiteSparse:GraphBLAS](https://github.com/DrTimothyAldenDavis/GraphBLAS), but should work with any GraphBLAS implementation.
+
+# Usage
+
+```c++
+#include <fast_matrix_market/app/GraphBLAS.hpp>
+```
+
+```c++
+GrB_Matrix A;
+// or
+GrB_Vector A;
+```
+### Reading
+```c++
+std::ifstream f("input.mtx");
+
+fast_matrix_market::read_matrix_market_graphblas(f, &A);
+```
+
+**A's GrB_Type** is determined by the first non-null condition:
+* Explicit `desired_type` argument to `read*()`. `A` will be read as this type.
+* `%%GraphBLAS type <ctype>` structured comment.
+* MatrixMarket header `field`.
+
+User types are not supported at this time. If you have a use case for user types stored in Matrix Market please get in touch.
+
+**Note:** if `A` is a `GrB_Vector` then the Matrix Market file must be either a `vector` or a
+1-by-N row or M-by-1 column matrix. Any other matrix will cause an exception to be thrown.
+
+### Writing
+
+```c++
+std::ofstream f("output.mtx");
+
+fast_matrix_market::write_matrix_market_graphblas(f, A);
+```
+
+**Compatibility Note:** `GrB_Vector` objects are written as `vector` Matrix Market files.
+Many MatrixMarket loaders cannot read vector files. If this is a problem write the vector
+as a 1-by-N or M-by-1 matrix.
+
+# GxB Extensions Used
+
+The following GraphBLAS extensions are used. All are individually switchable. All are present in SuiteSparse:GraphBLAS.
+* `GxB_FC32`/`GxB_FC64` complex types. Standard GraphBLAS does not include complex number types.
+  * Define `FMM_GXB_COMPLEX=0` to disable.
+* `GxB_Matrix_build_Scalar` is used to build iso-valued matrices from `pattern` files.
+  * Define `FMM_GXB_BUILD_SCALAR=0` to disable.
+* `GxB_Matrix_pack_FullC` is used to efficiently read `array` files.
+  * Define `FMM_GXB_PACK_UNPACK=0` to disable
+* `GxB_Iterator` for zero-copy writes.
+  * Define `FMM_GXB_ITERATORS=0` to disable.
+* `GxB_Matrix_type_name` to determine matrix type.
+  * Define `FMM_GXB_TYPE_NAME=0` to disable.
+
+All extensions can be disabled by defining `FMM_NO_GXB`.
+
diff --git a/README.md b/README.md
index d2c8a80..0dd5d73 100644
--- a/README.md
+++ b/README.md
@@ -4,6 +4,9 @@
 
 A fast and full-featured Matrix Market I/O library for C++ and [Python](python).
 
+Ready-to-use bindings for [GraphBLAS](README.GraphBLAS.md), [Eigen](README.Eigen.md), [CXSparse](README.CXSparse.md), [Blaze](README.Blaze.md), [SciPy](python/README.md), `std::vector`.  
+Easy to integrate with any datastructure.
+
 [Matrix Market](https://math.nist.gov/MatrixMarket/formats.html) is a simple, human-readable, and widely used sparse matrix file format that looks like this:
 ```
 %%MatrixMarket matrix coordinate real general
@@ -26,8 +29,6 @@ hit or miss as some valid matrices cannot be loaded.
 
 Use this library to fix these shortcomings.
 
-Ready-to-use bindings for [Eigen](README.Eigen.md), [CXSparse](README.CXSparse.md), [Blaze](README.Blaze.md), [SciPy](python/README.md).
-
 # Fast
 
 `fast_matrix_market` takes the fastest available sequential methods and parallelizes them.
@@ -112,6 +113,14 @@ fast_matrix_market::read_matrix_market_array(
 
 1D dense vectors supported by the same method.
 
+## GraphBLAS
+`GrB_Matrix` and `GrB_Vector`s are supported, with zero-copy where possible. See [GraphBLAS README](README.GraphBLAS.md).
+```c++
+GrB_Matrix A;
+fast_matrix_market::read_matrix_market_graphblas(input_stream, &A);
+```
+
+
 ## Eigen
 Sparse and dense matrices and vectors are supported. See [Eigen README](README.Eigen.md).
 ```c++
diff --git a/cmake/GraphBLAS.cmake b/cmake/GraphBLAS.cmake
new file mode 100644
index 0000000..6eb8589
--- /dev/null
+++ b/cmake/GraphBLAS.cmake
@@ -0,0 +1,43 @@
+# Attempt at a general-purpose GraphBLAS locator.
+#
+# Known to work on:
+# * macOS with `brew install suite-sparse` (laptop and GitHub Actions macos)
+# * Ubuntu with `sudo apt-get install -y libsuitesparse-dev` (GitHub Actions)
+#
+
+# Attempt to use FindGraphBLAS
+cmake_policy(SET CMP0074 NEW)
+find_package(GraphBLAS QUIET MODULE)
+
+# See if GraphBLAS is installed system-wide already
+if (NOT GraphBLAS_FOUND)
+    check_source_compiles(CXX "
+    #include <GraphBLAS.h>
+    int main(void) { return 0; }
+    " GRAPHBLAS_H_AVAILABLE)
+
+    if (GRAPHBLAS_H_AVAILABLE)
+        set(GraphBLAS_FOUND TRUE)
+        set(GRAPHBLAS_INCLUDE_DIR "")
+        set(GRAPHBLAS_LIBRARIES "graphblas")
+    #    set(GRAPHBLAS_STATIC "")
+    #    set(GRAPHBLAS_LIBRARY "")
+    endif()
+endif()
+
+# See if brew is installed, and if so check for suite-sparse
+if (NOT GraphBLAS_FOUND AND ${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+    find_program(BREW "brew")
+    if (BREW)
+        execute_process(COMMAND brew --prefix suite-sparse OUTPUT_VARIABLE PREFIX OUTPUT_STRIP_TRAILING_WHITESPACE ERROR_QUIET)
+        if (NOT ("${PREFIX}" STREQUAL "" ))
+            set(GraphBLAS_ROOT "${PREFIX}/lib/cmake/SuiteSparse")
+            message("suite-sparse found: ${PREFIX}")
+            message("checking: ${GraphBLAS_ROOT}")
+            set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH};${GraphBLAS_ROOT}")
+            find_package(GraphBLAS QUIET MODULE)
+        else()
+            message("Install GraphBLAS with 'brew install suite-sparse'")
+        endif()
+    endif()
+endif()
\ No newline at end of file
diff --git a/codecov.yml b/codecov.yml
index 511913b..e0fef4b 100644
--- a/codecov.yml
+++ b/codecov.yml
@@ -1,5 +1,12 @@
 coverage:
   range: 70..90
+  status:
+    project:
+      default:
+        informational: true
+    patch:
+      default:
+        informational: true
 parsers:
   gcov:
     branch_detection:
diff --git a/include/fast_matrix_market/app/GraphBLAS.hpp b/include/fast_matrix_market/app/GraphBLAS.hpp
new file mode 100644
index 0000000..e9768da
--- /dev/null
+++ b/include/fast_matrix_market/app/GraphBLAS.hpp
@@ -0,0 +1,1427 @@
+// Copyright (C) 2022-2023 Adam Lugowski. All rights reserved.
+// Use of this source code is governed by the BSD 2-clause license found in the LICENSE.txt file.
+
+#pragma once
+
+#include "../fast_matrix_market.hpp"
+
+extern "C" {
+#include <GraphBLAS.h>
+}
+
+// Define this to disable all extensions.
+//#define FMM_NO_GXB
+
+#ifdef FMM_NO_GXB
+// Disable all extensions
+#define FMM_GXB_COMPLEX 0
+#define FMM_GXB_BUILD_SCALAR 0
+#define FMM_GXB_ITERATORS 0
+#define FMM_GXB_PACK_UNPACK 0
+#define FMM_GXB_TYPE_NAME 0
+#endif
+
+// Switch to allow using the SuiteSparse:GraphBLAS complex number extension
+#ifndef FMM_GXB_COMPLEX
+#define FMM_GXB_COMPLEX 1
+#endif
+
+// Switch to allow using GxB_Matrix_build_Scalar
+#ifndef FMM_GXB_BUILD_SCALAR
+#define FMM_GXB_BUILD_SCALAR 1
+#endif
+
+// Switch to allow using the SuiteSparse:GraphBLAS iterator extension
+#ifndef FMM_GXB_ITERATORS
+#define FMM_GXB_ITERATORS 1
+#endif
+
+// Switch to allow using the SuiteSparse:GraphBLAS iterators for vectors.
+// The author is not fully confident on the correctness of this implementation, even though it does pass
+// tests, because the docs and behavior leave some ambiguities. For example, how can pmax (the largest index
+// of an element) be larger than nvals (the number of elements)?
+// Enable this if you've checked over the implementation and need the zero copy.
+// For everyone else, just use the extractTuples version.
+#ifndef FMM_GXB_VECTOR_ITERATORS
+#define FMM_GXB_VECTOR_ITERATORS 0
+#endif
+
+// Switch to allow using pack/unpack methods
+#ifndef FMM_GXB_PACK_UNPACK
+#define FMM_GXB_PACK_UNPACK 1
+#endif
+
+// Switch to allow using pack/unpack methods
+#ifndef FMM_GXB_TYPE_NAME
+#define FMM_GXB_TYPE_NAME 1
+#endif
+
+
+namespace fast_matrix_market {
+
+    /**
+     * Check that a GrB function call was successful and throw an exception if it wasn't.
+     */
+    void ok(GrB_Info ec) {
+        if (ec != GrB_SUCCESS) {
+            throw fmm_error("GraphBLAS method returned error code: " + std::to_string(ec));
+        }
+    }
+
+    /**
+     * Return map[key], throwing an exception with a helpful error message if key is not in map.
+     */
+    template <typename ENUM>
+    ENUM parse_key(const std::string& s, std::map<const std::string, ENUM> mp, const std::string& what) {
+        auto iter = mp.find(s);
+        if (iter != mp.end()) {
+            return iter->second;
+        }
+
+        std::string acceptable;
+        std::string delim;
+        for (const auto& [key, value] : mp) {
+            acceptable += delim + std::string(key);
+            delim = ", ";
+        }
+        throw invalid_argument(std::string("Invalid ") + what + ". Must be one of: " + acceptable);
+    }
+
+    /**
+     * Map ctype strings to GraphBLAS GrB_Type.
+     */
+    inline std::map<const std::string, GrB_Type> str_ctype_to_GrB_Type() {
+        return std::map<const std::string, GrB_Type>{
+                {"bool", GrB_BOOL},
+                {"int8_t", GrB_INT8},
+                {"int16_t", GrB_INT16},
+                {"int32_t", GrB_INT32},
+                {"int64_t", GrB_INT64},
+                {"uint8_t", GrB_UINT8},
+                {"uint16_t", GrB_UINT16},
+                {"uint32_t", GrB_UINT32},
+                {"uint64_t", GrB_UINT64},
+                {"float", GrB_FP32},
+                {"double", GrB_FP64},
+#if FMM_GXB_COMPLEX
+                {"float complex", GxB_FC32},
+                {"double complex", GxB_FC64},
+                {"GxB_FC32_t", GxB_FC32},
+                {"GxB_FC64_t", GxB_FC64},
+#endif
+        };
+    }
+
+    /**
+     * Map GrB_Type to a ctype string. This is the reverse of `str_ctype_to_GrB_Type`.
+     */
+    inline std::map<GrB_Type, std::string> GrB_Type_to_header_type() {
+        std::map<GrB_Type, std::string> ret;
+        for (const auto& [str, type] : str_ctype_to_GrB_Type()) {
+            ret[type] = str;
+        }
+
+#if FMM_GXB_COMPLEX
+        // Explicitly define complex values
+        ret[GxB_FC32] = "float complex";
+        ret[GxB_FC64] = "double complex";
+#endif
+        return ret;
+    }
+
+    /**
+     * Parse a type string into a GrB_Type, with error checking.
+     */
+    GrB_Type parse_type(const std::string& s) {
+        return parse_key(s, str_ctype_to_GrB_Type(), "type");
+    }
+
+    /**
+     * GraphBLAS typed functions are duplicated with their type part of the function name. This makes any call
+     * to these functions require a large switch statement. This template specializes each explicit GraphBLAS type
+     * so the logic methods can just use regular C++ templates.
+     *
+     * @tparam T C++ type that has a mapping to GraphBLAS.
+     */
+    template <typename T>
+    struct GraphBLAS_typed {
+        static GrB_Type type() { return nullptr; }
+
+        static GrB_Info build_matrix([[maybe_unused]] GrB_Matrix mat, [[maybe_unused]] const GrB_Index* rows, [[maybe_unused]] const GrB_Index* cols, [[maybe_unused]] const T* vals, [[maybe_unused]] GrB_Index nvals) {
+            return GrB_PANIC;
+        }
+
+        static GrB_Info set_element([[maybe_unused]] GrB_Scalar scalar, [[maybe_unused]] const T& x) {
+            return GrB_PANIC;
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples([[maybe_unused]] GrB_Index *I, [[maybe_unused]] T *X, [[maybe_unused]] GrB_Index *nvals, [[maybe_unused]] const GrB_Matrix& A) {
+            return GrB_PANIC;
+        }
+
+#if FMM_GXB_ITERATORS
+        static T GxB_Iterator_get([[maybe_unused]] GxB_Iterator iterator) {
+            throw fmm_error("FMM Bug: GxB_Iterator_get<T> called.");
+        }
+#endif
+        static GrB_Info build_vector([[maybe_unused]] GrB_Vector vec, [[maybe_unused]] const GrB_Index* indices, [[maybe_unused]] const T* vals, [[maybe_unused]] GrB_Index nvals) {
+            return GrB_PANIC;
+        }
+
+        static GrB_Info GrB_Vector_extractTuples([[maybe_unused]] GrB_Index *I, [[maybe_unused]] T *X, [[maybe_unused]] GrB_Index *nvals, [[maybe_unused]] const GrB_Vector& A) {
+            return GrB_PANIC;
+        }
+    };
+    template <> struct GraphBLAS_typed<bool> {
+        static GrB_Type type() { return GrB_BOOL; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const bool* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_BOOL(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, bool x) {
+            return GrB_Scalar_setElement_BOOL(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, bool *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_BOOL(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static bool GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_BOOL(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const bool* vals, GrB_Index nvals) {
+            return GrB_Vector_build_BOOL(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, bool *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_BOOL(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<int8_t> {
+        static GrB_Type type() { return GrB_INT8; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const int8_t* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_INT8(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, int8_t x) {
+            return GrB_Scalar_setElement_INT8(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, int8_t *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_INT8(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static int8_t GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_INT8(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const int8_t* vals, GrB_Index nvals) {
+            return GrB_Vector_build_INT8(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, int8_t *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_INT8(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<int16_t> {
+        static GrB_Type type() { return GrB_INT16; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const int16_t* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_INT16(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, int16_t x) {
+            return GrB_Scalar_setElement_INT16(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, int16_t *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_INT16(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static int16_t GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_INT16(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const int16_t* vals, GrB_Index nvals) {
+            return GrB_Vector_build_INT16(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, int16_t *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_INT16(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<int32_t> {
+        static GrB_Type type() { return GrB_INT32; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const int32_t* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_INT32(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, int32_t x) {
+            return GrB_Scalar_setElement_INT32(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, int32_t *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_INT32(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static int32_t GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_INT32(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const int32_t* vals, GrB_Index nvals) {
+            return GrB_Vector_build_INT32(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, int32_t *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_INT32(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<int64_t> {
+        static GrB_Type type() { return GrB_INT64; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const int64_t* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_INT64(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, int64_t x) {
+            return GrB_Scalar_setElement_INT64(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, int64_t *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_INT64(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static int64_t GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_INT64(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const int64_t* vals, GrB_Index nvals) {
+            return GrB_Vector_build_INT64(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, int64_t *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_INT64(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<uint8_t> {
+        static GrB_Type type() { return GrB_UINT8; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const uint8_t* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_UINT8(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, uint8_t x) {
+            return GrB_Scalar_setElement_UINT8(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, uint8_t *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_UINT8(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static uint8_t GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_UINT8(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const uint8_t* vals, GrB_Index nvals) {
+            return GrB_Vector_build_UINT8(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, uint8_t *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_UINT8(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<uint16_t> {
+        static GrB_Type type() { return GrB_UINT16; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const uint16_t* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_UINT16(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, uint16_t x) {
+            return GrB_Scalar_setElement_UINT16(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, uint16_t *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_UINT16(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static uint16_t GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_UINT16(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const uint16_t* vals, GrB_Index nvals) {
+            return GrB_Vector_build_UINT16(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, uint16_t *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_UINT16(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<uint32_t> {
+        static GrB_Type type() { return GrB_UINT32; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const uint32_t* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_UINT32(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, uint32_t x) {
+            return GrB_Scalar_setElement_UINT32(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, uint32_t *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_UINT32(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static uint32_t GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_UINT32(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const uint32_t* vals, GrB_Index nvals) {
+            return GrB_Vector_build_UINT32(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, uint32_t *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_UINT32(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<uint64_t> {
+        static GrB_Type type() { return GrB_UINT64; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const uint64_t* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_UINT64(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, uint64_t x) {
+            return GrB_Scalar_setElement_UINT64(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, uint64_t *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_UINT64(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static uint64_t GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_UINT64(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const uint64_t* vals, GrB_Index nvals) {
+            return GrB_Vector_build_UINT64(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, uint64_t *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_UINT64(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<float> {
+        static GrB_Type type() { return GrB_FP32; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const float* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_FP32(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, float x) {
+            return GrB_Scalar_setElement_FP32(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, float *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_FP32(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static float GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_FP32(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const float* vals, GrB_Index nvals) {
+            return GrB_Vector_build_FP32(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, float *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_FP32(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<double> {
+        static GrB_Type type() { return GrB_FP64; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const double* vals, GrB_Index nvals) {
+            return GrB_Matrix_build_FP64(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, double x) {
+            return GrB_Scalar_setElement_FP64(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, double *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GrB_Matrix_extractTuples_FP64(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static double GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_FP64(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const double* vals, GrB_Index nvals) {
+            return GrB_Vector_build_FP64(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, double *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GrB_Vector_extractTuples_FP64(I, X, nvals, A);
+        }
+    };
+#if FMM_GXB_COMPLEX
+    template <> struct GraphBLAS_typed<std::complex<float>> {
+        static GrB_Type type() { return GxB_FC32; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const std::complex<float>* vals, GrB_Index nvals) {
+            return GxB_Matrix_build_FC32(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, const std::complex<float>& x) {
+            return GxB_Scalar_setElement_FC32(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, std::complex<float> *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GxB_Matrix_extractTuples_FC32(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static std::complex<float> GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_FC32(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const std::complex<float>* vals, GrB_Index nvals) {
+            return GxB_Vector_build_FC32(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, std::complex<float> *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GxB_Vector_extractTuples_FC32(I, X, nvals, A);
+        }
+    };
+    template <> struct GraphBLAS_typed<std::complex<double>> {
+        static GrB_Type type() { return GxB_FC64; }
+
+        static GrB_Info build_matrix(GrB_Matrix mat, const GrB_Index* rows, const GrB_Index* cols, const std::complex<double>* vals, GrB_Index nvals) {
+            return GxB_Matrix_build_FC64(mat, rows, cols, vals, nvals, nullptr);
+        }
+
+        static GrB_Info set_element(GrB_Scalar scalar, const std::complex<double>& x) {
+            return GxB_Scalar_setElement_FC64(scalar, x);
+        }
+
+        static GrB_Info GrB_Matrix_extractTuples(GrB_Index *I, GrB_Index *J, std::complex<double> *X, GrB_Index *nvals, const GrB_Matrix& A) {
+            return GxB_Matrix_extractTuples_FC64(I, J, X, nvals, A);
+        }
+
+#if FMM_GXB_ITERATORS
+        static std::complex<double> GxB_Iterator_get(GxB_Iterator iterator) {
+            return GxB_Iterator_get_FC64(iterator);
+        }
+#endif
+        static GrB_Info build_vector(GrB_Vector vec, const GrB_Index* indices, const std::complex<double>* vals, GrB_Index nvals) {
+            return GxB_Vector_build_FC64(vec, indices, vals, nvals, nullptr);
+        }
+
+        static GrB_Info GrB_Vector_extractTuples(GrB_Index *I, std::complex<double> *X, GrB_Index *nvals, const GrB_Vector& A) {
+            return GxB_Vector_extractTuples_FC64(I, X, nvals, A);
+        }
+    };
+#endif
+
+    /**
+     * Read a Matrix Market coordinate body into a sparse GraphBLAS matrix using triplets.
+     */
+    template <typename T>
+    void read_body_graphblas_coordinate(std::istream &instream,
+                                        matrix_market_header &header,
+                                        GrB_Matrix mat,
+                                        const read_options& options) {
+        size_t storage_nnz = get_storage_nnz(header, options);
+
+        // Read into triplets
+        std::vector<GrB_Index> rows(storage_nnz);
+        std::vector<GrB_Index> cols(storage_nnz);
+
+#if FMM_GXB_BUILD_SCALAR
+        if (header.field == pattern) {
+            // read the indices
+            auto handler = triplet_pattern_parse_handler(rows.begin(), cols.begin());
+            read_matrix_market_body_no_pattern(instream, header, handler, options);
+
+            // create the scalar
+            GrB_Scalar one;
+            ok(GrB_Scalar_new(&one, GraphBLAS_typed<T>::type()));
+            ok(GraphBLAS_typed<T>::set_element(one, pattern_default_value(static_cast<T*>(nullptr))));
+
+            // Build iso matrix
+            ok(GxB_Matrix_build_Scalar(mat, rows.data(), cols.data(), one, storage_nnz));
+
+            // clean up
+            ok(GrB_Scalar_free(&one));
+        } else
+#endif
+        {
+            // Allocate values. Cannot use std::vector due to bool specialization
+            auto vals = std::make_unique<T[]>(storage_nnz);
+
+            // read indices and values
+            auto handler = triplet_parse_handler(rows.begin(), cols.begin(), vals.get());
+            read_matrix_market_body(instream, header, handler, pattern_default_value(static_cast<T*>(nullptr)), options);
+
+            // build matrix from triplets
+            ok(GraphBLAS_typed<T>::build_matrix(mat, rows.data(), cols.data(), vals.get(), storage_nnz));
+        }
+    }
+
+
+    /**
+     * Read a Matrix Market array body into a full GraphBLAS matrix.
+     */
+    template <typename T>
+    void read_body_graphblas_array(std::istream &instream,
+                                   matrix_market_header &header,
+                                   GrB_Matrix mat,
+                                   const read_options& options) {
+#if FMM_GXB_PACK_UNPACK
+        // Allocate the dense matrix. GraphBLAS takes ownership of this pointer.
+        auto arr = (T*)calloc( header.nrows*header.ncols, sizeof(T));
+
+        auto handler = dense_adding_parse_handler(arr, col_major, header.nrows, header.ncols);
+        read_matrix_market_body(instream, header, handler, 1, options);
+        ok(GxB_Matrix_pack_FullC(mat, (void**)&arr, (header.nrows*header.ncols*sizeof(T)), false, nullptr));
+#else
+        // No efficient way to read dense arrays in standard GraphBLAS.
+        read_body_graphblas_coordinate<T>(instream, header, mat, options);
+#endif
+    }
+
+    /**
+     * Read a Matrix Market body into a GraphBLAS matrix.
+     */
+    template <typename T>
+    void read_body_graphblas(std::istream &instream,
+                             matrix_market_header &header,
+                             GrB_Matrix mat,
+                             const read_options& options) {
+        if (header.format == array) {
+            read_body_graphblas_array<T>(instream, header, mat, options);
+        } else {
+            read_body_graphblas_coordinate<T>(instream, header, mat, options);
+        }
+    }
+
+    /**
+     * Parse any '%%GraphBLAS type <ctype>' comments.
+     *
+     * If not present, use the MatrixMarket header field.
+     */
+    inline GrB_Type get_type_from_header(const matrix_market_header& header) {
+        const std::string kStructuredCommentBanner = "%GraphBLAS"; // The first % was consumed by the comment parser
+
+        GrB_Type desired_type = nullptr;
+
+        std::istringstream comment_iss(header.comment);
+        std::string line;
+        while (std::getline(comment_iss, line)) {
+            if (!fast_matrix_market::starts_with(line, kStructuredCommentBanner)) {
+                continue;
+            }
+
+            std::istringstream iss(line);
+            std::string banner, mode, type_string;
+
+            iss >> banner >> mode;
+            if (mode == "type") {
+                // The type is the C type. Note that in C the complex types contain whitespace:
+                // float complex, double complex
+                std::getline(iss, type_string);
+            } else {
+                // Unknown structured comment. Ignore.
+                continue;
+            }
+
+            desired_type = parse_type(trim(type_string));
+            break;
+        }
+
+        if (desired_type == nullptr) {
+            // The type is not specified in the file or by the caller. Choose appropriate type based on
+            // header field
+            switch (header.field) {
+                case pattern:
+                    desired_type = GrB_BOOL;
+                    break;
+                case integer:
+                    desired_type = GrB_INT64;
+                    break;
+                case real:
+                    desired_type = GrB_FP64;
+                    break;
+                case complex:
+#if FMM_GXB_COMPLEX
+                    desired_type = GxB_FC64;
+                    break;
+#else
+                    throw invalid_mm("Matrix is complex but no complex support in GraphBLAS.");
+#endif
+                default:
+                    desired_type = GrB_FP64;
+                    break;
+            }
+
+        }
+
+        return desired_type;
+    }
+
+    /**
+     * Read Matrix Market file into a sparse GraphBLAS matrix and a header struct.
+     */
+    inline void read_matrix_market_graphblas(std::istream &instream,
+                                             matrix_market_header &header,
+                                             GrB_Matrix* mat,
+                                             const read_options& options = {},
+                                             GrB_Type desired_type = nullptr) {
+        read_header(instream, header);
+
+        // Figure out the matrix type
+        if (desired_type == nullptr) {
+            desired_type = get_type_from_header(header);
+        }
+
+        // sanity check
+        if (header.symmetry == skew_symmetric && (desired_type == GrB_BOOL ||
+                desired_type == GrB_UINT8 || desired_type == GrB_UINT16 ||
+                desired_type == GrB_UINT32 || desired_type == GrB_UINT64)) {
+            throw invalid_mm("Skew-symmetric matrix cannot have an unsigned type");
+        }
+
+        ok(GrB_Matrix_new(mat, desired_type, header.nrows, header.ncols));
+
+        if (desired_type == GrB_BOOL) {
+            read_body_graphblas<bool>(instream, header, *mat, options);
+        } else if (desired_type == GrB_INT8) {
+            read_body_graphblas<int8_t>(instream, header, *mat, options);
+        } else if (desired_type == GrB_INT16) {
+            read_body_graphblas<int16_t>(instream, header, *mat, options);
+        } else if (desired_type == GrB_INT32) {
+            read_body_graphblas<int32_t>(instream, header, *mat, options);
+        } else if (desired_type == GrB_INT64) {
+            read_body_graphblas<int64_t>(instream, header, *mat, options);
+        } else if (desired_type == GrB_UINT8) {
+            read_body_graphblas<uint8_t>(instream, header, *mat, options);
+        } else if (desired_type == GrB_UINT16) {
+            read_body_graphblas<uint16_t>(instream, header, *mat, options);
+        } else if (desired_type == GrB_UINT32) {
+            read_body_graphblas<uint32_t>(instream, header, *mat, options);
+        } else if (desired_type == GrB_UINT64) {
+            read_body_graphblas<uint64_t>(instream, header, *mat, options);
+        } else if (desired_type == GrB_FP32) {
+            read_body_graphblas<float>(instream, header, *mat, options);
+        } else if (desired_type == GrB_FP64) {
+            read_body_graphblas<double>(instream, header, *mat, options);
+#if FMM_GXB_COMPLEX
+        } else if (desired_type == GxB_FC32) {
+            read_body_graphblas<std::complex<float>>(instream, header, *mat, options);
+        } else if (desired_type == GxB_FC64) {
+            read_body_graphblas<std::complex<double>>(instream, header, *mat, options);
+#endif
+        } else {
+            throw invalid_mm("Unknown matrix type.");
+        }
+    }
+
+    /**
+     * Read Matrix Market file into a GraphBLAS matrix.
+     */
+    inline void read_matrix_market_graphblas(std::istream &instream,
+                                             GrB_Matrix* mat,
+                                             const read_options& options = {},
+                                             GrB_Type desired_type = nullptr) {
+
+        matrix_market_header header;
+        read_matrix_market_graphblas(instream, header, mat, options, desired_type);
+    }
+
+    /**
+     * Write a GraphBLAS matrix to a coordinate MatrixMarket body by using GrB_Matrix_extractTuples.
+     */
+    template <typename T>
+    void write_body_graphblas_triplet(std::ostream &os,
+                                      matrix_market_header &header,
+                                      const GrB_Matrix& mat,
+                                      const write_options& options) {
+        std::vector<GrB_Index> rows(header.nnz);
+        std::vector<GrB_Index> cols(header.nnz);
+        auto vals = std::make_unique<T[]>(header.nnz); // Cannot use vector due to bool specialization
+
+        GrB_Index nvals = header.nnz;
+        auto ec = GraphBLAS_typed<T>::GrB_Matrix_extractTuples(rows.data(), cols.data(), vals.get(), &nvals, mat);
+        if (ec != GrB_SUCCESS) {
+            throw fast_matrix_market::invalid_argument("GrB_Matrix_extractTuples returned " + std::to_string(ec));
+        }
+        line_formatter<GrB_Index, T> lf(header, options);
+        auto formatter = triplet_formatter(lf,
+                                           rows.cbegin(), rows.cend(),
+                                           cols.cbegin(), cols.cend(),
+                                           vals.get(), header.field == pattern ? vals.get() : vals.get() + header.nnz);
+        write_body(os, formatter, options);
+    }
+
+#if FMM_GXB_ITERATORS
+    struct gblas_row_iter_impl {
+        typedef GrB_Matrix container;
+
+        static std::pair<GrB_Index, GrB_Index> setup(const GrB_Matrix& mat) {
+            GrB_Matrix_wait(mat, GrB_MATERIALIZE);
+
+            GrB_Index nvals;
+            GrB_Matrix_nvals(&nvals, mat);
+
+            GxB_Iterator iterator = attach(mat);
+            GrB_Index kount = GxB_rowIterator_kount(iterator);
+            GxB_Iterator_free(&iterator);
+
+            return std::make_pair(kount, nvals);
+        }
+
+        static GxB_Iterator attach(const GrB_Matrix& mat) {
+            GxB_Iterator iterator;
+            ok(GxB_Iterator_new(&iterator));
+            GrB_Info info = GxB_rowIterator_attach(iterator, mat, nullptr);
+            if (info < 0) {
+                throw fmm_error("GxB_rowIterator_attach returned: " + std::to_string(info));
+            }
+            return iterator;
+        }
+
+        static GrB_Info kseek(GxB_Iterator iterator, GrB_Index k) {
+            return GxB_rowIterator_kseek(iterator, k);
+        }
+
+        static GrB_Index getMajorIndex(GxB_Iterator iterator) {
+            return GxB_rowIterator_getRowIndex(iterator);
+        }
+
+        static GrB_Index getMinorIndex(GxB_Iterator iterator) {
+            return GxB_rowIterator_getColIndex(iterator);
+        }
+
+        static GrB_Info nextMinor(GxB_Iterator iterator) {
+            return GxB_rowIterator_nextCol(iterator);
+        }
+
+        static GrB_Info nextMajor(GxB_Iterator iterator) {
+            return GxB_rowIterator_nextRow(iterator);
+        }
+
+        static std::pair<GrB_Index, GrB_Index> to_row_col(GrB_Index major, GrB_Index minor) {
+            return {major, minor};
+        }
+    };
+
+    struct gblas_col_iter_impl {
+        typedef GrB_Matrix container;
+
+        static std::pair<GrB_Index, GrB_Index> setup(const GrB_Matrix& mat) {
+            GrB_Matrix_wait(mat, GrB_MATERIALIZE);
+
+            GrB_Index nvals;
+            GrB_Matrix_nvals(&nvals, mat);
+
+            GxB_Iterator iterator = attach(mat);
+            GrB_Index kount = GxB_colIterator_kount(iterator);
+            GxB_Iterator_free(&iterator);
+
+            return std::make_pair(kount, nvals);
+        }
+
+        static GxB_Iterator attach(const GrB_Matrix& mat) {
+            GxB_Iterator iterator;
+            ok(GxB_Iterator_new(&iterator));
+            GrB_Info info = GxB_colIterator_attach(iterator, mat, nullptr);
+            if (info < 0) {
+                throw fmm_error("GxB_colIterator_attach returned: " + std::to_string(info));
+            }
+            return iterator;
+        }
+
+        static GrB_Info kseek(GxB_Iterator iterator, GrB_Index k) {
+            return GxB_colIterator_kseek(iterator, k);
+        }
+
+        static GrB_Index getMinorIndex(GxB_Iterator iterator) {
+            return GxB_colIterator_getRowIndex(iterator);
+        }
+
+        static GrB_Index getMajorIndex(GxB_Iterator iterator) {
+            return GxB_colIterator_getColIndex(iterator);
+        }
+
+        static GrB_Info nextMinor(GxB_Iterator iterator) {
+            return GxB_colIterator_nextRow(iterator);
+        }
+
+        static GrB_Info nextMajor(GxB_Iterator iterator) {
+            return GxB_colIterator_nextCol(iterator);
+        }
+
+        static std::pair<GrB_Index, GrB_Index> to_row_col(GrB_Index major, GrB_Index minor) {
+            return {minor, major};
+        }
+    };
+
+#if FMM_GXB_VECTOR_ITERATORS
+    struct gblas_vec_iter_impl {
+        typedef GrB_Vector container;
+
+        static std::pair<GrB_Index, GrB_Index> setup(const GrB_Vector& vec) {
+            GrB_Vector_wait(vec, GrB_MATERIALIZE);
+
+            GrB_Index nvals;
+            GrB_Vector_nvals(&nvals, vec);
+
+            // getpmax() sometimes returns a value > nvals. No idea how that makes sense.
+            // reproduced with a dense vector with a zero value.
+//            GxB_Iterator iterator = attach(vec);
+//            GrB_Index kount = GxB_Vector_Iterator_getpmax(iterator);
+//            GxB_Iterator_free(&iterator);
+
+            return std::make_pair(nvals, nvals);
+        }
+        static GxB_Iterator attach(const GrB_Vector& vec) {
+            GxB_Iterator iterator;
+            ok(GxB_Iterator_new(&iterator));
+            GrB_Info info = GxB_Vector_Iterator_attach(iterator, vec, nullptr);
+            if (info < 0) {
+                throw fmm_error("GxB_Vector_Iterator_attach returned: " + std::to_string(info));
+            }
+            return iterator;
+        }
+
+        static GrB_Info kseek(GxB_Iterator iterator, GrB_Index k) {
+            return GxB_Vector_Iterator_seek(iterator, k);
+        }
+
+        static GrB_Index getMinorIndex([[maybe_unused]] GxB_Iterator iterator) {
+            return 0; // vectors don't have a minor axis
+        }
+
+        static GrB_Index getMajorIndex(GxB_Iterator iterator) {
+            return GxB_Vector_Iterator_getIndex(iterator);
+        }
+
+        static GrB_Info nextMinor([[maybe_unused]] GxB_Iterator iterator) {
+            return GrB_INVALID_VALUE; // vectors don't have a minor axis
+        }
+
+        static GrB_Info nextMajor(GxB_Iterator iterator) {
+            return GxB_Vector_Iterator_next(iterator);
+        }
+
+        static std::pair<GrB_Index, GrB_Index> to_row_col(GrB_Index major, GrB_Index minor) {
+            return {major, minor};
+        }
+    };
+#endif
+
+    /**
+    * Format GrB_Matrix using iterator.
+    */
+    template<typename LF, typename T, typename IMPL>
+    class GrB_Matrix_Iterator_formatter {
+    public:
+        explicit GrB_Matrix_Iterator_formatter(LF lf, const typename IMPL::container& mat) : line_formatter(lf), mat(mat) {
+            auto pair = IMPL::setup(mat);
+            kount = pair.first;
+            GrB_Index nvals = pair.second;
+            nnz_per_kount = ((double)nvals) / (double)kount;
+        }
+
+        [[nodiscard]] bool has_next() const {
+            return kount_iter < kount;
+        }
+
+        class chunk {
+        public:
+            explicit chunk(LF lf, const typename IMPL::container& mat, GrB_Index kount_iter, GrB_Index kount_end) :
+                    line_formatter(lf), mat(mat), kount_iter(kount_iter), kount_end(kount_end) {}
+
+            std::string operator()() {
+                std::string chunk;
+                chunk.reserve((kount_end - kount_iter)*250);
+
+                GxB_Iterator iterator = IMPL::attach(mat);
+
+                // seek to first assigned row or column
+                GrB_Info info = IMPL::kseek(iterator, kount_iter);
+                for (; info != GxB_EXHAUSTED && kount_iter < kount_end; ++kount_iter) {
+                    GrB_Index major = IMPL::getMajorIndex (iterator);
+
+                    // iterate over the columns in this row
+                    while (info == GrB_SUCCESS) {
+                        // get the entry A(i,j)
+                        GrB_Index minor = IMPL::getMinorIndex(iterator);
+                        auto[row, col] = IMPL::to_row_col(major, minor);
+
+                        const T& value = GraphBLAS_typed<T>::GxB_Iterator_get(iterator);
+                        chunk += line_formatter.coord_matrix(row, col, value);
+
+                        // move to the next entry
+                        info = IMPL::nextMinor(iterator);
+                    }
+                    // move to the next row or column
+                    info = IMPL::nextMajor(iterator);
+                }
+                GxB_Iterator_free(&iterator);
+
+                return chunk;
+            }
+
+            LF line_formatter;
+            const typename IMPL::container& mat;
+            GrB_Index kount_iter, kount_end;
+        };
+
+        chunk next_chunk(const write_options& options) {
+            auto kount_this_chunk = (GrB_Index)(nnz_per_kount * (double)options.chunk_size_values + 1);
+            kount_this_chunk = std::min(kount_this_chunk, kount - kount_iter);
+
+            GrB_Index kount_end = kount_iter + kount_this_chunk;
+            chunk c(line_formatter, mat, kount_iter, kount_end);
+            kount_iter = kount_end;
+
+            return c;
+        }
+
+    protected:
+        LF line_formatter;
+        const typename IMPL::container& mat;
+        double nnz_per_kount;
+        GrB_Index kount_iter = 0;
+        GrB_Index kount;
+    };
+
+#if 0
+    /**
+     * Write a GraphBLAS matrix to an array MatrixMarket body by using a column GrB_Iterator.
+     */
+    template <typename T>
+    void write_body_graphblas_array_iterator(std::ostream &os,
+                                             matrix_market_header &header,
+                                             const GrB_Matrix& mat,
+                                             const write_options& options) {
+        throw fmm_error("GraphBLAS throws GrB_NOT_IMPLEMENTED.");
+        line_formatter<GrB_Index, T> lf(header, options);
+        auto formatter = GrB_Matrix_Iterator_formatter<decltype(lf), T, gblas_col_iter_impl>(lf, mat);
+        write_body(os, formatter, options);
+    }
+#endif
+
+    /**
+     * Write a GraphBLAS matrix to a coordinate MatrixMarket body by using GrB_Iterator.
+     * This method is zero-copy, but uses a GraphBLAS extension.
+     */
+    template <typename T>
+    void write_body_graphblas_iterator(std::ostream &os,
+                                       matrix_market_header &header,
+                                       const GrB_Matrix& mat,
+                                       const write_options& options) {
+        line_formatter<GrB_Index, T> lf(header, options);
+
+        int32_t format;
+        ok(GxB_Matrix_Option_get_INT32(mat, GxB_FORMAT, &format));
+
+        if (format == GxB_BY_ROW) {
+            auto formatter = GrB_Matrix_Iterator_formatter<decltype(lf), T, gblas_row_iter_impl>(lf, mat);
+            write_body(os, formatter, options);
+        } else if (format == GxB_BY_COL) {
+            auto formatter = GrB_Matrix_Iterator_formatter<decltype(lf), T, gblas_col_iter_impl>(lf, mat);
+            write_body(os, formatter, options);
+        } else {
+            // shouldn't happen
+            write_body_graphblas_triplet<T>(os, header, mat, options);
+        }
+    }
+#endif
+
+    /**
+     * Write a GraphBLAS matrix to a MatrixMarket body.
+     */
+    template <typename T>
+    void write_body_graphblas(std::ostream &os,
+                              matrix_market_header &header,
+                              const GrB_Matrix& mat,
+                              const write_options& options) {
+        if (header.format != coordinate) {
+            throw fmm_error("FMM bug: Only coordinate format implemented.");
+        }
+
+#if FMM_GXB_ITERATORS
+        if (header.format == coordinate) {
+            write_body_graphblas_iterator<T>(os, header, mat, options);
+        } // else {
+//            write_body_graphblas_array_iterator<T>(os, header, mat, options);
+//        }
+#else
+        write_body_graphblas_triplet<T>(os, header, mat, options);
+#endif
+    }
+
+    /**
+     * Deduce Matrix Market header field type from GrB_Type.
+     */
+    inline field_type get_field(GrB_Type type) {
+        if (type == GrB_FP32 || type == GrB_FP64) {
+            return real;
+#if FMM_GXB_COMPLEX
+        } else if (type == GxB_FC32 || type == GxB_FC64) {
+            return complex;
+#endif
+        } else if (type == GrB_BOOL   || type == GrB_INT8   || type == GrB_INT16  ||
+                   type == GrB_INT32  || type == GrB_INT64  || type == GrB_UINT8  ||
+                   type == GrB_UINT16 || type == GrB_UINT32 || type == GrB_UINT64) {
+            return integer;
+        } else {
+            throw fast_matrix_market::invalid_argument("Unable to deduce Matrix Market type from GrB_Type type.");
+        }
+    }
+
+    /**
+     * Prepend a '%%GraphBLAS type <ctype>' structured comment to the actual comment
+     */
+    inline void add_structured_comment(matrix_market_header& header, GrB_Type type, const std::string& str_type) {
+        if (type == nullptr) {
+            return;
+        }
+
+        std::string header_type;
+        auto type_map = GrB_Type_to_header_type();
+        auto type_map_iter = type_map.find(type);
+        if (type_map_iter != type_map.end()) {
+            // known type
+            header_type = type_map_iter->second;
+        } else {
+            // unknown type, use the type string
+            if (str_type.empty()) {
+                return;
+            } else {
+                header_type = str_type;
+            }
+        }
+        std::string newline = header.comment.empty() ? "" : "\n";
+        header.comment = std::string("%GraphBLAS type ") + header_type + newline + header.comment;
+    }
+
+    /**
+     * Write a GraphBLAS matrix to MatrixMarket.
+     */
+    inline void write_matrix_market_graphblas(std::ostream &os,
+                                              const GrB_Matrix& mat,
+                                              const write_options& options = {}, matrix_market_header header = {}) {
+        GrB_Type type;
+        {
+            GrB_Index temp;
+
+            ok(GrB_Matrix_nrows(&temp, mat));
+            header.nrows = (int64_t)temp;
+
+            ok(GrB_Matrix_ncols(&temp, mat));
+            header.ncols = (int64_t)temp;
+
+            ok(GrB_Matrix_nvals(&temp, mat));
+            header.nnz = (int64_t)temp;
+
+            // get the matrix type
+#if FMM_GXB_TYPE_NAME
+            std::string mat_type(GxB_MAX_NAME_LEN, '\0');
+            ok(GxB_Matrix_type_name(mat_type.data(), mat));
+            mat_type.resize(mat_type.find('\0')); // truncate to the actual type name length
+            type = parse_type(mat_type);
+
+            // Write the '%%GraphBLAS type <ctype>' line
+            add_structured_comment(header, type, mat_type);
+#else
+            type = GrB_FP64; // most compatible type from the standard set.
+#endif
+        }
+
+        header.object = matrix;
+        if (header.field != pattern) {
+            header.field = get_field(type);
+        }
+        header.format = coordinate;
+#if FMM_GXB_ITERATORS
+        // GraphBLAS throws GrB_NOT_IMPLEMENTED.
+//        if (header.nnz == header.nrows * header.ncols) {
+//            // Full matrix. Only efficient way to write arrays is using a column iterator, which requires a GraphBLAS extension.
+//            header.format = array;
+//        }
+#endif
+
+        write_header(os, header);
+
+        if (type == GrB_BOOL) {
+            write_body_graphblas<bool>(os, header, mat, options);
+        } else if (type == GrB_INT8) {
+            write_body_graphblas<int8_t>(os, header, mat, options);
+        } else if (type == GrB_INT16) {
+            write_body_graphblas<int16_t>(os, header, mat, options);
+        } else if (type == GrB_INT32) {
+            write_body_graphblas<int32_t>(os, header, mat, options);
+        } else if (type == GrB_INT64) {
+            write_body_graphblas<int64_t>(os, header, mat, options);
+        } else if (type == GrB_UINT8) {
+            write_body_graphblas<uint8_t>(os, header, mat, options);
+        } else if (type == GrB_UINT16) {
+            write_body_graphblas<uint16_t>(os, header, mat, options);
+        } else if (type == GrB_UINT32) {
+            write_body_graphblas<uint32_t>(os, header, mat, options);
+        } else if (type == GrB_UINT64) {
+            write_body_graphblas<uint64_t>(os, header, mat, options);
+        } else if (type == GrB_FP32) {
+            write_body_graphblas<float>(os, header, mat, options);
+        } else if (type == GrB_FP64) {
+            write_body_graphblas<double>(os, header, mat, options);
+#if FMM_GXB_COMPLEX
+        } else if (type == GxB_FC32) {
+            write_body_graphblas<std::complex<float>>(os, header, mat, options);
+        } else if (type == GxB_FC64) {
+            write_body_graphblas<std::complex<double>>(os, header, mat, options);
+#endif
+        } else {
+            throw invalid_mm("Unknown matrix type.");
+        }
+    }
+
+
+    //////////////////////////////////////////////////////
+    // Vector methods
+
+
+    /**
+     * Read a Matrix Market body into a GraphBLAS vector.
+     */
+    template <typename T>
+    void read_body_graphblas(std::istream &instream,
+                             matrix_market_header &header,
+                             GrB_Vector vec,
+                             const read_options& options) {
+        // Read into doublets
+        std::vector<GrB_Index> inds(header.nnz);
+        // Allocate values. Cannot use std::vector due to bool specialization
+        auto vals = std::make_unique<T[]>(header.nnz);
+
+        auto handler = doublet_parse_handler(inds.begin(), vals.get());
+        read_matrix_market_body(instream, header, handler, pattern_default_value(static_cast<T*>(nullptr)), options);
+
+        // build vector from doublets
+        ok(GraphBLAS_typed<T>::build_vector(vec, inds.data(), vals.get(), header.nnz));
+    }
+
+    /**
+     * Read Matrix Market file into a sparse GraphBLAS vector and a header struct.
+     */
+    inline void read_matrix_market_graphblas(std::istream &instream,
+                                             matrix_market_header &header,
+                                             GrB_Vector* vec,
+                                             const read_options& options = {},
+                                             GrB_Type desired_type = nullptr) {
+        read_header(instream, header);
+
+        // Figure out the matrix type
+        if (desired_type == nullptr) {
+            desired_type = get_type_from_header(header);
+        }
+
+        // sanity check
+        if (header.symmetry != general) {
+            throw invalid_mm("Vectors cannot have non-general symmetry.");
+        }
+        if (std::min(header.nrows, header.ncols) > 1) {
+            throw invalid_argument("Cannot load a matrix into a vector structure");
+        }
+
+        ok(GrB_Vector_new(vec, desired_type, header.vector_length));
+
+        if (desired_type == GrB_BOOL) {
+            read_body_graphblas<bool>(instream, header, *vec, options);
+        } else if (desired_type == GrB_INT8) {
+            read_body_graphblas<int8_t>(instream, header, *vec, options);
+        } else if (desired_type == GrB_INT16) {
+            read_body_graphblas<int16_t>(instream, header, *vec, options);
+        } else if (desired_type == GrB_INT32) {
+            read_body_graphblas<int32_t>(instream, header, *vec, options);
+        } else if (desired_type == GrB_INT64) {
+            read_body_graphblas<int64_t>(instream, header, *vec, options);
+        } else if (desired_type == GrB_UINT8) {
+            read_body_graphblas<uint8_t>(instream, header, *vec, options);
+        } else if (desired_type == GrB_UINT16) {
+            read_body_graphblas<uint16_t>(instream, header, *vec, options);
+        } else if (desired_type == GrB_UINT32) {
+            read_body_graphblas<uint32_t>(instream, header, *vec, options);
+        } else if (desired_type == GrB_UINT64) {
+            read_body_graphblas<uint64_t>(instream, header, *vec, options);
+        } else if (desired_type == GrB_FP32) {
+            read_body_graphblas<float>(instream, header, *vec, options);
+        } else if (desired_type == GrB_FP64) {
+            read_body_graphblas<double>(instream, header, *vec, options);
+#if FMM_GXB_COMPLEX
+        } else if (desired_type == GxB_FC32) {
+            read_body_graphblas<std::complex<float>>(instream, header, *vec, options);
+        } else if (desired_type == GxB_FC64) {
+            read_body_graphblas<std::complex<double>>(instream, header, *vec, options);
+#endif
+        } else {
+            throw invalid_mm("Unknown vector type.");
+        }
+    }
+
+    /**
+     * Read Matrix Market file into a GraphBLAS vector.
+     */
+    inline void read_matrix_market_graphblas(std::istream &instream,
+                                             GrB_Vector* vec,
+                                             const read_options& options = {},
+                                             GrB_Type desired_type = nullptr) {
+
+        matrix_market_header header;
+        read_matrix_market_graphblas(instream, header, vec, options, desired_type);
+    }
+
+    /**
+     * Write a GraphBLAS vector to a coordinate MatrixMarket body by using GrB_Vector_extractTuples.
+     */
+    template <typename T>
+    void write_body_graphblas_doublet(std::ostream &os,
+                                      matrix_market_header &header,
+                                      const GrB_Vector& vec,
+                                      const write_options& options) {
+        std::vector<GrB_Index> indices(header.nnz);
+        auto vals = std::make_unique<T[]>(header.nnz); // Cannot use vector due to bool specialization
+
+        GrB_Index nvals = header.nnz;
+        ok(GraphBLAS_typed<T>::GrB_Vector_extractTuples(indices.data(), vals.get(), &nvals, vec));
+
+        // write the body
+        vector_line_formatter<GrB_Index, T> lf(header, options);
+        auto formatter = triplet_formatter(lf,
+                                           indices.begin(), indices.end(),
+                                           indices.begin(), indices.end(),
+                                           vals.get(), vals.get() + header.nnz);
+        write_body(os, formatter, options);
+    }
+
+#if FMM_GXB_ITERATORS && FMM_GXB_VECTOR_ITERATORS
+    /**
+     * Write a GraphBLAS vector to a coordinate MatrixMarket body by using GrB_Iterator.
+     * This method is zero-copy, but uses a GraphBLAS extension.
+     */
+    template <typename T>
+    void write_body_graphblas_iterator(std::ostream &os,
+                                       matrix_market_header &header,
+                                       const GrB_Vector& vec,
+                                       const write_options& options) {
+        vector_line_formatter<GrB_Index, T> lf(header, options);
+        auto formatter = GrB_Matrix_Iterator_formatter<decltype(lf), T, gblas_vec_iter_impl>(lf, vec);
+        write_body(os, formatter, options);
+    }
+#endif
+
+    /**
+     * Write a GraphBLAS vector to a MatrixMarket body
+     */
+    template <typename T>
+    void write_body_graphblas(std::ostream &os,
+                              matrix_market_header &header,
+                              const GrB_Vector& vec,
+                              const write_options& options) {
+#if FMM_GXB_ITERATORS && FMM_GXB_VECTOR_ITERATORS
+        write_body_graphblas_iterator<T>(os, header, vec, options);
+#else
+        write_body_graphblas_doublet<T>(os, header, vec, options);
+#endif
+    }
+
+    /**
+     * Write a GraphBLAS vector to MatrixMarket.
+     */
+    inline void write_matrix_market_graphblas(std::ostream &os,
+                                              const GrB_Vector& vec,
+                                              const write_options& options = {}, matrix_market_header header = {}) {
+        GrB_Type type;
+        {
+            GrB_Index temp;
+
+            ok(GrB_Vector_size(&temp, vec));
+            header.vector_length = (int64_t)temp;
+
+            ok(GrB_Vector_nvals(&temp, vec));
+            header.nnz = (int64_t)temp;
+
+            // get the matrix type
+#if FMM_GXB_TYPE_NAME
+            std::string mat_type(GxB_MAX_NAME_LEN, '\0');
+            ok(GxB_Vector_type_name(mat_type.data(), vec));
+            mat_type.resize(mat_type.find('\0')); // truncate to the actual type name length
+            type = parse_type(mat_type);
+
+            // Write the '%%GraphBLAS type <ctype>' line
+            add_structured_comment(header, type, mat_type);
+#else
+            type = GrB_FP64; // most compatible type from the standard set.
+#endif
+        }
+
+        header.object = vector;
+        if (header.field != pattern) {
+            header.field = get_field(type);
+        }
+        header.format = coordinate;
+
+        write_header(os, header);
+
+        if (type == GrB_BOOL) {
+            write_body_graphblas<bool>(os, header, vec, options);
+        } else if (type == GrB_INT8) {
+            write_body_graphblas<int8_t>(os, header, vec, options);
+        } else if (type == GrB_INT16) {
+            write_body_graphblas<int16_t>(os, header, vec, options);
+        } else if (type == GrB_INT32) {
+            write_body_graphblas<int32_t>(os, header, vec, options);
+        } else if (type == GrB_INT64) {
+            write_body_graphblas<int64_t>(os, header, vec, options);
+        } else if (type == GrB_UINT8) {
+            write_body_graphblas<uint8_t>(os, header, vec, options);
+        } else if (type == GrB_UINT16) {
+            write_body_graphblas<uint16_t>(os, header, vec, options);
+        } else if (type == GrB_UINT32) {
+            write_body_graphblas<uint32_t>(os, header, vec, options);
+        } else if (type == GrB_UINT64) {
+            write_body_graphblas<uint64_t>(os, header, vec, options);
+        } else if (type == GrB_FP32) {
+            write_body_graphblas<float>(os, header, vec, options);
+        } else if (type == GrB_FP64) {
+            write_body_graphblas<double>(os, header, vec, options);
+#if FMM_GXB_COMPLEX
+        } else if (type == GxB_FC32) {
+            write_body_graphblas<std::complex<float>>(os, header, vec, options);
+        } else if (type == GxB_FC64) {
+            write_body_graphblas<std::complex<double>>(os, header, vec, options);
+#endif
+        } else {
+            throw invalid_mm("Unknown vector type.");
+        }
+    }
+}
diff --git a/include/fast_matrix_market/fast_matrix_market.hpp b/include/fast_matrix_market/fast_matrix_market.hpp
index 8aa0daa..68d0b67 100644
--- a/include/fast_matrix_market/fast_matrix_market.hpp
+++ b/include/fast_matrix_market/fast_matrix_market.hpp
@@ -6,6 +6,7 @@
 
 #pragma once
 
+#include <algorithm>
 #include <complex>
 #include <future>
 #include <iostream>
@@ -193,12 +194,36 @@ namespace fast_matrix_market {
         return (flags & flag) == flag;
     }
 
+    inline bool starts_with(const std::string &str, const std::string& prefix) {
+        if (prefix.size() > str.size()) {
+            return false;
+        }
+        return std::equal(prefix.begin(), prefix.end(), str.begin());
+    }
+
     inline bool ends_with(const std::string &str, const std::string& suffix) {
         if (suffix.size() > str.size()) {
             return false;
         }
         return std::equal(suffix.rbegin(), suffix.rend(), str.rbegin());
     }
+
+    /**
+     * Trim the whitespace from both ends of a string. Returns a copy.
+     */
+    inline std::string trim(std::string s) {
+        // ltrim
+        s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](unsigned char ch) {
+            return !std::isspace(ch);
+        }));
+
+        // rtrim
+        s.erase(std::find_if(s.rbegin(), s.rend(), [](unsigned char ch) {
+            return !std::isspace(ch);
+        }).base(), s.end());
+
+        return s;
+    }
 }
 
 #include "field_conv.hpp"
diff --git a/include/fast_matrix_market/formatters.hpp b/include/fast_matrix_market/formatters.hpp
index ab07335..eab8e6b 100644
--- a/include/fast_matrix_market/formatters.hpp
+++ b/include/fast_matrix_market/formatters.hpp
@@ -20,6 +20,10 @@ namespace fast_matrix_market {
                                                                                            options(options) {}
 
         std::string coord_matrix(const IT& row, const IT& col, const VT& val) {
+            if (header.format == array) {
+                return array_matrix(row, col, val);
+            }
+
             std::string line{};
             line += int_to_string(row + 1);
             line += kSpace;
diff --git a/include/fast_matrix_market/header.hpp b/include/fast_matrix_market/header.hpp
index 4e623b5..cbf9e36 100644
--- a/include/fast_matrix_market/header.hpp
+++ b/include/fast_matrix_market/header.hpp
@@ -195,7 +195,12 @@ namespace fast_matrix_market {
                 } else {
                     header.nnz = header.nrows * header.ncols;
                 }
-                header.vector_length = -1;
+                if (std::min(header.nrows, header.ncols) == 1) {
+                    // row or column matrix. Either one can be loaded into a vector data structure.
+                    header.vector_length = std::max(header.nrows, header.ncols);
+                } else {
+                    header.vector_length = -1;
+                }
             }
         }
 


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index c224f8a..8ffe9f5 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -21,9 +21,23 @@ jobs:
       with:
         python-version: ${{ matrix.python-version }}
 
+# The GraphBLAS that is in the Ubuntu repos is so old it doesn't even have GrB_Scalar.
+#    - name: Install GraphBLAS (Ubuntu)
+#      if: contains(matrix.os, 'ubuntu')
+#        # install suitesparse-dev which installs GraphBLAS.h and libgraphblas.a and .so
+#        # These are very old, so disable all the fun features.
+#      run: |
+#        sudo apt-get install -y libsuitesparse-dev
+#        echo "FMM_CMAKE_FLAGS=-DFMM_NO_GXB=ON" >> $GITHUB_ENV
+
+    - name: Install GraphBLAS (macOS)
+      if: contains(matrix.os, 'macos')
+        # install suite-sparse which includes GraphBLAS and also FindGraphBLAS.cmake
+      run: brew install suite-sparse
+
     - name: Configure
       run: |
-        cmake -S . -B build -DCMAKE_BUILD_TYPE=Debug -DFAST_MATRIX_MARKET_TEST=ON -DFAST_MATRIX_MARKET_TEST_COVERAGE=ON
+        cmake -S . -B build -DCMAKE_BUILD_TYPE=Debug -DFAST_MATRIX_MARKET_TEST=ON -DFAST_MATRIX_MARKET_TEST_COVERAGE=ON ${{ env.FMM_CMAKE_FLAGS }}
 
     - name: Build
       run: |
@@ -60,4 +74,4 @@ jobs:
       with:
         gcov: true
         gcov_include: include/*   # C++ coverage
-#        directory: python/tests   # Python coverage. This seems to confuse codedov, disable for now.
\ No newline at end of file
+#        directory: python/tests   # Python coverage. This seems to confuse codecov, disable for now.
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index e60fc10..2c87da1 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -80,4 +80,34 @@ if (FAST_MATRIX_MARKET_TEST_EXTERNAL_APPS)
     target_include_directories(blaze_test PUBLIC ${BLAZE_INCLUDE_DIR})
 
     gtest_discover_tests(blaze_test)
+
+    # GraphBLAS
+    option(FAST_MATRIX_MARKET_TEST_GRAPHBLAS "Enable tests for GraphBLAS. Requires GraphBLAS installed on system such that find_package will find it." ON)
+    if (FAST_MATRIX_MARKET_TEST_GRAPHBLAS)
+        include(../cmake/GraphBLAS.cmake)
+
+        if (GraphBLAS_FOUND)
+            message("GraphBLAS found.")
+            message("GRAPHBLAS_INCLUDE_DIR: ${GRAPHBLAS_INCLUDE_DIR}")
+            message("GRAPHBLAS_LIBRARIES: ${GRAPHBLAS_LIBRARIES}")
+            message("GRAPHBLAS_STATIC: ${GRAPHBLAS_STATIC}")
+            message("GRAPHBLAS_LIBRARY: ${GRAPHBLAS_LIBRARY}")
+
+            add_executable(graphblas_test graphblas_test.cpp)
+
+            if (NOT ("${GRAPHBLAS_INCLUDE_DIR}" STREQUAL "" ))
+                target_include_directories(graphblas_test PUBLIC ${GRAPHBLAS_INCLUDE_DIR})
+            endif()
+            target_link_libraries(graphblas_test GTest::gtest_main fast_matrix_market::fast_matrix_market ${GRAPHBLAS_LIBRARIES})
+
+            option(FMM_NO_GXB "Disable GxB_* extensions" OFF)
+            if (FMM_NO_GXB)
+                target_compile_definitions(graphblas_test PUBLIC FMM_NO_GXB)
+            endif()
+
+            gtest_discover_tests(graphblas_test)
+        else()
+            message("GraphBLAS not found, skipping GraphBLAS tests.")
+        endif()
+    endif()
 endif()
diff --git a/tests/basic_test.cpp b/tests/basic_test.cpp
index 2c16fc2..b6d5d5b 100644
--- a/tests/basic_test.cpp
+++ b/tests/basic_test.cpp
@@ -154,6 +154,16 @@ TEST(Utils, misc) {
     EXPECT_TRUE(fast_matrix_market::ends_with("foobar", "bar"));
     EXPECT_FALSE(fast_matrix_market::ends_with("", "bar"));
 
+    EXPECT_FALSE(fast_matrix_market::starts_with("foo", "bar"));
+    EXPECT_TRUE(fast_matrix_market::starts_with("foobar", "foo"));
+    EXPECT_FALSE(fast_matrix_market::starts_with("", "bar"));
+
+    EXPECT_EQ(fast_matrix_market::trim("foo"), "foo");
+    EXPECT_EQ(fast_matrix_market::trim(" foo"), "foo");
+    EXPECT_EQ(fast_matrix_market::trim(" \nfoo"), "foo");
+    EXPECT_EQ(fast_matrix_market::trim("foo "), "foo");
+    EXPECT_EQ(fast_matrix_market::trim("foo\n  \n"), "foo");
+
     std::string msg("error");
     EXPECT_EQ(fast_matrix_market::fmm_error(msg).what(), msg);
 }
diff --git a/tests/fmm_tests.hpp b/tests/fmm_tests.hpp
index e84c770..513dc80 100644
--- a/tests/fmm_tests.hpp
+++ b/tests/fmm_tests.hpp
@@ -49,8 +49,39 @@ struct array_matrix {
             return vals[col * nrows + row];
         }
     }
+
+    typename std::vector<VT>::const_reference operator()(int64_t row, int64_t col) const {
+        if (order == fast_matrix_market::row_major) {
+            return vals[row * ncols + col];
+        } else {
+            return vals[col * nrows + row];
+        }
+    }
 };
 
+template <typename FT>
+std::ostream& operator<<(std::ostream& os, const std::complex<FT>& c) {
+    os << c.real();
+    if (c.imag() != 0) {
+        os << "+" << c.imag() << "i";
+    }
+    return os;
+}
+
+template <typename VT>
+std::ostream& operator<<(std::ostream& os, const array_matrix<VT>& arr) {
+    std::string order = (arr.order == fast_matrix_market::row_major ? "row_major" : "col_major");
+    os << arr.nrows << "-by-" << arr.ncols << " " << order << " dense array" << std::endl;
+
+    for (int64_t row = 0; row < arr.nrows; ++row) {
+        for (int64_t col = 0; col < arr.ncols; ++col) {
+            os << " " << std::setw(5) << arr(row, col);
+        }
+        os << std::endl;
+    }
+    return os;
+}
+
 /**
  * Sort the elements of a triplet matrix for consistent equality checks.
  */
diff --git a/tests/graphblas_test.cpp b/tests/graphblas_test.cpp
new file mode 100644
index 0000000..a13d7cc
--- /dev/null
+++ b/tests/graphblas_test.cpp
@@ -0,0 +1,571 @@
+// Copyright (C) 2022-2023 Adam Lugowski. All rights reserved.
+// Use of this source code is governed by the BSD 2-clause license found in the LICENSE.txt file.
+
+#include <algorithm>
+#include <fstream>
+#include <sstream>
+
+#if defined(__clang__)
+// GraphBLAS index type, GrB_Index, is unsigned whereas all fast_matrix_market types are signed.
+//#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wsign-compare"
+#endif
+
+#include "fmm_tests.hpp"
+
+#if defined(__clang__)
+// Disable warnings from GraphBLAS headers. They must be included with extern "C" but the complex types are defined as
+// std::complex.
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wreturn-type-c-linkage"
+#endif
+
+#include <fast_matrix_market/app/GraphBLAS.hpp>
+
+#if defined(__clang__)
+#pragma clang diagnostic pop
+#endif
+
+
+/**
+ * GraphBLAS needs GrB_init() to be called before any other methods, else you get a GrB_PANIC.
+ * Use the GoogleTest Environment feature to call the setup and finalizers appropriately.
+ */
+class GraphBLASEnvironment : public ::testing::Environment {
+public:
+    ~GraphBLASEnvironment() override = default;
+
+    // Override this to define how to set up the environment.
+    void SetUp() override {
+        GrB_init(GrB_BLOCKING);
+    }
+
+    // Override this to define how to tear down the environment.
+    void TearDown() override {
+        GrB_finalize();
+    }
+};
+
+[[maybe_unused]] testing::Environment* const graphblas_env =
+        ::testing::AddGlobalTestEnvironment(new GraphBLASEnvironment);
+
+
+GrB_Matrix read_mtx(const std::string& source) {
+    fast_matrix_market::read_options options;
+    options.chunk_size_bytes = 1;
+
+    std::istringstream iss(source);
+    GrB_Matrix ret;
+    fast_matrix_market::read_matrix_market_graphblas(iss, &ret, options);
+    return ret;
+}
+
+GrB_Vector read_mtx_vec(const std::string& source) {
+    fast_matrix_market::read_options options;
+    options.chunk_size_bytes = 1;
+
+    std::istringstream iss(source);
+    GrB_Vector ret;
+    fast_matrix_market::read_matrix_market_graphblas(iss, &ret, options);
+    return ret;
+}
+
+template <typename MAT_OR_VEC>
+std::string write_mtx(MAT_OR_VEC mat, bool pattern_only=false) {
+    fast_matrix_market::write_options options;
+    options.chunk_size_values = 1;
+
+    fast_matrix_market::matrix_market_header header{};
+    if (pattern_only) header.field = fast_matrix_market::pattern;
+
+    std::ostringstream oss;
+    fast_matrix_market::write_matrix_market_graphblas(oss, mat, options, header);
+    return oss.str();
+}
+
+template <typename ARRAY>
+void read_array_file(const std::string& matrix_filename, ARRAY& array, fast_matrix_market::read_options options = {}) {
+    std::ifstream f(matrix_filename);
+    options.chunk_size_bytes = 1;
+
+    fast_matrix_market::read_matrix_market_array(f, array.nrows, array.ncols, array.vals, array.order, options);
+}
+
+GrB_Index nrows(const GrB_Matrix& m) {
+    GrB_Index ret;
+    GrB_Matrix_nrows(&ret, m);
+    return ret;
+}
+
+GrB_Index ncols(const GrB_Matrix& m) {
+    GrB_Index ret;
+    GrB_Matrix_ncols(&ret, m);
+    return ret;
+}
+
+GrB_Index nnz(const GrB_Matrix& m) {
+    GrB_Index ret;
+    GrB_Matrix_nvals(&ret, m);
+    return ret;
+}
+
+GrB_Index nnz(const GrB_Vector& m) {
+    GrB_Index ret;
+    GrB_Vector_nvals(&ret, m);
+    return ret;
+}
+
+GrB_Index size(const GrB_Vector& m) {
+    GrB_Index ret;
+    GrB_Vector_size(&ret, m);
+    return ret;
+}
+
+bool element_equals(const GrB_Matrix& lhs, const GrB_Matrix& rhs, GrB_Index row, GrB_Index col) {
+    // every tested type can be converted to std::complex<double>
+    std::complex<double> l, r;
+    auto l_ret = GxB_Matrix_extractElement_FC64(&l, lhs, row, col);
+    auto r_ret = GxB_Matrix_extractElement_FC64(&r, rhs, row, col);
+    if (l_ret != r_ret) {
+        return false;
+    }
+
+    if (l_ret == GrB_NO_VALUE) {
+        return true;
+    }
+
+    if (l_ret == GrB_SUCCESS) {
+        return l == r;
+    }
+
+    return false;
+}
+
+bool element_equals(const GrB_Vector& lhs, const GrB_Vector& rhs, GrB_Index row) {
+    // every tested type can be converted to std::complex<double>
+    std::complex<double> l, r;
+    auto l_ret = GxB_Vector_extractElement_FC64(&l, lhs, row);
+    auto r_ret = GxB_Vector_extractElement_FC64(&r, rhs, row);
+    if (l_ret != r_ret) {
+        return false;
+    }
+
+    if (l_ret == GrB_NO_VALUE) {
+        return true;
+    }
+
+    if (l_ret == GrB_SUCCESS) {
+        return l == r;
+    }
+
+    return false;
+}
+
+bool element_equals(const array_matrix<std::complex<double>>& lhs, const GrB_Matrix& rhs, GrB_Index row, GrB_Index col) {
+    // every tested type can be converted to std::complex<double>
+    std::complex<double> l, r;
+    l = lhs((int64_t)row, (int64_t)col);
+    auto r_ret = GxB_Matrix_extractElement_FC64(&r, rhs, row, col);
+
+    if (r_ret == GrB_NO_VALUE && l.real() == 0) {
+        return true;
+    }
+
+    if (r_ret == GrB_SUCCESS) {
+        return l == r;
+    }
+
+    return false;
+}
+
+/**
+ * Equality check that prints the matrices if they're unequal. Useful for debugging when remote tests fail.
+ */
+bool is_equal(const GrB_Matrix& lhs, const GrB_Matrix& rhs) {
+    if (nrows(lhs) != nrows(rhs) || ncols(lhs) != ncols(rhs)) {
+        return false;
+    }
+
+    bool equal = true;
+    for (GrB_Index row = 0; row < nrows(rhs); ++row) {
+        for (GrB_Index col = 0; col < ncols(rhs); ++col) {
+            if (!element_equals(lhs, rhs, row, col)) {
+                equal = false;
+            }
+        }
+    }
+
+    if (!equal) {
+        GxB_Matrix_fprint(lhs, "Left", GxB_COMPLETE, nullptr);
+        GxB_Matrix_fprint(rhs, "Right", GxB_COMPLETE, nullptr);
+        std::cout << std::endl << std::endl;
+    }
+
+    return equal;
+}
+
+/**
+ * Equality check that prints the matrices if they're unequal. Useful for debugging when remote tests fail.
+ */
+bool is_equal(const array_matrix<std::complex<double>>& lhs, const GrB_Matrix& rhs) {
+    if (lhs.nrows != nrows(rhs) || lhs.ncols != ncols(rhs)) {
+        return false;
+    }
+
+    bool equal = true;
+    for (GrB_Index row = 0; row < nrows(rhs); ++row) {
+        for (GrB_Index col = 0; col < ncols(rhs); ++col) {
+            if (!element_equals(lhs, rhs, row, col)) {
+                equal = false;
+            }
+        }
+    }
+
+    if (!equal) {
+        std::cout << "Left:" << std::endl << lhs << std::endl;
+        GxB_Matrix_fprint(rhs, "Right", GxB_COMPLETE, nullptr);
+        std::cout << std::endl << std::endl;
+    }
+
+    return equal;
+}
+
+/**
+ * Equality check that prints the vectors if they're unequal. Useful for debugging when remote tests fail.
+ */
+bool is_equal(const GrB_Vector& lhs, const GrB_Vector& rhs) {
+    bool equal = true;
+    if (size(lhs) != size(rhs)) {
+        equal = false;
+    }
+
+    if (equal) {
+        for (GrB_Index row = 0; row < size(rhs); ++row) {
+            if (!element_equals(lhs, rhs, row)) {
+                equal = false;
+            }
+        }
+    }
+
+    if (!equal) {
+        GxB_Vector_fprint(lhs, "Left", GxB_COMPLETE, nullptr);
+        GxB_Vector_fprint(rhs, "Right", GxB_COMPLETE, nullptr);
+        std::cout << std::endl << std::endl;
+    }
+
+    return equal;
+}
+
+
+bool element_equals(const array_matrix<std::complex<double>>& lhs, const GrB_Vector& rhs, GrB_Index row) {
+    // every tested type can be converted to std::complex<double>
+    std::complex<double> l, r;
+    l = lhs((int64_t)row, (int64_t)0);
+    auto r_ret = GxB_Vector_extractElement_FC64(&r, rhs, row);
+
+    if (r_ret == GrB_NO_VALUE && l.real() == 0) {
+        return true;
+    }
+
+    if (r_ret == GrB_SUCCESS) {
+        return l == r;
+    }
+
+    return false;
+}
+
+/**
+ * Equality check that prints the vectors if they're unequal. Useful for debugging when remote tests fail.
+ */
+bool is_equal(const array_matrix<std::complex<double>>& lhs, const GrB_Vector& rhs) {
+    bool equal = true;
+    if (lhs.nrows != size(rhs)) {
+        equal = false;
+    }
+
+    if (equal) {
+        for (GrB_Index row = 0; row < size(rhs); ++row) {
+            if (!element_equals(lhs, rhs, row)) {
+                equal = false;
+            }
+        }
+    }
+
+    if (!equal) {
+        std::cout << "Left:" << std::endl << lhs << std::endl;
+        GxB_Vector_fprint(rhs, "Right", GxB_COMPLETE, nullptr);
+        std::cout << std::endl << std::endl;
+    }
+
+    return equal;
+}
+
+class GraphBLASMatrixTest : public testing::TestWithParam<std::string> {
+public:
+    ~GraphBLASMatrixTest() override {
+        GrB_Matrix_free(&full);
+        GrB_Matrix_free(&bitmap_col);
+        GrB_Matrix_free(&bitmap_row);
+        GrB_Matrix_free(&csc);
+        GrB_Matrix_free(&csr);
+        GrB_Matrix_free(&dcsc);
+        GrB_Matrix_free(&dcsr);
+    }
+
+    void load(const std::string& param) {
+        std::string load_path = kTestMatrixDir + "/" + param;
+
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_header(f, header);
+        }
+
+#if FMM_GXB_COMPLEX
+#else
+        if (header.field == fast_matrix_market::complex) {
+            // no complex value support in this version of GraphBLAS.
+            return;
+        }
+#endif
+        // Load dense array for correctness tests
+        read_array_file(load_path, array);
+
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &full);
+            GxB_Matrix_Option_set(full, GxB_SPARSITY_CONTROL, GxB_FULL);
+        }
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &bitmap_col);
+            GxB_Matrix_Option_set(bitmap_col, GxB_SPARSITY_CONTROL, GxB_BITMAP);
+            GxB_Matrix_Option_set(bitmap_col, GxB_FORMAT, GxB_BY_COL);
+        }
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &bitmap_row);
+            GxB_Matrix_Option_set(bitmap_row, GxB_SPARSITY_CONTROL, GxB_BITMAP);
+            GxB_Matrix_Option_set(bitmap_row, GxB_FORMAT, GxB_BY_ROW);
+        }
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &csc);
+            GxB_Matrix_Option_set(csc, GxB_SPARSITY_CONTROL, GxB_SPARSE);
+            GxB_Matrix_Option_set(csc, GxB_FORMAT, GxB_BY_COL);
+        }
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &csr);
+            GxB_Matrix_Option_set(csr, GxB_SPARSITY_CONTROL, GxB_SPARSE);
+            GxB_Matrix_Option_set(csr, GxB_FORMAT, GxB_BY_ROW);
+        }
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &dcsc);
+            GxB_Matrix_Option_set(dcsc, GxB_SPARSITY_CONTROL, GxB_HYPERSPARSE);
+            GxB_Matrix_Option_set(dcsc, GxB_FORMAT, GxB_BY_COL);
+        }
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &dcsr);
+            GxB_Matrix_Option_set(dcsr, GxB_SPARSITY_CONTROL, GxB_HYPERSPARSE);
+            GxB_Matrix_Option_set(dcsr, GxB_FORMAT, GxB_BY_ROW);
+        }
+    }
+
+    protected:
+    fast_matrix_market::matrix_market_header header;
+    array_matrix<std::complex<double>> array;
+
+    GrB_Matrix full = nullptr;
+    GrB_Matrix bitmap_col = nullptr;
+    GrB_Matrix bitmap_row = nullptr;
+    GrB_Matrix csc = nullptr;
+    GrB_Matrix csr = nullptr;
+    GrB_Matrix dcsc = nullptr;
+    GrB_Matrix dcsr = nullptr;
+
+    std::map<std::string, GrB_Matrix&> sparse_mats = {
+            {"full", full},
+            {"bitmap_col", bitmap_col},
+            {"bitmap_row", bitmap_row},
+            {"csc", csc},
+            {"csr", csr},
+            {"dcsc", dcsc},
+            {"dcsr", dcsr},
+    };
+};
+
+TEST_P(GraphBLASMatrixTest, SmallMatrices) {
+    EXPECT_NO_THROW(load(GetParam()));
+
+#if FMM_GXB_COMPLEX
+#else
+    if (header.field == fast_matrix_market::complex) {
+        GTEST_SKIP() << "No complex value support in this GraphBLAS.";
+    }
+#endif
+
+    EXPECT_TRUE(is_equal(csc, csr));
+
+    // Make sure all the variants are loaded correctly.
+    for (const auto& [label, mat] : sparse_mats) {
+        EXPECT_EQ(header.nrows, nrows(mat));
+        EXPECT_EQ(header.ncols, ncols(mat));
+        EXPECT_EQ(header.nnz, nnz(mat));
+        EXPECT_TRUE(is_equal(full, mat));
+
+        // compare to non-GraphBLAS array
+        EXPECT_TRUE(is_equal(array, mat));
+
+        // Test write
+        std::string s = write_mtx(mat);
+        auto loaded = read_mtx(s);
+        EXPECT_TRUE(is_equal(mat, loaded));
+    }
+
+    // write pattern
+    std::string pattern_mtx = write_mtx(csc, true);
+    GrB_Matrix pat_mat = nullptr;
+    EXPECT_TRUE(pattern_mtx.find("pattern") > 0); // pattern should appear in the header
+    {
+        std::istringstream pat_iss(pattern_mtx);
+        fast_matrix_market::read_matrix_market_array(pat_iss, array.nrows, array.ncols, array.vals, array.order);
+    }
+    {
+        std::istringstream pat_iss(pattern_mtx);
+        fast_matrix_market::read_matrix_market_graphblas(pat_iss, &pat_mat);
+    }
+
+    EXPECT_TRUE(is_equal(array, pat_mat));
+    GrB_Matrix_free(&pat_mat);
+}
+
+INSTANTIATE_TEST_SUITE_P(
+        GraphBLASMatrixTest,
+        GraphBLASMatrixTest,
+        ::testing::Values("eye3.mtx"
+                , "row_3by4.mtx"
+                , "2col_array.mtx"
+                , "kepner_gilbert_graph.mtx"
+                , "vector_array.mtx"
+                , "vector_coordinate.mtx"
+                , "eye3_pattern.mtx"
+                , "eye3_complex.mtx"
+                , "vector_coordinate_complex.mtx"
+        ));
+
+// Test %%GraphBLAS type <ctype>
+TEST(GraphBLASMatrixTest, TypeComment) {
+#if !FMM_GXB_TYPE_NAME
+    GTEST_SKIP() << "No type interrogation support in this GraphBLAS.";
+#endif
+
+    auto filenames = std::vector{
+            "matrix_bool.mtx"
+            , "matrix_double.mtx"
+            , "matrix_double_complex.mtx"
+            , "matrix_float.mtx"
+            , "matrix_float_complex.mtx"
+            , "matrix_int16.mtx"
+            , "matrix_int32.mtx"
+            , "matrix_int64.mtx"
+            , "matrix_int8.mtx"
+            , "matrix_uint16.mtx"
+            , "matrix_uint32.mtx"
+            , "matrix_uint64.mtx"
+            , "matrix_uint8.mtx"
+    };
+    for (auto param : filenames) {
+        std::string load_path = kTestMatrixDir + "graphblas/" + param;
+        fast_matrix_market::matrix_market_header header;
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_header(f, header);
+        }
+#if FMM_GXB_COMPLEX
+#else
+        if (header.field == fast_matrix_market::complex) {
+            // No complex value support in this GraphBLAS.
+            continue;
+        }
+#endif
+        // test ctype for matrices
+        {
+            GrB_Matrix mat = nullptr;
+
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &mat);
+
+            // Test write
+            std::string s = write_mtx(mat);
+            auto loaded = read_mtx(s);
+            EXPECT_TRUE(is_equal(mat, loaded));
+
+            // Verify the written matrix has the same type written
+            std::string type_comment = header.comment;
+            EXPECT_NE(s.find(type_comment), std::string::npos);
+        }
+        // test ctype for vectors
+        {
+            GrB_Vector vec = nullptr;
+
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &vec);
+
+            // Test write
+            std::string s = write_mtx(vec);
+            auto loaded = read_mtx_vec(s);
+            EXPECT_TRUE(is_equal(vec, loaded));
+
+            // Verify the written matrix has the same type written
+            std::string type_comment = header.comment;
+            EXPECT_NE(s.find(type_comment), std::string::npos);
+        }
+   }
+}
+
+
+/////////////////////////////////////////////////////////
+///// Vectors
+/////////////////////////////////////////////////////////
+
+TEST(GraphBLASMatrixTest, Vector) {
+    auto filenames = std::vector{
+            "vector_array.mtx"
+            , "vector_coordinate.mtx"
+            , "vector_coordinate_complex.mtx"
+    };
+    for (auto param : filenames) {
+        std::string load_path = kTestMatrixDir + param;
+        fast_matrix_market::matrix_market_header header;
+        GrB_Vector vec = nullptr;
+
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_header(f, header);
+        }
+#if FMM_GXB_COMPLEX
+#else
+        if (header.field == fast_matrix_market::complex) {
+            // No complex value support in this GraphBLAS.
+            continue;
+        }
+#endif
+        array_matrix<std::complex<double>> array;
+        read_array_file(load_path, array);
+
+        {
+            std::ifstream f(load_path);
+            fast_matrix_market::read_matrix_market_graphblas(f, &vec);
+        }
+
+        // test for correct load against a dense matrix
+        EXPECT_TRUE(is_equal(array, vec));
+
+        // Test write
+        std::string s = write_mtx(vec);
+        auto loaded = read_mtx_vec(s);
+        EXPECT_TRUE(is_equal(vec, loaded));
+    }
+}
\ No newline at end of file
diff --git a/tests/matrices/2col_array.mtx b/tests/matrices/2col_array.mtx
new file mode 100644
index 0000000..f1897c7
--- /dev/null
+++ b/tests/matrices/2col_array.mtx
@@ -0,0 +1,9 @@
+%%MatrixMarket matrix array real general
+%
+3 2
+10
+11
+12
+20
+21
+22
diff --git a/tests/matrices/graphblas/matrix_bool.mtx b/tests/matrices/graphblas/matrix_bool.mtx
new file mode 100644
index 0000000..2008e47
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_bool.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type bool
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_double.mtx b/tests/matrices/graphblas/matrix_double.mtx
new file mode 100644
index 0000000..744059d
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_double.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type double
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_double_complex.mtx b/tests/matrices/graphblas/matrix_double_complex.mtx
new file mode 100644
index 0000000..fb13920
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_double_complex.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate complex general
+%%GraphBLAS type double complex
+3 1 3
+1 1 1 0
+2 1 1 0
+3 1 1 0
diff --git a/tests/matrices/graphblas/matrix_float.mtx b/tests/matrices/graphblas/matrix_float.mtx
new file mode 100644
index 0000000..e613512
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_float.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type float
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_float_complex.mtx b/tests/matrices/graphblas/matrix_float_complex.mtx
new file mode 100644
index 0000000..07e42ac
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_float_complex.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate complex general
+%%GraphBLAS type float complex
+3 1 3
+1 1 1 0
+2 1 1 0
+3 1 1 0
diff --git a/tests/matrices/graphblas/matrix_int16.mtx b/tests/matrices/graphblas/matrix_int16.mtx
new file mode 100644
index 0000000..9061e49
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_int16.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type int16_t
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_int32.mtx b/tests/matrices/graphblas/matrix_int32.mtx
new file mode 100644
index 0000000..48bdc60
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_int32.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type int32_t
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_int64.mtx b/tests/matrices/graphblas/matrix_int64.mtx
new file mode 100644
index 0000000..f323e9c
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_int64.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type int64_t
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_int8.mtx b/tests/matrices/graphblas/matrix_int8.mtx
new file mode 100644
index 0000000..12093ff
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_int8.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type int8_t
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_uint16.mtx b/tests/matrices/graphblas/matrix_uint16.mtx
new file mode 100644
index 0000000..e093da5
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_uint16.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type uint16_t
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_uint32.mtx b/tests/matrices/graphblas/matrix_uint32.mtx
new file mode 100644
index 0000000..68473c5
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_uint32.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type uint32_t
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_uint64.mtx b/tests/matrices/graphblas/matrix_uint64.mtx
new file mode 100644
index 0000000..21d7fcf
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_uint64.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type uint64_t
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/graphblas/matrix_uint8.mtx b/tests/matrices/graphblas/matrix_uint8.mtx
new file mode 100644
index 0000000..18b3352
--- /dev/null
+++ b/tests/matrices/graphblas/matrix_uint8.mtx
@@ -0,0 +1,6 @@
+%%MatrixMarket matrix coordinate real general
+%%GraphBLAS type uint8_t
+3 1 3
+1 1 1
+2 1 1
+3 1 1
diff --git a/tests/matrices/row_3by4.mtx b/tests/matrices/row_3by4.mtx
index 80e3ba3..f748e93 100644
--- a/tests/matrices/row_3by4.mtx
+++ b/tests/matrices/row_3by4.mtx
@@ -1,4 +1,4 @@
-%%MatrixMarket matrix coordinate real general
+%%MatrixMarket matrix coordinate integer general
 % 3-by-4 matrix with just one row set
 3 4 4
 1 1 1