With the release of <strong>ArrayFire v3.4</strong>, we have added a number of features that will help our vast range of users. Our users have been asking us for a way to reduce the memory footprint of their data on the GPU. And now we have answered our users' requests.
Introducing Sparse Matrices in ArrayFire!

<h4>Sparse Data Structure</h4>
A sparse data structure is one where all the non-zero elements are not stored. Sparse matrices are useful when the number of zero-values elements are &gt;&gt; number of non-zero elements (sparsity of the matrix is high).
The sparse data structure in general is stored as 3 arrays.
<ul>
<li>A data or values array containing all the non-zero elements</li>
<li>A vector for row indices (based on storage format).</li>
<li>A vector for column indices (based on storage format).</li>
</ul>

There are many ways to store sparse matrices, the most prominent of which are
<ul>
<li><a href="https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_.28CSR.2C_CRS_or_Yale_format.29">Compressed Sparse Row</a></li>
<li><a href="https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_.28CSC_or_CCS.29">Compressed Storage Column</a></li>
<li><a href="https://en.wikipedia.org/wiki/Sparse_matrix#Coordinate_list_.28COO.29">Coordinate Format</a>.</li>
</ul>

ArrayFire supports CSR and COO format, with CSC support coming in future releases.

A more detailed description of Sparse Matrix can be found on <a href="https://en.wikipedia.org/wiki/Sparse_matrix">Wikipedia</a>.

<h4>Sparse in ArrayFire</h4>
In ArrayFire, Sparse arrays use the same user facing object - af::array. This makes working with sparse arrays all the more seamless.
We have added a function <a href="http://arrayfire.org/docs/group__method__mat.htm#ga5945f95549bae25f75694a69d98ca584">af::array::issparse()</a> that allows users to check if an array is sparse or not.
Although externally, there is no difference between a sparse and dense array, internally, the sparse array is stored as 3 independent arrays along with other meta data.
These arrays, as well as other properties can be fetched using the <a href="http://arrayfire.org/docs/group__sparse__func.htm">sparse functions</a>.

<h4>Sparse Functionality</h4>
As of v3.4, the following features are supported for sparse arrays:
<ul>
<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm">Creation for sparse arrays.</a>
<ul>
<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm#gae4621e58bfaf6d80668309824fde86c0">From 3 (data, row, column) sparse host or device memory vectors.</a></li>
<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm#gae9a843c7e93b05e6b090ca32b271ebe9">From 3 (data, row, column) sparse af::arrays.</a></li>
<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm#gac0d367d8c9daf27b36a6981bbfcdb970">From dense af::array.</a></li>
<li>To create a sparse array from host or device dense array, first create a dense af::array, and then pass it to create sparse af::array from dense af::array.</li>
</ul>
</li>
<li><a href="http://arrayfire.org/docs/group__sparse__func__convert__to.htm">Storage format conversions.</a>
<ul>
<li><a href="http://arrayfire.org/docs/group__sparse__func__create.htm#gac0d367d8c9daf27b36a6981bbfcdb970">Dense to CSR/COO.</a></li>
<li><a href="http://arrayfire.org/docs/group__sparse__func__dense.htm">CSR/COO to dense.</a></li>
</ul>
</li>
<li><a href="http://arrayfire.org/docs/group__blas__func__matmul.htm#ga63306b6ed967bd1055086db862fe885b">Sparse-Dense Matrix Multiplication (SPMM/SPMV)</a>
<ul>
<li>Sparse Matrix - Dense Matrix</li>
<li>Sparse Matrix - Dense Vector</li>
<li>Options for Sparse Matrix to be transposed.</li>
</ul>
</li>
</ul>

The Sparse-Dense Matrix Multiplcation is a implemented as a part of <a href="http://arrayfire.org/docs/group__blas__func__matmul.htm#ga63306b6ed967bd1055086db862fe885b">af::matmul()</a>. This is part of making the usage of sparse arrays seamless with regular arrays.

<h4>Example</h4>
<pre>
/*******************************************************
 * Copyright (c) 2014, ArrayFire
 * All rights reserved.
 *
 * This file is distributed under 3-clause BSD license.
 * The complete license agreement can be obtained at:
 * http://sandbox.arrayfire.com/licenses/BSD-3-Clause
 ********************************************************/
#include <arrayfire.h>
#include <iostream>

using namespace af;

static size_t dimension = 16 * 1024;
static const int maxIter = 10;
static const int sparsityFactor = 7;

static array A;
static array As;  // Sparse A
static array x0;
static array b;

void setupInputs()
{
    // Generate a random input: A
    array T = randu(dimension, dimension, f32);
    // Create 0s in input.
    // Anything that is no divisible by sparsityFactor will become 0.
    A = floor(T * 1000);
    A = A * ((A % sparsityFactor) == 0) / 1000;
    // Make it positive definite
    A = transpose(A) + A + A.dims(0)*identity(A.dims(0), A.dims(0), f32);

    // Make A sparse as As
    As = sparse(A);

    // Generate x0: Random guess
    x0 = randu(A.dims(0), f32);

    //Generate b
    b = matmul(A, x0);

    // Run evals on inputs
    A .eval();
    sparseGetValues(As).eval();
    sparseGetRowIdx(As).eval();
    sparseGetColIdx(As).eval();
    x0.eval();
    b .eval();
    af::sync();
}

void sparseConjugateGradient(void)
{
    array x = constant(0, b.dims(), f32);
    array r = b - matmul(As, x);
    array p = r;

    for (int i = 0; i < maxIter; ++i) {
        array Ap = matmul(As, p);
        array alpha_num = dot(r, r);
        array alpha_den = dot(p, Ap);
        array alpha = alpha_num/alpha_den;
        r = r - tile(alpha, Ap.dims())*Ap;
        x = x + tile(alpha, Ap.dims())*p;
        array beta_num = dot(r, r);
        array beta = beta_num/alpha_num;
        p = r + tile(beta, p.dims())*p;
    }
}

void denseConjugateGradient(void)
{
    array x = constant(0, b.dims(), f32);
    array r = b - matmul(A, x);
    array p = r;

    for (int i = 0; i < maxIter; ++i) {
        array Ap = matmul(A, p);
        array alpha_num = dot(r, r);
        array alpha_den = dot(p, Ap);
        array alpha = alpha_num/alpha_den;
        r = r - tile(alpha, Ap.dims())*Ap;
        x = x + tile(alpha, Ap.dims())*p;
        array beta_num = dot(r, r);
        array beta = beta_num/alpha_num;
        p = r + tile(beta, p.dims())*p;
    }
}

void checkConjugateGradient(const af::array in)
{
    array x = constant(0, b.dims(), f32);
    array r = b - matmul(in, x);
    array p = r;

    for (int i = 0; i < maxIter; ++i) {
        array Ap = matmul(in, p);
        array alpha_num = dot(r, r);
        array alpha_den = dot(p, Ap);
        array alpha = alpha_num/alpha_den;
        r = r - tile(alpha, Ap.dims())*Ap;
        x = x + tile(alpha, Ap.dims())*p;
        array beta_num = dot(r, r);
        array beta = beta_num/alpha_num;
        p = r + tile(beta, p.dims())*p;
    }
    array res = x0 - x;

    std::cout<<"Final difference in solutions:\n";
    af_print(dot(res, res));
}

int main(int argc, char *argv[])
{
    af::info();
    setupInputs();

    std::cout << "Verifying Dense Conjugate Gradient:" << std::endl;
    checkConjugateGradient(A);

    std::cout << "Verifying Sparse Conjugate Gradient:" << std::endl;
    checkConjugateGradient(As);

    af::sync();

    std::cout << "Dense Conjugate Gradient Time: "
              << timeit(denseConjugateGradient) * 1000 << "ms" << std::endl;
    af::sync();

    std::cout << "Sparse Conjugate Gradient Time: "
              << timeit(sparseConjugateGradient) * 1000 << "ms" << std::endl;
    af::sync();

    return 0;
}
</pre>
<h4>Performance and Memory Footprint</h4>

As with all ArrayFire features, Sparse support is available across all backends and all architectures.