Skip to content

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

7 messages · Dirk Eddelbuettel, Douglas Bates, Whit Armstrong +1 more

Original
Douglas Bates
# 
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
/* 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(), ...

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.
Dirk Eddelbuettel
# 
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.

Dirk
Douglas Bates
#
On Sat, Mar 27, 2010 at 6:12 PM, Dirk Eddelbuettel <edd at debian.org> wrote:
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.
Whit Armstrong
# 
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:

http://osl.iu.edu/research/mtl/mtl4/doc/matrix__vector__expr.html
http://arma.sourceforge.net/docs.html#operators

My secret desire is to implement R using Armadillo as the core data
structures...

just my 2c.

-Whit
On Sat, Mar 27, 2010 at 7:22 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
Dirk Eddelbuettel
#
On 27 March 2010 at 19:56, Whit Armstrong wrote:
| 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:
| 
| http://osl.iu.edu/research/mtl/mtl4/doc/matrix__vector__expr.html
| http://arma.sourceforge.net/docs.html#operators

And as Romain pointed out before the template programming means that we're
saving temporary objects.  So it can be cleaner to read and be more
efficient, now that compilers have caught up.
 
| My secret desire is to implement R using Armadillo as the core data
| structures...

Ha :)  Talk to Andrew Runnalls about how his CXXR reimplementation is coming
along.  He has a new object structure in there but alas no sign of Armadillo.

Dirk
Douglas Bates
# 
On Sat, Mar 27, 2010 at 6:56 PM, Whit Armstrong
<armstrong.whit at gmail.com> wrote:
I would plan to use the matrix templates available in Rcpp or
Armadillo but I still don't like the approach taken by Armadillo for
matrix decompositions.  From my reading they pattern their operations
on Matlab and, perhaps because my involvement in numerical linear
algebra predates Matlab, I tend to think in terms of more basic
operations.  For something like the QR decomposition I prefer to get
access to the Householder transformations and not matrices derived
from them.

I would write wrappers for all the Lapack routines that I use.  What I
am considering here is how the wrappers can be written.
1 day later
Romain Francois
# 
Le 27/03/10 23:10, Douglas Bates a ?crit :
no. you have to do what you do above. references are not pointers. I 
guess you can do thing like :

F77_CALL(dgemm)(transa, transb,&(A.nrow()),&(A.ncol()), ...)

but I'd advise to stick to your current code which is easier to understand.
I don't think so (I have no clue about fortran however).