CODEWITHSUNDEEP
rmatrixrcpp
aenima
aenima

Reputation: 367

RcppArmadillo: Complex matrix inverse compilation error

I am trying to invert a complex square matrix using RcppArmadillo:

# include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
cx_mat fn(cx_mat x) {
  return x.i();
}

However, this throws an error when I sourceCpp it: "undefined reference to zgetri_'". It compiles and works fine if I just replace cx_mat with mat, but then it will only work with real matrices. Using inv() throws the same compilation error. Interestingly, using the pseudo-inverse pinv() passes compilation, but the result would be slightly different compared with R's solve():

# include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
cx_mat fn(cx_mat x) {
  return pinv(x);
}

> a<-matrix(c(3+0.1i,7+2i,4+0i,2+0.5i),ncol=2)
> a
       [,1]   [,2]
[1,] 3+0.1i 4+0.0i
[2,] 7+2.0i 2+0.5i
> identical(solve(a), fn(a))
[1] FALSE
> solve(a) - fn(a)
                           [,1]                       [,2]
[1,] -6.938894e-17-7.80626e-18i  8.326673e-17-6.93889e-18i
[2,]  1.665335e-16+1.95156e-17i -1.665335e-16+4.16334e-17i

I know the difference is flirting with machine precision in this case, but I am still wondering if there is anyway to get inv() or .i() to work on complex matrices. Thanks.

Upvotes: 3

Views: 601

Answers (1)

Dirk is no longer here
Dirk is no longer here

Reputation: 368251

That is a known issue if and when you use RcppArmadillo with an R installation using Rlapack.so as for example on Windows or some RedHat systems.

The best answers are to either

  • not use complex-valued functions not included as Rlapack has some but sadly not all (but you may not have a choice)
  • reconfigure you R installation to use an full external liblapack.so

We have three open tickets at the RcppArmadillo repo (and in fact I even added one today listing the twelve missing complex functions used by Armadillo but missing in Rlapack.so), and I just pleaded with R Core to add more complex-valued functions to Rlapack.

And just to stress the second point, your example works fine here as I do not use Rlapack on the Debian/Ubuntu builds:

R> library(Rcpp)
R> sourceCpp("/tmp/aenima.cpp")

R> a <- matrix(c(3+0.1i,7+2i,4+0i,2+0.5i),ncol=2)

R> fn(a)
                      [,1]                [,2]
[1,] -0.0898473+0.0029949i  0.167715-0.047919i
[2,]  0.3174603+0.0000000i -0.126984+0.031746i
R>

using a slightly modified version of your file with the example at the end:

# include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
cx_mat fn(cx_mat x) {
  return x.i();
}

/*** R
a <- matrix(c(3+0.1i,7+2i,4+0i,2+0.5i),ncol=2)
fn(a)
*/

Upvotes: 7

Related Questions