Accelerate’s SparseMultiply randomly crashes with EXC_BAD_ACCESS

I am using SparseMultiply of a SparseMatrix_Double and a DenseVector_Double, using the Accelerate framework’s Sparse Solvers, and it is randomly crashing:

Sometimes, it works. Sometimes, it crashes with EXC_BAD_ACCESS. Sometimes, it returns with incorrect results.

I’ve tried turning on various “address sanitizer” and/or various memory-checking options on the “Diagnostics” tab of my scheme settings, and they do not report any issues.

What is going on?

func matrixProductExperiment() {
    // Given sparse matrix A, and dense vector X, calculate product of dense vector Y, i.e., y = Ax:
    //
    //              A                X    =      Y
    //   ( 10.0  1.0      2.5 )  ( 2.20 ) = ( 32.025 )
    //   (  1.0 12.0 -0.3 1.1 )  ( 2.85 ) = ( 38.720 )
    //   (      -0.3  9.5     )  ( 2.79 ) = ( 25.650 )
    //   (  2.5  1.1      6.0 )  ( 2.87 ) = ( 25.855 )

    // We use a format known as Compressed Sparse Column (CSC) to store the data. Further,
    // as the matrix is symmetric, we only need to store half the data. The CSC format
    // stores the matrix as a series of column vectors where only the non-zero entries are
    // specified, stored as the pair of (row index, value), although in separate arrays:

    let rowCount: Int32     = 4
    let columnCount: Int32  = 4
    var matrixValues        = [ 10.0, 1.0, 2.5, 12.0, -0.3, 1.1, 9.5, 6.0 ]
    var rowIndices: [Int32] = [    0,   1,   3,    1,    2,   3,   2,   3 ]
    var columnStarts        = [    0,              3,              6,   7 ]

    // vector X

    var xValues = [ 2.20, 2.85, 2.79, 2.87 ]

    // In this library, this raw information is all wrapper into a flexible data type
    // that allows for more complex use cases in other situations.

    rowIndices.withUnsafeMutableBufferPointer { rowIndicesPointer in
        columnStarts.withUnsafeMutableBufferPointer { columnStartsPointer in
            matrixValues.withUnsafeMutableBufferPointer { valuesPointer in
                xValues.withUnsafeMutableBufferPointer { xPointer in
                    let a = SparseMatrix_Double(
                        // Structure of the matrix, without any values
                        structure: SparseMatrixStructure(
                            rowCount:     rowCount,
                            columnCount:  columnCount,
                            columnStarts: columnStartsPointer.baseAddress!,
                            rowIndices:   rowIndicesPointer.baseAddress!,
                            // Matrix meta-data
                            attributes: SparseAttributes_t(
                                transpose: false,
                                triangle: SparseLowerTriangle,
                                kind: SparseSymmetric,
                                _reserved: 0,
                                _allocatedBySparse: false
                            ),
                            blockSize: 1),
                        // Numerical values of the matrix
                        data: valuesPointer.baseAddress!
                    )

                    let x = DenseVector_Double(count: columnCount, data: xPointer.baseAddress!)

                    let y = [Double](unsafeUninitializedCapacity: Int(rowCount)) { resultBuffer, count in
                        let y = DenseVector_Double(count: rowCount, data: resultBuffer.baseAddress!)
                        SparseMultiply(a, x, y)
                        count = Int(rowCount)
                    }

                    print(y) // [32.025, 38.72, 25.65, 25.855] – Correct
                }
            }
        }
    }
}

The problem is columnStarts. This array must include, after all the actual “column starts”, the number of elements in the array of matrix values (so it can calculate how many elements are in that last column). Note the extra value at the end of the columnStarts array:

func matrixProductExperiment() {
    …

    let rowCount: Int32     = 4
    let columnCount: Int32  = 4
    var matrixValues        = [ 10.0, 1.0, 2.5, 12.0, -0.3, 1.1, 9.5, 6.0 ]
    var rowIndices: [Int32] = [    0,   1,   3,    1,    2,   3,   2,   3 ]
    var columnStarts        = [    0,              3,              6,   7,   8 ]
    assert(columnStarts.last == matrixValues.count, "\(String(describing: columnStarts.last)) ≠ \(matrixValues.count)")

    …
}

If you use SparseConvertFromCoordinate, as shown in the SparseMultiply documentation, you supply the blockCount. But the SparseMatrixStructure initializer has no “block count” parameter, and it extracts this (somewhat unobviously, IMHO) from the last value in columnStarts.

See Creating Sparse Matrices, which says:

This array requires an additional, final entry that defines the final column's length.

To my eye, it was not entirely intuitive that the “starts” array must end with the “count” of the matrix values. And, annoyingly, the SparseMatrixStructure documentation neglected to point this out at the time (but now does). But, in retrospect, we can see why we must.

That is a few hours of my life that I will never get back. The solution may be obvious to many, but for those stumbling with random crashes with the sparse matrices in Accelerate, take a hard look at the columnStarts array. This value does not enjoy meaningful validation and/or error messages. One will receive intermittent/cryptic errors if the count is omitted.

It was a stupid mistake on my part, but tracking down the problem took a lot of work. Hopefully, this will save some future readers a few hours of their time.