Re: [Rcpp-devel] RcppArmadillo, solve always requires zgels?
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@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. | 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_bRNS4_IT_EES8_RKNS_4BaseIS6_T0_EE | [bool arma::auxlib::solve_ud, arma::Mat > >(arma::Mat >&, arma::Mat | >&, arma::Base, arma::Mat > > 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@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 ___ 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
Re: [Rcpp-devel] RcppArmadillo, solve always requires zgels?
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 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). On 1 June 2015 at 01:12, Dirk Eddelbuettel 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_bRNS4_IT_EES8_RKNS_4BaseIS6_T0_EE > | [bool arma::auxlib::solve_ud, arma::Mat | > >(arma::Mat >&, > arma::Mat > | >&, arma::Base, arma::Mat > > > 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,
Re: [Rcpp-devel] RcppArmadillo, solve always requires zgels?
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 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_bRNS4_IT_EES8_RKNS_4BaseIS6_T0_EE | | [bool arma::auxlib::solve_ud, arma::Mat > >(arma::Mat >&, arma::Mat | | >&, arma::Base, arma::Mat > > | const | | &)]+0x3a5): undefined reference to `zgels_' | | | | So it thinks it shoul
Re: [Rcpp-devel] RcppArmadillo, solve always requires zgels?
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 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 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_z