[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 <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