[Rcpp-devel] Is declaring an argument as a reference the same as declaring it a pointer?
armstrong.whit at gmail.com
Sun Mar 28 00:56:53 CET 2010
All the newer c++ matrix template libraries are very good. I like
Armadillo in particular with MTL as a close second.
As much of a pain as it is to learn a new library, you code will be
infinitely more readable if you use a matrix library.
here are the operator examples from MTL4 and Armadillo for reference:
My secret desire is to implement R using Armadillo as the core data
just my 2c.
On Sat, Mar 27, 2010 at 7:22 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
> On Sat, Mar 27, 2010 at 6:12 PM, Dirk Eddelbuettel <edd at debian.org> wrote:
>> That's a long with a number of questions, I'll just cherry-pick:
>> On 27 March 2010 at 17:10, Douglas Bates wrote:
>> | I am writing code to interface from Rcpp to the Lapack linear algebra
>> | code which is written in Fortran.
>> | Because Lapack is written in Fortran all arguments must be passed from
>> | C as pointers. There are many examples of the C versions of the
>> | declarations in the file R_HOME/include/R_ext/Lapack.h where you will
>> | see that a common idiom is to pass a pointer to the contents of a
>> | matrix, a pointer to an integer giving the number of columns, the same
>> | for the number of rows, and a third pointer which is the "leading
>> | dimension of the array as declared in the calling program". (If all
>> | this seems bizarre, be grateful that you never needed to learn to
>> | program in Fortran.)
>> | For example, the C declaration of dgemm, which computes C := alpha * A
>> | %*% B + beta * C, is
>> (Everybody's favourite Blas-3 function :)
>> | /* DGEMM - perform one of the matrix-matrix operations */
>> | /* C := alpha*op( A )*op( B ) + beta*C */
>> | void
>> | F77_NAME(dgemm)(const char *transa, const char *transb, const int *m,
>> | const int *n, const int *k, const double *alpha,
>> | const double *a, const int *lda,
>> | const double *b, const int *ldb,
>> | const double *beta, double *c, const int *ldc);
>> | The arguments m, n, k, alpha, lda, ldb, and beta are all int's or
>> | double's that must be passed as pointers. If A, B and C are
>> | NumericMatrix objects then I need to write code like
>> | const char *transa = "N", transb = "N";
>> | const int m = A.nrow(), n = A.ncol(), ...
>> | const double one = 1.0, zero = 0.0;
>> | F77_CALL(dgemm)(transa, transb, &m, &n, &one, A.begin(), &m, ...)
>> | which is somewhat roundabout because I need to declare int's and
>> | double's and initialize them to fixed values so I can pass those fixed
>> | values as pointers for the Fortran routine.
>> | Can I declare a C++ interface replacing the const int *m, const double
>> | *alpha, etc. by
>> | const int &m, const double &alpha, etc. so that I can call the Fortran
>> | subroutine as
>> | F77_CALL(dgemm)(transa, transb, A.nrow(), A.ncol(), 1.0,
>> | A.begin(), A.nrow(), ...
>> I played with something related the other day when I (re-)worked the two MPI
>> examples contributed to RInside. I didn't end up getting this sorted out. It
>> seems one really needs the "C approach" here with 'const int m = N.nrow()'
>> and passing &m.
>> | In other words, do the formal argument declarations const int &m and
>> | const int *m both end up passing an address of an int, with the only
>> | difference being in what the form of the actual argument is?
>> | If so, is there any way to avoid the dance with the character strings?
>> | It turns out that only the first character is ever examined in the
>> | Lapack code by I still need to pass an address to that character, not
>> | the character itself.
>> I fear it may be the same issue, so you may need the char.
>> Would the C++ interface to Lapack offer help? Otherwise, you could play
>> "just don't do it" and use a modern shim like Armadillo, Eigen, ... or even
>> uBlas from Boost which all end up calling Lapack functions for you. Or you
>> sweat the distateful bits, write'em once and hide them behind other functions
>> you call.
> I think a C++ interface would help in that it makes it easier to read.
> Because everything in Fortran is passed as an address there is no
> real distinction between scalars and vectors (or pointers in C). I
> wrote the R_HOME/include/R_ext/Lapack.h file and went to some trouble
> to differentiate between a declaration of int *m and int iwork. The
> first would be a scalar (passed as a pointer) and the second would be
> a vector, even though the two declarations are equivalent in C.
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
More information about the Rcpp-devel