[Rcpp-commits] r2354 - in pkg/Rcpp: . inst/include/Rcpp/vector inst/unitTests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Oct 21 22:10:23 CEST 2010


Author: romain
Date: 2010-10-21 22:10:23 +0200 (Thu, 21 Oct 2010)
New Revision: 2354

Modified:
   pkg/Rcpp/ChangeLog
   pkg/Rcpp/inst/include/Rcpp/vector/Matrix.h
   pkg/Rcpp/inst/include/Rcpp/vector/MatrixColumn.h
   pkg/Rcpp/inst/unitTests/runit.Matrix.R
Log:
faster matrix indexing

Modified: pkg/Rcpp/ChangeLog
===================================================================
--- pkg/Rcpp/ChangeLog	2010-10-21 19:10:37 UTC (rev 2353)
+++ pkg/Rcpp/ChangeLog	2010-10-21 20:10:23 UTC (rev 2354)
@@ -10,6 +10,9 @@
     sugar expression. Faster indexing (caching the pointer instead of calling
     nrow many times)
     
+    * inst/include/Rcpp/vector/Matrix.h: first pass at making matrix indexing 
+    faster
+    
 2010-10-21  Dirk Eddelbuettel  <edd at debian.org>
 
 	* inst/include/Rcpp/Module.h: Reorder instantiation for class_ 

Modified: pkg/Rcpp/inst/include/Rcpp/vector/Matrix.h
===================================================================
--- pkg/Rcpp/inst/include/Rcpp/vector/Matrix.h	2010-10-21 19:10:37 UTC (rev 2353)
+++ pkg/Rcpp/inst/include/Rcpp/vector/Matrix.h	2010-10-21 20:10:23 UTC (rev 2354)
@@ -36,29 +36,36 @@
     
 	Matrix() : VECTOR() {}
 	
-	Matrix(SEXP x) throw(not_compatible) : VECTOR(){
+	Matrix(SEXP x) throw(not_compatible) : VECTOR(), nrows(0) {
 		if( ! ::Rf_isMatrix(x) ) throw not_compatible("not a matrix") ;
 		SEXP y = r_cast<RTYPE>( x ) ;
 		VECTOR::setSEXP( y );
+		nrows = VECTOR::dims()[0] ;
 	}
 	
-	Matrix( const Dimension& dims) throw(not_compatible) : VECTOR() {
+	Matrix( const Dimension& dims) throw(not_compatible) : VECTOR(), nrows(dims[0]) {
 		if( dims.size() != 2 ) throw not_compatible("not a matrix") ;
 		VECTOR::setSEXP( Rf_allocVector( RTYPE, dims.prod() ) ) ;
 		VECTOR::init() ;
 		VECTOR::attr( "dim" ) = dims ;
 	}
 	
-	Matrix( const int& nrows, const int& ncols) : VECTOR( Dimension( nrows, ncols ) ) {}
+	Matrix( const int& nrows_, const int& ncols) : 
+	    VECTOR( Dimension( nrows_, ncols ) ), 
+	    nrows(nrows_)
+	{}
 	
 	template <typename Iterator>
-	Matrix( const int& nrows, const int& ncols, Iterator start ) : VECTOR( start, start + (nrows*ncols) ) {
+	Matrix( const int& nrows_, const int& ncols, Iterator start ) : 
+	    VECTOR( start, start + (nrows_*ncols) ),
+	    nrows(nrows_)
+	{
 		VECTOR::attr( "dim" ) = Dimension( nrows, ncols ) ; 
 	}
 	
-	Matrix( const int& n) : VECTOR( Dimension( n, n ) ) {}
+	Matrix( const int& n) : VECTOR( Dimension( n, n ) ), nrows(n) {}
 	
-	Matrix( const Matrix& other) throw(not_compatible) : VECTOR() {
+	Matrix( const Matrix& other) throw(not_compatible) : VECTOR(), nrows(other.nrows) {
 		SEXP x = other.asSexp() ;
 		if( ! ::Rf_isMatrix(x) ) throw not_compatible("not a matrix") ;
 		VECTOR::setSEXP( x ) ;
@@ -68,21 +75,21 @@
 		SEXP x = other.asSexp() ;
 		if( ! ::Rf_isMatrix(x) ) not_compatible("not a matrix") ;
 		VECTOR::setSEXP( x ) ;
+		nrows = other.nrows ;
 		return *this ;
 	}
 	
 	template <bool NA, typename MAT>
-    Matrix( const MatrixBase<RTYPE,NA,MAT>& other ) : VECTOR() {
-    	int nr = other.nrow() ;
+    Matrix( const MatrixBase<RTYPE,NA,MAT>& other ) : VECTOR(), nrows(other.nrow()) {
     	int nc = other.ncol() ;
-    	SEXP x = PROTECT( Rf_allocVector( RTYPE, nr * nc ) ) ;
+    	SEXP x = PROTECT( Rf_allocVector( RTYPE, nrows * nc ) ) ;
     	SEXP d = PROTECT( Rf_allocVector( INTSXP, 2) ) ;
-    	INTEGER(d)[0] = nr ;
+    	INTEGER(d)[0] = nrows ;
     	INTEGER(d)[1] = nc ;
     	Rf_setAttrib( x, R_DimSymbol, d ) ;
     	RObject::setSEXP( x ) ;
     	UNPROTECT( 2 ) ;
-    	import_matrix_expression<NA,MAT>( other, nr, nc ) ;
+    	import_matrix_expression<NA,MAT>( other, nrows, nc ) ;
 	}
    
 private:
@@ -121,11 +128,25 @@
     inline Proxy operator[]( int i ) const {
     	return static_cast< const Vector<RTYPE>* >( this )->operator[]( i ) ;
     }
-    inline Proxy operator()( int i, int j ) const {
-    	return static_cast< const Vector<RTYPE>* >( this )->operator()( i, j ) ;
+    
+    inline Proxy operator()( const size_t& i, const size_t& j) throw(not_a_matrix,index_out_of_bounds){
+		return static_cast< Vector<RTYPE>* >( this )->operator[]( 
+		    offset( i, j )
+		    ) ;
+	}
+	
+	inline Proxy operator()( const size_t& i, const size_t& j) const throw(not_a_matrix,index_out_of_bounds){
+		return static_cast< const Vector<RTYPE>* >( this )->operator[]( 
+		    offset( i, j )
+		    ) ;
+	}
+	
+private:
+    
+    inline int offset( int i, int j) const {
+        return i + nrows * j ;
     }
     
-private:
 	virtual void update(){
 		RCPP_DEBUG_1( "%s::update", DEMANGLE(Matrix) ) ;
 		VECTOR::update_vector() ;
@@ -161,13 +182,13 @@
     	return VECTOR::dims()[1]; 
     }
     inline int nrow() const throw(not_a_matrix){
-    	return VECTOR::dims()[0]; 
+    	return nrows ; 
     }
     inline int cols() const throw(not_a_matrix){ 
     	return VECTOR::dims()[1]; 
     }
     inline int rows() const throw(not_a_matrix){ 
-    	return VECTOR::dims()[0]; 
+    	return nrows ; 
     }
 	    
 	typedef MatrixRow<RTYPE> Row ;
@@ -178,8 +199,10 @@
 	
 	inline iterator begin() const{ return VECTOR::begin() ; }
 	inline iterator end() const{ return VECTOR::end() ; }
-	 
     
+private:
+    int nrows ;	
+	
 } ;
 
 #endif

Modified: pkg/Rcpp/inst/include/Rcpp/vector/MatrixColumn.h
===================================================================
--- pkg/Rcpp/inst/include/Rcpp/vector/MatrixColumn.h	2010-10-21 19:10:37 UTC (rev 2353)
+++ pkg/Rcpp/inst/include/Rcpp/vector/MatrixColumn.h	2010-10-21 20:10:23 UTC (rev 2354)
@@ -31,7 +31,7 @@
 	typedef typename MATRIX::iterator iterator ;
 	
 	MatrixColumn( MATRIX& object, int i ) : 
-	    parent(object), index(i), start(parent.begin())
+	    parent(object), index(i), start(parent.begin() + i * parent.nrow() )
 	{
 		if( i < 0 || i >= parent.ncol() ) throw index_out_of_bounds() ;
 	}

Modified: pkg/Rcpp/inst/unitTests/runit.Matrix.R
===================================================================
--- pkg/Rcpp/inst/unitTests/runit.Matrix.R	2010-10-21 19:10:37 UTC (rev 2353)
+++ pkg/Rcpp/inst/unitTests/runit.Matrix.R	2010-10-21 20:10:23 UTC (rev 2354)
@@ -199,6 +199,26 @@
     }
 }
 
+
+test.List.column <- function(){
+	funx <- .rcpp.Matrix$runit_Row_Column_sugar
+	x <- matrix( 1:16+.5, nc = 4 )
+	res <- funx( x )
+	target <- list( 
+	    x[1,], 
+	    x[,1], 
+	    x[2,],
+	    x[,2], 
+	    x[2,] + x[,2]
+	    )
+	cat( "\n" )
+	print( x )
+	print( target[[4]] )
+	print( res[[4]] )       
+	checkEquals( res, target, msg = "column and row as sugar" )
+	
+}
+
 test.NumericMatrix <- function(){
 	funx <- .rcpp.Matrix$matrix_numeric
 	x <- matrix( 1:16 + .5, ncol = 4 )
@@ -308,19 +328,3 @@
 	
 }
 
-test.List.column <- function(){
-	funx <- .rcpp.Matrix$runit_Row_Column_sugar
-	x <- matrix( 1:16+.5, nc = 4 )
-	res <- funx( x )
-	target <- list( 
-	    x[1,], 
-	    x[,1], 
-	    x[2,],
-	    x[,2], 
-	    x[2,] + x[,2]
-	    )
-	checkEquals( res, target, msg = "colum and row as sugar" )
-	
-}
-
-



More information about the Rcpp-commits mailing list