/**
* \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 <vector>
#include <numeric>
#include <Eigen/Sparse>
/**
* 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]);
}
}