CSC Matrix Operations -- Code I wrote before I realised that Eigen can be made to directly interface with CSC matrices via SparseMatrix and makeCompressed()

csc_matrix.cpp

/*
 *  libbioneuronqp -- Library solving for synaptic weights
 *  Copyright (C) 2020  Andreas Stöckel
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Affero General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Affero General Public License for more details.
 *
 *  You should have received a copy of the GNU Affero General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */

#include "csc_matrix.hpp"

// Do nothing here, just make sure the header compiles

// Unit test: In the build directory, compile using
// g++ -g -O3 ../bioneuronqp/csc_matrix.cpp -DCSC_MATRIX_TEST -I../extern -I/usr/include/eigen3 -otest_csc_matrix && ./test_csc_matrix

#ifdef CSC_MATRIX_TEST

#include <iostream>
#include <random>

int main()
{
	std::default_random_engine rng(48910);
	size_t n_test = 1000;
	for (size_t test = 0; test < n_test; test++) {
		// Create a randomly sized matrix
		std::uniform_int_distribution<c_int> dist_row_cols(1, 20);
		c_int m = dist_row_cols(rng), n = dist_row_cols(rng);
		CSCMatrix mat(m, n);
		Eigen::MatrixXd mat_ref = Eigen::MatrixXd::Zero(m, n);

		std::uniform_int_distribution<c_int> dist_rows(0, m - 1);
		std::uniform_int_distribution<c_int> dist_cols(0, n - 1);
		std::uniform_int_distribution<c_int> dist_x(-2, 2);
		for (size_t it = 0; it < 10000; it++) {
			// Draw random update locations
			c_int ui = dist_rows(rng), uj = dist_cols(rng);

			// Draw a value to update to
			c_float x = dist_x(rng);

			// Update the matrix
			mat.set(ui, uj, x);
			mat_ref(ui, uj) = x;

			/*std::cout << "Setting " << ui << ", " << uj << " to " << x <<
			std::endl; std::cout << "[" << std::endl; for (c_int i = 0; i < m;
			i++) { std::cout << " ["; for (c_int j = 0; j < n; j++) { std::cout
			<< mat.get(i, j) << ", ";
			    }
			    std::cout << "]," << std::endl;
			}
			std::cout << "]" << std::endl; */

			// Compare the CSC matrix instance and the reference matrix instance
			for (c_int i = 0; i < m; i++) {
				for (c_int j = 0; j < n; j++) {
					if (mat.get(i, j) != mat_ref(i, j)) {
						std::cout << std::endl
						          << "Error comparing values at " << i << ", "
						          << j << "; expected " << mat_ref(i, j)
						          << " but got " << mat.get(i, j) << std::endl;
						return 1;
					}
				}
			}
		}
		std::cout << "\r" << int(100.0 * double(test + 1) / double(n_test))
		          << "% done...";
	}
	std::cout << std::endl << "OK" << std::endl;
	return 0;
}

#endif /* CSC_MATRIX_TEST */

csc_matrix.hpp

/*
 *  libbioneuronqp -- Library solving for synaptic weights
 *  Copyright (C) 2020  Andreas Stöckel
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Affero General Public License as
 *  published by the Free Software Foundation, either version 3 of the
 *  License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Affero General Public License for more details.
 *
 *  You should have received a copy of the GNU Affero General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */

#ifndef BIONEURONQP_CSC_MATRIX_HPP
#define BIONEURONQP_CSC_MATRIX_HPP

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <vector>

#include <osqp/osqp.h>

#include <Eigen/Dense>

/**
 * This class represents a matrix in the Compressed Sparse Column (CSC) format
 * using zero-based indexing.
 *
 * See https://people.sc.fsu.edu/~jburkardt/data/cc/cc.html for a description of
 * the format.
 */
struct CSCMatrix {
private:
	csc m_csc;
	std::vector<c_int> m_col;
	std::vector<c_int> m_row;
	std::vector<c_float> m_x;

public:
	CSCMatrix(c_int rows, c_int cols)
	{
		assert(rows >= 0 && cols >= 0);

		// Fill the column vector
		m_col.push_back(0);
		for (c_int j = 0; j < cols; j++) {
			m_col.push_back(0);
		}

		// Set all properties of the csc matrix to initial values
		m_csc.nzmax = 0;
		m_csc.m = rows;
		m_csc.n = cols;
		m_csc.p = m_col.data();
		m_csc.i = nullptr;
		m_csc.x = nullptr;
		m_csc.nzmax = 0;
		m_csc.nz = -1;
	}

	CSCMatrix(const Eigen::MatrixXd &mat, bool upper_triangle = false)
	{
		c_int m = mat.rows(), n = mat.cols();

		m_col.push_back(0);
		for (c_int j = 0; j < n; j++) {
			for (c_int i = 0; i < (upper_triangle ? (j + 1) : m); i++) {
				if (mat(i, j) != 0.0) {
					m_row.push_back(i);
					m_x.push_back(mat(i, j));
				}
			}
			m_col.push_back(m_x.size());
		}

		m_csc.nzmax = m * n;
		m_csc.m = m;
		m_csc.n = n;
		m_csc.p = m_col.data();
		m_csc.i = m_row.data();
		m_csc.x = m_x.data();
		m_csc.nzmax = m_x.size();
		m_csc.nz = -1;
	}

	c_int rows() const { return m_csc.m; }

	c_int cols() const { return m_csc.n; }

	c_float get(c_int i, c_int j) const noexcept
	{
		// Make sure the given indices are not out of bounds
		assert(!(i < 0 || j < 0 || i >= rows() || j >= cols()));

		// Fetch the column pointers for the current and the next column
		c_int c0 = m_col[j], c1 = m_col[j + 1];
		c_int const *p0 = &m_row[c0], *p1 = &m_row[c1];

		// Search the requested row
		c_int const *p = std::lower_bound(p0, p1, i);
		if (p == p1 || *p != i) {
			return 0.0;  // Entry not found, return zero
		}

		// Compute the index the entry is pointing at
		ptrdiff_t idx = p - m_row.data();

		// Return the data associated with this row
		return m_x[idx];
	}

	void set(c_int i, c_int j, double x)
	{
		// Make sure the given indices are not out of bounds
		assert(!(i < 0 || j < 0 || i >= rows() || j >= cols()));

		// Fetch the column pointers for the current and the next column
		c_int c0 = m_col[j], c1 = m_col[j + 1];
		c_int const *p0 = &m_row[c0], *p1 = &m_row[c1];

		// Search the requested row
		c_int const *p = std::lower_bound(p0, p1, i);

		// Compute the index the entry is pointing at
		ptrdiff_t idx = p - m_row.data();

		if (p != p1 && *p == i) {
			// First case, the entry already exists
			if (x == 0.0) {
				// If we're setting the entry to 0, erase the corresponding
				// entry and decrease the column pointer for the next column.
				m_row.erase(m_row.begin() + idx);
				m_x.erase(m_x.begin() + idx);
				for (c_int k = j + 1; k <= cols(); k++) {
					m_col[k]--;
				}
				m_csc.nzmax--;
			} else {
				// Otherwise just overwrite the corresponding entry in the x
				// matrix
				m_x[idx] = x;
				return; // No changes to any pointers
			}
		} else {
			// Second case, the entry does not exist
			if (x == 0.0) {
				return; // No-op, do not create entries just containing zeros
			}

			// Otherwise we need to insert a new entry into m_row and m_x
			m_row.insert(m_row.begin() + idx, i);
			m_x.insert(m_x.begin() + idx, x);
			for (c_int k = j + 1; k <= cols(); k++) {
				m_col[k]++;
			}
			m_csc.nzmax++;
		}

		// Update the data pointers since we might have resized the underlying
		// matrix
		m_csc.i = m_row.data();
		m_csc.x = m_x.data();
	}

	operator csc *() { return &m_csc; }
};


#endif /* BIONEURONQP_CSC_MATRIX_HPP */