GSoC ode15{i,s} - type conversion and approximation problems

[Francesco Faccio]

Dear all,

I wrote a minimal version of ode15i that can solve the first benchmark problem I described here:
http://gsoc2016ode15s.blogspot.it/

The code I wrote is in this repository:
https://bitbucket.org/Francesco_Faccio/octave/src/f6fe4b50fab06b810245638362085a6a9245cd69/libinterp/dldfcn/ode15i.cc?at=default&fileviewer=file-view-default

I would  like to ask you some suggestions regarding some problems on the implementation.

Sundials uses its own types: realtype and N_Vector. An N_Vector is a vector of realtype, while a realtype can be both a float, a double or a long double, depending on how Sundials has been built (default type is double).

How can I deal with this? In the documentation, it seems that octave_value doesn't have a type like long double, so any conversion from N_Vector to ColumnVector would loose precision. Maybe I should just check sizeof (realtype) == sizeof (double) and print an error message if it's not true. 

The second thing I would like to ask regards dense/sparse Jacobian matrices.

Using dense methods, if Jacobian matrix is not supplied, Sundials approximates it with different quotients. If I want to pass a Jacobian to Sundials' solvers, I have to pass both DF/DY and DF/DYP. If the user gives in input only  [[ ], DF/DYP] or [DF/DY, [ ]], I have to approximate DF/DY or DF/DYP before passing them to Sundials' solver. Matlab in this case approximates them by finite differences: should I do the same or maybe use another external library to do it (Sundials IDA doesn't give a direct interface for its approximation functions)?

Using sparse methods, in order to use Sundials' sparse solvers, I must provide both DF/DY and DF/DYP in CSC format (it cannot be done automatically with different quotient by Sundials), so I need to approximate one of them if the user passes only one matrix. Again, should I write my own approximation or use an external library (and which could be better to use for this problem)?

Finally, I have to provide to Sundials's solvers a function which computes F(t, y(t), y'(t)), given N_Vector y, y' at time realtype t. I wrote it using a global pointer to the function F, converting y, y', t from N_Vector to ColumnVector, aggregating them in an octave_value_list variable, calling feval, and finally converting the output to an N_Vector. This works, but I think it would waste a lot of time for solving large problems. Should be a good idea to overload feval in order to use N_Vector types without any conversion?

Thank you very much.

Best regards,

Francesco Faccio