RcppEigen/inst/tinytest/cpp/sparse.cpp at master · RcppDeepState/RcppEigen · GitHub

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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
Eigen::SparseMatrix<double> wrapSparseDouble() {
    Eigen::SparseMatrix<double>  mm(9,3);
    mm.reserve(9);
    for (int j = 0; j < 3; ++j) {
        mm.startVec(j);
        for (int i = 3 * j; i < 3 * (j + 1); ++i)
            mm.insertBack(i, j) = 1.;
    }
    mm.finalize();
    return mm;
}

// [[Rcpp::export]]
Eigen::SparseMatrix<double, Eigen::ColMajor> wrapSparseDoubleColumnMajor() {
    Eigen::SparseMatrix<double, Eigen::ColMajor>  mm(9,3);
    mm.reserve(9);
    for (int j = 0; j < 3; ++j) {
        mm.startVec(j);
        for (int i = 3 * j; i < 3 * (j + 1); ++i)
            mm.insertBack(i, j) = 1.;
    }
    mm.finalize();
    return mm;
}

// [[Rcpp::export]]
Eigen::SparseMatrix<double, Eigen::RowMajor> wrapSparseDoubleRowMajor() {
    Eigen::SparseMatrix<double, Eigen::RowMajor>  mm(9,3);
    mm.reserve(9);
    for (int irow = 0; irow < 9; ++irow) {
        mm.startVec(irow);
        mm.insertBack(irow, irow / 3) = static_cast<double>( 9 - irow );
    }
    mm.finalize();
    return mm;
}

// [[Rcpp::export]]
Eigen::SparseMatrix<double, Eigen::ColMajor> asSparseDoubleColumnMajor(Eigen::SparseMatrix<double, Eigen::ColMajor> mm) {
    return mm;
}

// [[Rcpp::export]]
double asMappedSparseDoubleColMajor(Eigen::Map<Eigen::SparseMatrix<double, Eigen::ColMajor> > mm) {
    double s = mm.sum();          // access instantiated sparse matrix
    return s;
}

// [[Rcpp::export]]
double asMappedSparseDeprecatedDoubleColMajor(Eigen::MappedSparseMatrix<double, Eigen::ColMajor> mm) {
    // Deprecated
    double s = mm.sum();          // access instantiated sparse matrix
    return s;
}

// [[Rcpp::export]]
double asSparseDoubleRowMajor(Eigen::SparseMatrix<double, Eigen::RowMajor> mm) {
    double s = mm.sum();          // access instantiated sparse matrix
    return s;
}

// [[Rcpp::export]]
double asMappedSparseDoubleRowMajor(Eigen::Map<Eigen::SparseMatrix<double, Eigen::RowMajor> > mm) {
    double s = mm.sum();          // access instantiated sparse matrix
    return s;
}

// [[Rcpp::export]]
double asMappedSparseDeprecatedDoubleRowMajor(Eigen::MappedSparseMatrix<double, Eigen::RowMajor> mm) {
    double s = mm.sum();          // access instantiated sparse matrix
    return s;
}

// [[Rcpp::export]]
Rcpp::List sparseCholesky(Rcpp::List input) {
    using Eigen::VectorXd;
    using Eigen::MatrixXd;
    using Eigen::Lower;
    using Eigen::Map;
    using Eigen::SparseMatrix;
    using Eigen::SimplicialLDLT;
    using Eigen::Success;

    const Map<SparseMatrix<double> > m1 = input[0];
    const Map<VectorXd>              v1 = input[1];
    SparseMatrix<double>             m2(m1.cols(), m1.cols());
    m2.selfadjointView<Lower>().rankUpdate(m1.adjoint());

    SimplicialLDLT<SparseMatrix<double> > ff(m2);
    VectorXd                        res = ff.solve(m1.adjoint() * v1);

    return Rcpp::List::create(Rcpp::Named("res")   = res,
                              Rcpp::Named("rows")  = double(ff.rows()),
                              Rcpp::Named("cols")  = double(ff.cols()));

}