<div dir="ltr">Hi,<div><br></div><div>I attach a minimal sourceCpp file:</div><div><br></div><div>Rcpp::sourceCpp('testzgels.cpp')<br></div><div><br></div><div>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). </div><div><br></div><div>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. </div><div>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?</div><div><br></div><div>Thanks,</div><div><br></div><div>baptiste</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On 1 June 2015 at 10:13, Dirk Eddelbuettel <span dir="ltr"><<a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class=""><br>
On 1 June 2015 at 09:06, Baptiste Auguie wrote:<br>
| Hi,<br>
|<br>
| Thanks Dirk for your suggestions. Let me try to clarify a few things below.<br>
| Basically,<br>
|<br>
| - zgels should not be needed (probably, and if I was doing things correctly)<br>
|<br>
| - I believe the template mechanism is simply unaware of the square format of my<br>
| matrix, and defaults to the nrow < ncol case for lack of better knowledge. I do<br>
<br>
</span>Per this line<br>
<br>
  <a href="https://github.com/RcppCore/RcppArmadillo/blob/master/inst/include/armadillo_bits/glue_solve_meat.hpp#L26" target="_blank">https://github.com/RcppCore/RcppArmadillo/blob/master/inst/include/armadillo_bits/glue_solve_meat.hpp#L26</a><br>
<br>
I somehow have my doubts.<br>
<span class=""><br>
| not know how to specify that the matrix is square, it's passed on from another<br>
| function, i.e.<br>
|<br>
| <a href="https://github.com/baptiste/cda/blob/master/src/cd.cpp#L49" target="_blank">https://github.com/baptiste/cda/blob/master/src/cd.cpp#L49</a> <br>
| is using A from<br>
| <a href="https://github.com/baptiste/cda/blob/master/src/cda.cpp#L129" target="_blank">https://github.com/baptiste/cda/blob/master/src/cda.cpp#L129</a><br>
|<br>
| The way I picture it, I need to give a hint of the square size, and armadillo<br>
| will pick the right function/template/glue (I don't know the mechanism at<br>
| play).<br>
<br>
</span>Can you minimize the problem?  "Minimally reproducible example" which still<br>
gets to zgels when maybe it shouldn't?<br>
<span class="HOEnZb"><font color="#888888"><br>
Dirk<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
|<br>
| On 1 June 2015 at 01:12, Dirk Eddelbuettel <<a href="mailto:edd@debian.org">edd@debian.org</a>> wrote:<br>
|<br>
|<br>
|     Hi Baptiste,<br>
|<br>
|     On 31 May 2015 at 20:47, Baptiste Auguie wrote:<br>
|     | Dear list,<br>
|     |<br>
|     | My cda package (<a href="https://github.com/baptiste/cda" target="_blank">https://github.com/baptiste/cda</a>) solves a linear system<br>
|     Ax=b<br>
|     | with Armadillo; win-builder and CRAN complain about a missing "zgels_"<br>
|     when<br>
|     | trying to build on windows or mac, as it's not part of the subset of<br>
|     LAPACK<br>
|     | provided by R libraries (I asked about this in 2011). I have included a<br>
|     copy of<br>
|<br>
|     Right. When no system lapack/blas is found, R uses its own. Which is a<br>
|     subset.  We were bitten that once and now check in configure for<br>
|     availability<br>
|     of zgesdd.  I suppose this case is similar.<br>
|<br>
|     Look at the logic in<br>
|<br>
|       <a href="https://github.com/RcppCore/RcppArmadillo/blob/master/configure" target="_blank">https://github.com/RcppCore/RcppArmadillo/blob/master/configure</a><br>
|<br>
|     You may need to do something similar for zgels.  So this part of the<br>
|     problem<br>
|     ("do we have zgels ?") really may be part of "which type of R installation<br>
|     do<br>
|     we have".<br>
|<br>
|     Worst case we need new branching / #ifdef logic --- but right now I think<br>
|     maybe more at the level of your package then at RcppArmadillo.<br>
|<br>
|<br>
|<br>
| The easiest workaround, I find, is to include zgels.f, as I've been doing since<br>
| 2011. Now, my premise here is that this function is not needed in my package.<br>
|  <br>
|<br>
|     | this file (zgels.f) from Netlib for the past few years just to make the<br>
|     problem<br>
|     | go away, but I've recently realised that zgels is not supposed to be<br>
|     called! It<br>
|     | builds fine on Linux, or any machine with a full LAPACK installed, but<br>
|     I'm now<br>
|     | wondering if solve() always calls zgels (that would be inefficient for my<br>
|     | use-case). <br>
|<br>
|     From the top-level, when I just use ag (a fantastic and much faster grep<br>
|     alternative; and mail will drop the ansi code for colour highlihting from<br>
|     the<br>
|     console)<br>
|<br>
|<br>
| Yes, I've tried to look at where zgels is called in the original Armadillo;<br>
| it's not crystal clear (for me) with the templates and glue etc.  <br>
|<br>
|<br>
|     edd@max:~/git/rcpparmadillo(master)$ ag zgels<br>
|     inst/include/armadillo_bits/lapack_bones.hpp<br>
|     97:  #define arma_zgels  zgels<br>
|     205:  #define arma_zgels  ZGELS<br>
|     340:  void arma_fortran(arma_zgels)(char* trans, blas_int* m, blas_int* n,<br>
|     blas_int* nrhs, void*   a, blas_int* lda, void*   b, blas_int* ldb, void* <br>
|      work, blas_int* lwork, blas_int* info);<br>
|<br>
|     inst/include/armadillo_bits/lapack_wrapper.hpp<br>
|     689:      arma_fortran(arma_zgels)(trans, m, n, nrhs, (T*)a, lda, (T*)b,<br>
|     ldb, (T*)work, lwork, info);<br>
|     edd@max:~/git/rcpparmadillo(master)$<br>
|<br>
|     I think you should discuss the merits of what gets called with Conrad.<br>
|<br>
|<br>
| Looking at the code and docs, the intent for square matrices is probably not to<br>
| use zgels. I believe my code is simply failing to give enough hints for the<br>
| templates to work out which case is actually being required.<br>
|<br>
|<br>
|     | The matrix A is square: the system is not under/overdetermined, and I<br>
|     believe<br>
|     | Armadillo should select LU factorisation rather than (slower) QR. Here's<br>
|     the<br>
|<br>
|     It is hard to answer this in full generality --- and even harder to change<br>
|     an<br>
|     existing API.<br>
|<br>
|     If Armadillo has been using QR here for years, we cannot really change this<br>
|     now on the request of a single package. <br>
|<br>
|<br>
|     | error I get from win-builder when I remove zgels.f from my cda package,<br>
|     |<br>
|     | cd.o:cd.cpp:<br>
|     |<br>
|     (.text$_ZN4arma6auxlib8solve_udISt7complexIdENS_3MatIS3_EEEEbRNS4_IT_EES8_RKNS_4BaseIS6_T0_EE<br>
|     | [bool arma::auxlib::solve_ud<std::complex<double>, arma::Mat<std::complex<br>
|     | <double> > >(arma::Mat<std::complex<double> >&, arma::Mat<std::complex<br>
|     <double><br>
|     | >&, arma::Base<std::complex<double>, arma::Mat<std::complex<double> > ><br>
|     const<br>
|     | &)]+0x3a5): undefined reference to `zgels_'<br>
|     |<br>
|     | So it thinks it should call solve_ud, which is meant for under-determined<br>
|     | systems. Why would that be? <br>
|     | I tried to find the relevant piece in the template code, this looks<br>
|     relevant<br>
|     | but I don't understand it: <a href="https://github.com/RcppCore/RcppArmadillo/blob" target="_blank">https://github.com/RcppCore/RcppArmadillo/blob</a><br>
|     /<br>
|     | master/inst/include/armadillo_bits/glue_solve_meat.hpp#L26<br>
|<br>
|     Well that seems to be the call in the 'rows < cols' case.<br>
|<br>
|<br>
| I linked to the case which is returned in the win-builder log; the case of a<br>
| square matrix would be a few lines earlier (note that the solve_ud case is the<br>
| default, if previous conditions fail to apply).<br>
|  <br>
|<br>
|<br>
|     But if there is a function solve_ud, can you not just call solve_ud<br>
|     from your package?<br>
|<br>
|<br>
| That would be auxlib::solve(out, A, X, slow); Yes, it's probably an option,<br>
| although if I were to bypass the high-level solve() interface, I'd probably<br>
| feel more comfortable with an explicit call to lu() instead. <br>
|<br>
|<br>
|     I haven't had coffee yet so I may not need the full picture. <br>
|<br>
|     Cheers, Dirk<br>
|<br>
|  <br>
| Thanks!<br>
|<br>
| baptiste <br>
|<br>
|<br>
|<br>
|<br>
|     |<br>
|     | Best regards,<br>
|     |<br>
|     | baptiste<br>
|     | _______________________________________________<br>
|     | Rcpp-devel mailing list<br>
|     | <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
|     | <a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
|<br>
|     --<br>
|     <a href="http://dirk.eddelbuettel.com" target="_blank">http://dirk.eddelbuettel.com</a> | @eddelbuettel | <a href="mailto:edd@debian.org">edd@debian.org</a><br>
|<br>
|<br>
<br>
--<br>
<a href="http://dirk.eddelbuettel.com" target="_blank">http://dirk.eddelbuettel.com</a> | @eddelbuettel | <a href="mailto:edd@debian.org">edd@debian.org</a><br>
</div></div></blockquote></div><br></div>