I am calculating three-dimensional sparse hessian:

F3 <- F$jacfun()$jacfun(sparse = TRUE)$jacfun(sparse = TRUE)
res <- F3(f)

It returns a sparse matrix, which represents the 3-dimensional sparse matrix, but the first two dimensions are flattened into one.

This flattening procedure is very reasonable, but it would be good to know, how exactly is it done.

I had to dig into the code, and in .spjacfun in advector.R I found

P <- as(P, "CsparseMatrix") ## Permutes @x and calculates @i and @p
RangeProj(jac@tape, P@x) ## Permute tape output to match new pattern

So, is it that the input of $jacfun(sparse = TRUE) is flattened in column-by-column order?

Am I reading the code right? If so, please add this to the documentation of $jacfun(sparse = TRUE) to ?Tape help page.

Perhaps something like that applies to the non-sparse case?