Basic Datatypes#

MOCCA currently offers the following data types:

Name

Format

Size

Dimensions

Container

Vector

Dense

Dynamic

1

Array

Matrix

Dense

Dynamic

2

Array

Static Matrix

Dense

Static

1 or 2

std::array

Sparse Vector

Sparse

Dynamic

1

Array

CSR Matrix

Sparse

Dynamic

2

Array

All numeric types (32-bit and 64-bit integers, float and double) are supported. For now, MOCCA only has official support for real numbers.

Capacity:

  • The number of rows and columns for any vector or matrix can be retrieved by calling the rows() and cols() methods, respectively.

  • The size() routine returns either the number of entries when working with dense objects or the number of nonzero elements when working with sparse objects.

  • The size of a dynamic-sized object can be modified during the runtime by calling the resize() method. This routine is destructive (i.e., all previous allocation is removed before resizing) for all, except for the Array.

  • You can use reshape() to change the dimensions of the object without altering its content. Note that resizing the container may lead to data reallocation and copying.

Assignment:

  • All data types support the assignment operator = to either copy or move data between two objects. Assigning a vector to a matrix will trigger a Broadcast.

  • Trying to assign a dense object to a sparse one will result in an assertion failure.

  • When assigning an expression to a sparse object, it must have the same sparse format (i.e., CSR, BSR, etc.).

Array#

The Array object corresponds to a plain contiguous array that supports 64-bit indexing and is aligned to the SIMD width (e.g., 256-bit for AVX2), while also providing routines for easy data manipulation as well as a random access iterator for full compatibility with the STL and range-based loops. The Array can only hold POD`s (plain old data types, such as `float, int, std::array, etc.) as it will not call the object constructor when allocating the elements.

Similar to a std::vector, the Array have an allocation capacity and a size. The capacity stores the size of the current memory allocation, while the size stores the number of elements in the container. In this way, an Array can be resized at will without altering the memory allocation as long as the new size is not greater than the current capacity(). Otherwise, MOCCA will allocate a new block in memory and move the content of the container to this new block. To avoid constant reallocations when modifying the size of the container, its allocation capacity is always a multiple of the MOCCA_ALLOC_BLOCK (default: 1024 elements).

Static Matrix#

The Static Matrix represents a fixed-sized matrix. If one of the dimensions is set to 1, this object will behave like a vector. In this situation, you can use the aliases StaticRowVector and StaticVector.

Initialization#

A Static Matrix can be constructed in several ways:

using namespace mocca;

StaticMatrix<int, 5, 4> mat1;                                       // Creates an unitialized 5x4 integer matrix
StaticMatrix<int, 5, 5> mat2(7);                                    // Creates a 5x5 integer matrix filled with the value `7`
StaticMatrix<float, 5, 5> mat3 = mat2                               // Creates a copy of the mat3
StaticMatrix<float, 2, 3> mat4 = {1, 2, 3, 4, 5, 6};                // Creates a 2x3 matrix and fills it with the data from {...}
StaticMatrix<float, 5, 5> mat5 = mat3 + mat2;                       // Evalutes the sum between mat3 and mat2 and stores the results in mat5

Element Access#

The elements of the Static Matrix can be accessed with the subscript operator:

mocca::StaticMatrix<int, 4, 2> m;
m(2, 1) = 1;             // Access element at position (2, 1)
m[1 + 1 * m.cols()] = 3;    // Access element at position (1, 1) using the C-style notation

Vector#

The Vector class represent either a column or a row vector. All Vector are column vectors by default, unless the RowVector alias is used.

Initialization#

There are several ways to construct a Vector:

using namespace mocca;

Vector<int> vec1;                                       // Creates an empty column vector
Vector<double> vec2(10);                                // Creates a column vector with 10 doubles initialized at 0
RowVector<double> vec3(10, 3.4);                        // Creates a row vector with 10 doubles initialized at 3.4
Vector<double> vec4 = vec2;                             // Creates a copy of the vec2
Vector<int> vec5 = {1, 2, 3, 4};                        // Creates a column vector with the data specified by {...}
Vector<int> vec6 = 2 * vec5;                            // Evaluates 2 * vec5 and stores the results in vec6
Vector<int> vec7 = io::read_vec_txt<int>("f.vec");      // Creates a column vector from the file f.vec

Element Access#

The elements in the Vector can be accessed with the subscript operator:

v[1] = 0.1;     // Access element at position (1)
v(2) = 0.2;    // Access element at position (2)

Matrix#

The Matrix stores the data as a contiguous 2D array following the Row Major format. For example, the matrix

\[\begin{split}M = \begin{bmatrix} 5 & 7 & 8\\ 1 & 3 & 2\\ 7 & 9 & 1 \end{bmatrix} \quad\end{split}\]

is stored in the memory as

\[\begin{bmatrix} 5 & 7 & 8 & 1 & 3 & 2 & 7 & 9 & 1 \end{bmatrix}\]

Initialization#

A Matrix can be constructed in several ways:

using namespace mocca;

Matrix<int> mat1;                                    // Creates a empty integer matrix
Matrix<int> mat2(5, 5);                              // Creates a 5x5 integer matrix filled with zeroes
Matrix<float> mat3(3, 3, 2.2);                       // Creates a 3x3 floating-point matrix filled with 2.2
Matrix<float> mat4 = mat3;                           // Creates a copy of the mat3
Matrix<float> mat5 = {{1, 2, 3},                     // Creates a 2x3 matrix and fills it with the data from {...}
                      {4, 5, 6}}
Matrix<float> mat6 = mat3 + mat4;                    // Evalutes the sum between mat3 and mat4 and stores the results in mat6
Matrix<float> mat7(3, 3, mat3.begin(), mat3.end());  // Creates a 3x3 matrix and copies the content from the mat3 using iterators
Matrix<int> mat8 = io::read_dense_mtx<int>("f.mtx"); // Creates a matrix from the file f.mtx (MatrixMarket)

Element Access#

The elements of the Matrix can be accessed with the subscript operator:

mocca::Matrix<int> m(5, 5);
m(2, 1) = 1;             // Access element at position (2, 1)
m[3 + 4 * m.cols()] = 3;    // Access element at position (4, 3) using the C-style notation

CSR Matrix#

The CSR Matrix implements a variant of the widely-used Compressed Sparse Row (CSR) scheme. It consists of four buffers:

  • nz_value - Stores the values of each nonzero element

  • nz_column - Stores the column of each nonzero element

  • row_offset_ptr - Stores the index to the first element in the two previous for each row

  • nz_per_row - Stores the number of nonzero elements per row (optional)

Let us illustrate this scheme with an example:

\[\begin{split}M = \begin{bmatrix} 11 & 0 & 13 & 0 \\ 0 & 0 & 0 & 24 \\ 0 & 32 & 33 & 0 \\ 41 & 0 & 0 & 0 \end{bmatrix} \\\end{split}\]

Internally, the sparse matrix M is stored as follows:

\[\begin{split}&nz\_value = \begin{bmatrix} 11 & 13 & \_ & 24 & \_ & 32 & 33 & 41 & \_ \end{bmatrix} \\ &nz\_column = \begin{bmatrix} 0 & 2 & \_ & 3 & \_ & 1 & 2 & 0 & \_ \end{bmatrix} \\ &row\_offset\_ptr = \begin{bmatrix} 0 & 3 & 5 & 7 & 9 \end{bmatrix} \\ &nz\_per\_row = \begin{bmatrix} 2 & 1 & 2 & 1 \end{bmatrix}\end{split}\]

Note that all non-zeros entries within a single row are always sorted by their columns. The _ indicates an available space in the corresponding row for quickly inserting new elements. In this case, assuming no reallocation, inserting a random element in the i-th row only takes O(nz_per_row[i]) operations, while inserting elements sorted by columns only requires increasing nz_per_row[i] by one which is an O(1) operation.

If there are no empty spaces available, CSR Matrix have the same layout as the popular Compressed Sparse Row (CSR) format. This case is referred to as compressed mode and can be enabled by calling the compress() routine. In this mode, the nz_per_row buffer is redundant, and therefore, will be deallocated. Note that all MOCCA operations will always produce a compressed object and most external libraries require the matrix to be on the standard CSR layout.

Initialization#

A CSR Matrix can be created as follows:

using namespace mocca;

CSRMatrix<float> mat1;              //  Creates a empty CSRMatrix
CSRMatrix<int> mat2(5, 5);          //  Creates a 5x5 CSRMatrix filled with 0s
CSRMatrix<int> mat3(5, 4, 5);       //  Creates a 5x4 CSRMatrix and reserves space for 5 non-zero entries
CSRMatrix<int> mat4 = mat3;         //  Creates a copy of mat3
CSRMatrix<float> mat5 = mat3 + mat2 //  Evaluates the sum and assign the results into mat5

Due to the sparse storage scheme, inserting new non-zero entries to the CSR Matrix requires special care. For instance, inserting elements in a random order has a computational cost of \(O(nz)\) where \(nz\) is the total number of non-zeros before each insertion. A better and simpler way is to first build a list of Sparse Triplet objects (i.e., a simple struct containing the row, column and value of the non-zero entry) and then convert this list to CSR Matrix using either fill_sorted(), fill_sorted_by_col() or fill_unsorted() methods. For example,

using namespace mocca;

CSRMatrix<float> mat(nrows, ncols);
std::vector<SparseTriplet<float>> list;
list.reserve(nz);

for (...) {
    // ...
    list.push_back({.row = i, .col = j, .val = v});
}

mat.fill_unsorted(list.begin(), list.end());

Note

fill_sorted is the fastest and consumes less memory, but requires that list is fully sorted (i.e., the entries are sorted by row and then by column in ascending order), while fill_sorted_by_col handles the case where the entries are unsorted but the entries for the same row have ascending columns. For example,

(row, col, val)
(1, 2, 40)
(2, 6, 50)
(1, 5, 30)
(1, 10, 10)
...

Last, but not least, fill_unsorted handles unsorted ranges with duplicates at the cost of slower performance and possibly greater memory consumption.

Another way is to preallocate the memory space for the non-zeros entries and then perform a sorted insertion. For example,

using namespace mocca;

CSRMatrix<float> mat(nrows, ncols);
mat.reserve(n);                                // Preallocating for nz entries

for (int i = 0; i < nrows; ++i) {
    for (int j = 0; j < ncols; ++j) {
        // ...
        if (v != 0) mat.insert(i, j, v);        // Or mat.at_ref(i, j) = v;
    }
}

mat.compress();                                 // Optional

The reserve() method allocates space in memory for n non-zero entries in the CSR Matrix, which are typically known in advance. Since the entries are inserted in order, the cost of the operation is only \(O(1)\).

Element Access#

Some algorithms only require iterating over the non-zero entries with no regard to their position in the CSR Matrix. In this case, the elements of the matrix can be accessed using the C-style subscript operator:

using namespace mocca;

CSRMatrix<float> mat(nrows, ncols);
// ...

mat[0] = v;
mat[20] = u;

The information of a single non-zero entry at the i-th position in the matrix (counting only the non-zeros) can be retrieved using the nonzero_at method:

SparseTriplet t = mat.nonzero_at(i);

Note that this method will perform a binary search for determining the row of the non-zero element (a \(O(nrows)\) operation). If the row is already known, a better approach is to first create a row view and then retrieve the column and value of the non-zero from the view, avoiding the binary search:

auto csr_row = mat.row(r);
auto [col, val] = csr_row.nonzero_at(i);

An element at the (i, j) position in the matrix can be accessed in several ways:

std::cout << mat(i, j) << std::endl;
std::cout << mat.at(i, j) << std::endl;
mat.at_ref(i, j) = v;

Note that these operations require a binary search over all non-zeros entries in the i-th row to determine if the element in the (i, j) exist or not in the CSR Matrix. Only the at_ref routine will create a new element if it doesn’t. Assign values to mat(i, j) or mat.at(i, j) have no effect.

Sparse Vector#

The Sparse Vector class stores only the non-zero elements of a column or row vector. All Sparse Vector are column vectors by default, unless the SparseRowVector alias is used.

Initialization#

A Sparse Vector can be created as follows:

using namespace mocca;

SparseVector<float> vec1;                   //  Creates a empty sparse column vector
SparseRowVector<int> vec2(50);              //  Creates a sparse row vector with 50 integers, initialized with 0s
SparseVector<int> vec3(50, 10);             //  Creates a sparse column vector with 50 integers and reserves space for 10 non-zero entries
SparseRowVector<int> vec4 = vec3;           //  Creates a copy of the vec3
SparseVector<float> vec5 = vec3 + vec2      //  Evaluates the sum and assign the results into vec5

The best way to fill_sorted_by_col a Sparse Vector is to reserve space for all non-zero entries and then insert them in order:

using namespace mocca;

SparseVector<int> vec(n, nz);

for (int i = 0; i < n; ++i) {
    // ...
    if (v_i != 0) vec.insert(i, v_i);
}

Element Access#

The elements of the Sparse Vector can be accessed using two types of indexing. Counting only the nonzeros, the element at the i-th position can be accessed as

using namespace mocca;

SparseVector<double> vec(n, nz);

// ...

vec[i] = 10;
auto [col, val] = vec.nonzero_at(i);

The second type of indexing considers all the elements in the Sparse Vector, including zero entries. In this case, the element at the i-th position can be accessed as

std::cout << vec(i) << std::endl;
std::cout << vec.at(i) << std::endl;
vec.at_ref(i) = v;

Note that these operations require a binary search over all non-zeros entries to determine if the element in the (i) position exists or not in the Sparse Vector. Only the at_ref routine creates a new element if it doesn’t. Assign values to vec(i) or vec.at(i) have no effect.

Special Matrices and Vectors#

MOCCA also provided an easy way to initialize dense matrices and vectors that have special properties, such as the Identity Matrix \(\mathbf{I}\):

using namespace mocca;

Matrix<int> mat_ones = ones<int>(5, 4);             // Creates a 5x4 matrix filled with ones
Matrix<float> mat_zeros = zeros<float>(100, 100);   // Creates 100x100 matrix filled with zeros
Matrix<int> mat_const = const_matrix(20, 20, 5);    // Creates 20x20 matrix filled with the value `5`
Matrix<double> mat_eye = eye<double>(10, 10);       // Creates a 10x10 identity matrix I

Vector<float> seq1 = arange(-10.0f, 10.0f, 0.5f);   // Creates an evenly-spaced sequence in [-10, 10[ with a step = 0.5
Vector<int> seq2 = linspace(0, 100, 25);            // Creates an evenly-spaced sequence in [0, 100] with 25 elements (+ the end of the interval)

These special matrices can also be used as part of an expression and do not allocate additional memory (since they are generated from a mathematical formula). For example, the code below sums the 5 to the diagonal entries of mat1 and then assigns the results to the matrix mat2:

mocca::Matrix<int> mat2 = mat1 + 5 * eye<int>(5, 5);

You can also create a custom initializator:

using namespace mocca;

Matrix<float> A = init_matrix<float>(5, 5, [](index_t row, index_t col) {
    return row == col ? 5 : 1;
});

std::cout << "A = " << A << std::endl;

Vector<int> v = init_vector<int>(4, [](index_t i) {
    return 10 * i * M_PI;
});

std::cout << "v = " << v << std::endl;

Output:

A =
[[   5     1     1     1     1   ]
 [   1     5     1     1     1   ]
 [   1     1     5     1     1   ]
 [   1     1     1     5     1   ]
 [   1     1     1     1     5   ]]

v =
[   0
    31
    62
    94   ]