I am trying to implement some basic C++ code for big.matrix objects in R. I am using the Rcpp package, have read the demo here and even applied another simple function which I found on the rcpp-devel list:
#include "bigmemory/BigMatrix.h"
#include "bigmemory/MatrixAccessor.hpp"
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
void fun(SEXP A) {
Rcpp::XPtr<BigMatrix> bigMat(A);
MatrixAccessor<int> Am(*bigMat);
int nrows = bigMat->nrow();
int ncolumns = bigMat->ncol();
for (int j = 0; j < ncolumns; j++){
for (int i = 1; i < nrows; i++){
Am[j][i] = Am[j][i] + Am[j][i-1];
}
}
return;
}
// [[Rcpp::export]]
void BigTranspose(SEXP A)
{
Rcpp::XPtr<BigMatrix> pMat(A);
MatrixAccessor<int> mat(*pMat);
int r = pMat->nrow();
int c = pMat->ncol();
for(int i=0; i<r; ++i)
for(int j=0; j<c; ++j)
std::swap(mat[j][i], mat[i][j]);
return;
}
This fun function works perfectly fine, modifying the big.matrix object.
a <- matrix(seq(25), 5,5)
> a
[,1] [,2] [,3] [,4] [,5]
[1,] 1 6 11 16 21
[2,] 2 7 12 17 22
[3,] 3 8 13 18 23
[4,] 4 9 14 19 24
[5,] 5 10 15 20 25
> fun(b@address)
> head(b)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 6 11 16 21
[2,] 3 13 23 33 43
[3,] 6 21 36 51 66
[4,] 10 30 50 70 90
[5,] 15 40 65 90 115
However, when I try a simple square matrix transpose function the matrix is not modified. Why would the fun function work but not my 'BigTranspose`?
a <- matrix(seq(25), 5,5)
> a
[,1] [,2] [,3] [,4] [,5]
[1,] 1 6 11 16 21
[2,] 2 7 12 17 22
[3,] 3 8 13 18 23
[4,] 4 9 14 19 24
[5,] 5 10 15 20 25
b <- as.big.matrix(a)
BigTranspose(b@address)
> head(b)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 6 11 16 21
[2,] 2 7 12 17 22
[3,] 3 8 13 18 23
[4,] 4 9 14 19 24
[5,] 5 10 15 20 25
The problem is in your transpose algorithm; you are looping over all rows and columns, so your swaps are being re-swapped. You can fix this by just setting j's lower index to j = i+1 instead of j = 0:
#include "bigmemory/BigMatrix.h"
#include "bigmemory/MatrixAccessor.hpp"
#include <Rcpp.h>
// [[Rcpp::depends(BH, bigmemory)]]
using namespace Rcpp;
// [[Rcpp::export]]
void BigTranspose(SEXP A)
{
Rcpp::XPtr<BigMatrix> pMat(A);
MatrixAccessor<int> mat(*pMat);
int r = pMat->nrow();
int c = pMat->ncol();
for(int i=0; i<r; i++){
for(int j=0; j<c; j++){
std::swap(mat[j][i], mat[i][j]);
}
}
return;
}
// [[Rcpp::export]]
void BigTranspose2(SEXP A)
{
Rcpp::XPtr<BigMatrix> pMat(A);
MatrixAccessor<int> mat(*pMat);
int r = pMat->nrow();
int c = pMat->ncol();
for(int i=0; i<r; i++){
for(int j=(i+1); j<c; j++){
std::swap(mat[j][i], mat[i][j]);
}
}
return;
}
/*** R
##
b1 <- as.big.matrix(a)
head(b1)
##
BigTranspose(b1@address)
head(b1)
##
##
b2 <- as.big.matrix(a)
head(b2)
##
BigTranspose2(b2@address)
head(b2)
##
##
M <- matrix(1:25,ncol=5)
t(M)
##
*/
You can see that the second version works correctly by comparing the output of BigTranspose2 with that of t(M), where M is the matrix equivalent of b1 and b2:
> Rcpp::sourceCpp('bigMatTranspose.cpp')
> b1 <- as.big.matrix(a)
> head(b1)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 6 11 16 21
[2,] 2 7 12 17 22
[3,] 3 8 13 18 23
[4,] 4 9 14 19 24
[5,] 5 10 15 20 25
> ##
> BigTranspose(b1@address)
> head(b1)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 6 11 16 21
[2,] 2 7 12 17 22
[3,] 3 8 13 18 23
[4,] 4 9 14 19 24
[5,] 5 10 15 20 25
> ##
> ##
> b2 <- as.big.matrix(a)
> head(b2)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 6 11 16 21
[2,] 2 7 12 17 22
[3,] 3 8 13 18 23
[4,] 4 9 14 19 24
[5,] 5 10 15 20 25
> ##
> BigTranspose2(b2@address)
> head(b2)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 6 7 8 9 10
[3,] 11 12 13 14 15
[4,] 16 17 18 19 20
[5,] 21 22 23 24 25
> ##
> ##
> M <- matrix(1:25,ncol=5)
> t(M)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 6 7 8 9 10
[3,] 11 12 13 14 15
[4,] 16 17 18 19 20
[5,] 21 22 23 24 25
Just note that this will work correctly for square matrices, but you will have to make some modifications for it to work with matrices of arbitrary dimension.