How to make a Lapack routine accessible to the Csnippet

[etsagued]

Dear all,

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")

Call the C function

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

test the 'dmeasure':

replicate(n=10, rick4 |> pfilter(Np=10000) |> logLik() ) |> logmeanexp(se=TRUE)
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").

Thank you!
Edouard
sylvester.zip