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 <e...@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@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@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@lists.r-forge.r-project.org | | https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel | | -- | http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org | | -- http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org _______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
