Two other small points. I am assuming that the set of indices in the
second argument is sorted. Also, it would be best to declare the function
as
bool foo(const MSpMat X, const Eigen::Map<Eigen::VectorXi> idx) as in the
enclosed
On Tue, Mar 25, 2014 at 10:45 AM, Douglas Bates <[email protected]> wrote:
> The enclosed works on this example. The logic is to check each column in
> the index set for all the elements of the index set, except the one on the
> diagonal, being in the set of row indices. I have added some diagnostic
> output to demonstrate the flow.
>
>
>
>
> On Tue, Mar 25, 2014 at 9:48 AM, Douglas Bates <[email protected]>wrote:
>
>> On Tue, Mar 25, 2014 at 6:42 AM, Søren Højsgaard <[email protected]>wrote:
>>
>>> Dear all,
>>> I have a large sparse adjacency matrix X for an undirected graph. For a
>>> subset 'idx' of the vertices I want to find out if the subgraph defined by
>>> this subset is complete (i.e. has edges between all variables). So in R,
>>> one would check if X[idx,idx] has only ones outside the diagonal, but this
>>> is slow and I want to do it with Rcpp. Here is one attempt where I simply
>>> add up the elements above (or below) the diagonal, and to access the
>>> elements of X I use coeff which is relatively slow (because X is a sparse
>>> matrix).
>>>
>>> #include <RcppEigen.h>
>>> typedef Eigen::MappedSparseMatrix<double> MSpMat;
>>> typedef Eigen::SparseVector<double> SpVec;
>>> typedef SpVec::InnerIterator InIter;
>>>
>>> //[[Rcpp::export]]
>>> double foo (MSpMat X, Eigen::VectorXi idx){
>>> SpVec sidx = idx.sparseView();
>>> double out=0;
>>> int i2, j2;
>>> for (InIter i_(sidx); i_; ++i_){
>>> i2 = i_.index();
>>> for (InIter j_(sidx); j_; ++j_){
>>> j2 = j_.index();
>>> if (i2>j2)
>>> out += X.coeff( i2, j2);
>>> }
>>> }
>>> return out;
>>> }
>>>
>>> /*** R
>>> library(Matrix)
>>> M1 <- new("dgCMatrix"
>>> , i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L,
>>> 5L, 3L, 4L)
>>> , p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
>>> , Dim = c(6L, 6L)
>>> , Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c",
>>> "d", "e", "f"))
>>> , x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
>>> , factors = list()
>>> )
>>> M1
>>>
>>> foo(M1, c(2,3,4))
>>>
>>> */
>>> Can anyone see a better way of doing this?
>>>
>>> I was thinking about whether the sparse matrix iterators could be a
>>> solution, but I can't a hold on it.
>>
>>
>> Sparse matrix inner and outer iterators are definitely the way to go.
>> I'll write a version using them.
>>
>>
>>> I was also thinking about using sparse matrices in RcppArmadillo, but I
>>> guess the problem (speed) is the same (haven't tried!).
>>>
>>
>>
>>> Any advice will be greatly appreciated.
>>> Cheers
>>> Søren
>>>
>>>
>>>
>>> _______________________________________________
>>> Rcpp-devel mailing list
>>> [email protected]
>>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>>>
>>>
>
// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
typedef Eigen::MappedSparseMatrix<double> MSpMat;
typedef MSpMat::InnerIterator InIter;
//[[Rcpp::export]]
bool foo (const MSpMat X, const Eigen::Map<Eigen::VectorXi> idx){
int n = X.cols();
if (X.rows() != n) throw std::invalid_argument("Sparse matrix X must be square");
for (int i = 0; i < idx.size(); ++i) {
int jj = idx[i] - 1; // 0-based column index
InIter it(X,jj);
Rcpp::Rcout << "i = " << i << ", jj = " << jj << std::endl;
for (int k = 0; k < idx.size(); ++k) {
int kk = idx[k]-1; // 0-based row index to search for
if (kk == jj) continue; // skip positions on the diagonal
Rcpp::Rcout << " kk = " << kk << ", it.row =";
bool foundit = false;
for (; it; ++it) {
Rcpp::Rcout << " " << it.row();
if (it.row() == kk) {
foundit = true;
++it;
break;
}
if (it.row() > kk) return false;
}
if (!foundit) return false;
Rcpp::Rcout << std::endl;
}
}
return true;
}
/*** R
library(Matrix)
M1 <- new("dgCMatrix"
, i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L, 5L, 3L, 4L)
, p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
, Dim = c(6L, 6L)
, Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c", "d", "e", "f"))
, x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
, factors = list()
)
M1
foo(M1, c(2L,3L,4L))
***/
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
