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 <[email protected]> wrote: > On Sat, Mar 27, 2010 at 6:12 PM, Dirk Eddelbuettel <[email protected]> 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 > [email protected] > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel > _______________________________________________ Rcpp-devel mailing list [email protected] https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
