TestEigenCSR.cpp
/**
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#include <gtest/gtest.h>
#include <Eigen/Sparse>
#include <numeric>
#include <vector>
/**
* This test case checks if the internal Eigen::SparseMatrix compressed storage
* format is a conventional CSR matrix. Currently this is the case, but it is
* not guaranteed for all time.
*
* Cf. section "Sparse matrix format" on page
* http://eigen.tuxfamily.org/dox/group__TutorialSparse.html
*/
TEST(MathLibEigen, Eigen2CSR)
{
const int nrows = 16;
const int ncols = nrows;
Eigen::SparseMatrix<double, Eigen::RowMajor> mat(nrows, ncols);
// set up sparsity pattern
std::vector<int> pat(nrows);
for (std::size_t i = 0; i < nrows; ++i)
{
if (i == 0 || i == nrows - 1)
{
pat[i] = 2;
}
else
{
pat[i] = 3;
}
}
// CSR representation of the matrix
std::vector<double> values;
std::vector<int> ia(nrows + 1); // row offsets
std::vector<int> ja; // column indices
std::partial_sum(pat.begin(), pat.end(), ia.begin() + 1);
const int nnz = ia.back();
values.reserve(nnz);
ja.reserve(nnz);
mat.reserve(pat);
// init matrix, build CSR matrix in parallel
for (int row = 0; row < nrows; ++row)
{
for (int col = -1; col <= 1; ++col)
{
int cidx = row + col;
if (cidx < 0 || cidx >= ncols)
{
continue;
}
const double val = (col == 0) ? 2.0 : -1.0;
values.push_back(val);
mat.coeffRef(row, cidx) = val;
ja.push_back(cidx);
}
}
// change matrix
for (int row = 0; row < nrows; ++row)
{
for (int col = -1; col <= 1; ++col)
{
int cidx = row + col;
if (cidx < 0 || cidx >= ncols)
{
continue;
}
mat.coeffRef(row, cidx) = 4.0 * mat.coeff(row, cidx);
}
}
// adapt entries of CSR matrix
for (auto& v : values)
{
v *= 4.0;
}
mat.makeCompressed();
ASSERT_EQ(nnz, mat.nonZeros());
ASSERT_EQ(nnz, values.size());
ASSERT_EQ(nnz, ja.size());
int* ptr = mat.outerIndexPtr();
int* col = mat.innerIndexPtr();
double* data = mat.valuePtr();
for (int r = 0; r < nrows; ++r)
{
EXPECT_EQ(ia[r], ptr[r]);
}
for (int i = 0; i < nnz; ++i)
{
EXPECT_EQ(values[i], data[i]);
EXPECT_EQ(ja[i], col[i]);
}
}