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

Baptiste Auguie baptiste.auguie at gmail.com
Mon Jun 1 01:19:59 CEST 2015


Hi,

I attach a minimal sourceCpp file:

Rcpp::sourceCpp('testzgels.cpp')

provides extra debugging info, and shows my complete lack of understanding
of c++ templates. It would appear that upon compilation, all cases (nrow ==
ncol, nrow > ncol, nrow < ncol) are considered (of course: the matrix A has
no size information for the compiler). The size check is only performed
when calling the function, and dispatches to the correct LAPACK routine
(gels vs gesv).

Now, unfortunately this example does not help me solve a more specialised
case, where **I know** only square matrices will be passed. I do not need
the various possible routines, which in fact won't be all available within
R's subset of LAPACK.
zgels is assumed to be required (by CRAN and win-builder), presumably
because the templates cater for all possible cases and the code-checking
just looks up the list of symbol names, although in actuality it would
never be used. Does that make sense? How could I specialise the code to
restrict solve to square systems?

Thanks,

baptiste






On 1 June 2015 at 10:13, Dirk Eddelbuettel <edd at debian.org> wrote:

>
> On 1 June 2015 at 09:06, Baptiste Auguie wrote:
> | Hi,
> |
> | Thanks Dirk for your suggestions. Let me try to clarify a few things
> below.
> | Basically,
> |
> | - zgels should not be needed (probably, and if I was doing things
> correctly)
> |
> | - I believe the template mechanism is simply unaware of the square
> format of my
> | matrix, and defaults to the nrow < ncol case for lack of better
> knowledge. I do
>
> Per this line
>
>
> https://github.com/RcppCore/RcppArmadillo/blob/master/inst/include/armadillo_bits/glue_solve_meat.hpp#L26
>
> I somehow have my doubts.
>
> | not know how to specify that the matrix is square, it's passed on from
> another
> | function, i.e.
> |
> | https://github.com/baptiste/cda/blob/master/src/cd.cpp#L49
> | is using A from
> | https://github.com/baptiste/cda/blob/master/src/cda.cpp#L129
> |
> | The way I picture it, I need to give a hint of the square size, and
> armadillo
> | will pick the right function/template/glue (I don't know the mechanism at
> | play).
>
> Can you minimize the problem?  "Minimally reproducible example" which still
> gets to zgels when maybe it shouldn't?
>
> Dirk
>
> |
> | On 1 June 2015 at 01:12, Dirk Eddelbuettel <edd at debian.org> wrote:
> |
> |
> |     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.
> |
> |
> |
> | The easiest workaround, I find, is to include zgels.f, as I've been
> doing since
> | 2011. Now, my premise here is that this function is not needed in my
> package.
> |
> |
> |     | 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)
> |
> |
> | Yes, I've tried to look at where zgels is called in the original
> Armadillo;
> | it's not crystal clear (for me) with the templates and glue etc.
> |
> |
> |     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.
> |
> |
> | Looking at the code and docs, the intent for square matrices is probably
> not to
> | use zgels. I believe my code is simply failing to give enough hints for
> the
> | templates to work out which case is actually being required.
> |
> |
> |     | 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.
> |
> |
> | I linked to the case which is returned in the win-builder log; the case
> of a
> | square matrix would be a few lines earlier (note that the solve_ud case
> is the
> | default, if previous conditions fail to apply).
> |
> |
> |
> |     But if there is a function solve_ud, can you not just call solve_ud
> |     from your package?
> |
> |
> | That would be auxlib::solve(out, A, X, slow); Yes, it's probably an
> option,
> | although if I were to bypass the high-level solve() interface, I'd
> probably
> | feel more comfortable with an explicit call to lu() instead.
> |
> |
> |     I haven't had coffee yet so I may not need the full picture.
> |
> |     Cheers, Dirk
> |
> |
> | Thanks!
> |
> | baptiste
> |
> |
> |
> |
> |     |
> |     | 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
> |
> |
>
> --
> http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150601/c65586f2/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testzgels.cpp
Type: text/x-c++src
Size: 608 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150601/c65586f2/attachment-0001.cpp>


More information about the Rcpp-devel mailing list