[Rcpp-devel] RcppArmadillo, solve always requires zgels?

Dirk Eddelbuettel edd at debian.org
Sun May 31 15:12:34 CEST 2015


Hi Baptiste,

On 31 May 2015 at 20:47, Baptiste Auguie wrote:
| Dear list,
| 
| My cda package (https://github.com/baptiste/cda) solves a linear system Ax=b
| with Armadillo; win-builder and CRAN complain about a missing "zgels_" when
| trying to build on windows or mac, as it's not part of the subset of LAPACK
| provided by R libraries (I asked about this in 2011). I have included a copy of

Right. When no system lapack/blas is found, R uses its own. Which is a
subset.  We were bitten that once and now check in configure for availability
of zgesdd.  I suppose this case is similar.

Look at the logic in

  https://github.com/RcppCore/RcppArmadillo/blob/master/configure

You may need to do something similar for zgels.  So this part of the problem
("do we have zgels ?") really may be part of "which type of R installation do
we have".

Worst case we need new branching / #ifdef logic --- but right now I think
maybe more at the level of your package then at RcppArmadillo.

| this file (zgels.f) from Netlib for the past few years just to make the problem
| go away, but I've recently realised that zgels is not supposed to be called! It
| builds fine on Linux, or any machine with a full LAPACK installed, but I'm now
| wondering if solve() always calls zgels (that would be inefficient for my
| use-case). 

>From the top-level, when I just use ag (a fantastic and much faster grep
alternative; and mail will drop the ansi code for colour highlihting from the
console)

edd at max:~/git/rcpparmadillo(master)$ ag zgels
inst/include/armadillo_bits/lapack_bones.hpp
97:  #define arma_zgels  zgels
205:  #define arma_zgels  ZGELS
340:  void arma_fortran(arma_zgels)(char* trans, blas_int* m, blas_int* n, blas_int* nrhs, void*   a, blas_int* lda, void*   b, blas_int* ldb, void*   work, blas_int* lwork, blas_int* info);

inst/include/armadillo_bits/lapack_wrapper.hpp
689:      arma_fortran(arma_zgels)(trans, m, n, nrhs, (T*)a, lda, (T*)b, ldb, (T*)work, lwork, info);
edd at max:~/git/rcpparmadillo(master)$

I think you should discuss the merits of what gets called with Conrad.

| The matrix A is square: the system is not under/overdetermined, and I believe
| Armadillo should select LU factorisation rather than (slower) QR. Here's the

It is hard to answer this in full generality --- and even harder to change an
existing API.

If Armadillo has been using QR here for years, we cannot really change this
now on the request of a single package.

| error I get from win-builder when I remove zgels.f from my cda package,
| 
| cd.o:cd.cpp:
| (.text$_ZN4arma6auxlib8solve_udISt7complexIdENS_3MatIS3_EEEEbRNS4_IT_EES8_RKNS_4BaseIS6_T0_EE
| [bool arma::auxlib::solve_ud<std::complex<double>, arma::Mat<std::complex
| <double> > >(arma::Mat<std::complex<double> >&, arma::Mat<std::complex<double>
| >&, arma::Base<std::complex<double>, arma::Mat<std::complex<double> > > const
| &)]+0x3a5): undefined reference to `zgels_'
| 
| So it thinks it should call solve_ud, which is meant for under-determined
| systems. Why would that be? 
| I tried to find the relevant piece in the template code, this looks relevant
| but I don't understand it: https://github.com/RcppCore/RcppArmadillo/blob/
| master/inst/include/armadillo_bits/glue_solve_meat.hpp#L26

Well that seems to be the call in the 'rows < cols' case.

But if there is a function solve_ud, can you not just call solve_ud
from your package?

I haven't had coffee yet so I may not need the full picture.

Cheers, Dirk



| 
| Best regards,
| 
| baptiste
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

-- 
http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org


More information about the Rcpp-devel mailing list