File: SparseGenMatProd.cpp

package info (click to toggle)
spectra 1.2.0-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 1,788 kB
  • sloc: cpp: 23,044; ansic: 175; fortran: 131; makefile: 90
file content (53 lines) | stat: -rw-r--r-- 1,606 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Spectra/MatOp/SparseGenMatProd.h>

#include "catch.hpp"

using namespace Spectra;

using Eigen::Index;
template <typename TestType>
using Matrix = Eigen::Matrix<TestType, Eigen::Dynamic, Eigen::Dynamic>;

template <typename T>
Eigen::SparseMatrix<T> generate_random_sparse(Index rows, Index cols)
{
    std::default_random_engine gen;
    std::uniform_real_distribution<T> dist(0.0, 1.0);

    std::vector<Eigen::Triplet<T>> tripletVector;
    for (int i = 0; i < rows; ++i)
        for (int j = 0; j < cols; ++j)
        {
            auto v_ij = dist(gen);
            if (v_ij < 0.5)
            {
                // if larger than treshold, insert it
                tripletVector.push_back(Eigen::Triplet<T>(i, j, v_ij));
            }
        }
    Eigen::SparseMatrix<T> mat(rows, cols);
    // create the matrix
    mat.setFromTriplets(tripletVector.begin(), tripletVector.end());

    return mat;
}

TEMPLATE_TEST_CASE("matrix operations [100x100]", "[SparseGenMatProd]", float, double)
{
    std::srand(123);
    constexpr Index n = 100;

    Eigen::SparseMatrix<TestType> mat1 = generate_random_sparse<TestType>(n, n);
    Matrix<TestType> mat2 = Matrix<TestType>::Random(n, n);

    SparseGenMatProd<TestType> sparse1(mat1);
    Matrix<TestType> xs = sparse1 * mat2;
    Matrix<TestType> ys = mat1 * mat2;

    INFO("The matrix-matrix product must be the same as in eigen.")
    REQUIRE(xs.isApprox(ys));
    INFO("The accesor operator must produce the same element as in eigen")
    REQUIRE(mat1.coeff(45, 22) == sparse1(45, 22));
}