[Rcpp-devel] Is declaring an argument as a reference the same as declaring it a pointer?

[Douglas Bates]

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.