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

[Baptiste Auguie]

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 <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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150601/a91963b0/attachment-0001.html>