How to make a Lapack routine accessible to the Csnippet

[etsagued]

I am new to the pomp package and would like to use it to estimate a nonlinear state-space model where the state-space form is a solution to a Sylvester equation of the form AX + XB = C. The main difficulty so far has been to speed up the code by using the C snippets. It is not entirely obvious to me on how to solve the Sylvester equation within the C snippets. I have proceeded in two steps.

Step 1: Call a lapack routine and use it R
Here is the test code

library(pomp)
library(magrittr)
library(maotai)

system("R CMD SHLIB lapack_test.c -llapack")
dyn.load("lapack_test.so")

A <- matrix(c(1, 2, 3, 4), nrow = 2)
B <- matrix(c(1, 0, 0, 1), nrow = 2)
C <- matrix(c(5, 6, 7, 8), nrow = 2)

## solution of sylvester equation from package maotai

(X <- maotai::sylvester(A, B, C))

## check the solution from lapack directly
.C("solve_sylvester",
  A = A,
  B = B,
  C = C,
  m = nrow(A),
  n = nrow(A)
)

.C("solve_sylvester",
  A = as.double(A),
  B = as.double(B),
  C = as.double(C),
  m = as.integer(2),
  n = as.integer(2)
)

Step 2: Integrate Sylvester equation solver within the rprocess
For this I consider the following toy example

library(pomp)
library(magrittr)
library(maotai)

readLines("sylvester.h") |> paste(collapse="\n") -> global.h
readLines("sylvester.c") -> global.c
global.c[grepl("include \"my_sylvester.h\"",global.c)] <- global.h
extract(global.c,!grepl("include \"sylvester.h\"",global.c)) |>
  paste(collapse="\n") -> global.c

ricker() |>
  simulate(rmeasure=Csnippet("
    int n = 2, m=2;
    double A[4] = {1.0, 3.0, 2.0, 4.0};
    double B[4] = {1.0, 0.0, 0.0, 1.0};
    double C[4] = {5.0, 7.0, 6.0, 8.0};
    double X[3][3];
    solve_sylvester(A, B, C, &m, &n);
    y = rnorm(N+C[0],1);
    y = (y > 0) ? y : 0;"),
  dmeasure=Csnippet("
    lik = dnorm(y,N,1,0);
    if (give_log) lik = log(lik);"),
  statenames="N",
  globals=global.c
)  -> rick4

rick4 |>
  pfilter(Np=10000) |>
  logLik() |>
  replicate(n=10) |>
  logmeanexp(se=TRUE)

sylvester.zip

I have not been able to make this toy example to work. I am getting the following error:

Error: in ‘simulate’: unable to load shared object '/tmp/RtmpMQc4ey/833422/pomp_753876c7d31d99a58252ecbdc9d0ea9d.so':
  /tmp/RtmpMQc4ey/833422/pomp_753876c7d31d99a58252ecbdc9d0ea9d.so: undefined symbol: dtrsyl_

Is this error indicating the compilation was not successful? How to link the lapack library during compilation? In step 1 I use system("R CMD SHLIB lapack_test.c -llapack").

[kingaa]

On the one hand, @etsagued, the answer to your question is relatively straightforward. Have a look at FAQ 3.7, which gives instructions on linking C snippets with external libraries.

On the other hand, I wonder whether for such a model it is necessary to use simulation-based methods. Given the linear structure of the problem, more efficient approaches are available, no?

[etsagued]

Thank you very much your promptitude!

I actually use the FAQ 3.7 to create the toy example above. I will reread it again and try to make it work.

Sorry I was not clear in my previous exposition. My model have a linear and nonlinear part. The linear part is a solution to a Sylvester equation which is then combined in a nonlinear with nonadditive shocks to create the full model. That is why, I thought I could use the MIF algorithm.

Thanks!

[etsagued]

I finally made to work. As explained in the example FAQ 3.7, I create a directory "include" that contains "sylvester.h" and another one "libs" with "sylvester.so" in it and everything worked beautifully.

Thank you again for your help!

[kingaa]

Sure thing! I look forward to hearing about your results!