Functional meets Fortran

I decided my life needed a little bit more math, and so to decided to add a few differential equations in some F# code. Math in code means applied math and lots of arrays, but arrays are a pain functional style. So heres my appoach.

I have written many numerial applications, mostly in C++, yet also some in Fortran. I mention Fortran because it has the notion of equivalence, where you can map arrays on to each other.  Equivalence is like C unions for array and there to  reminds us is that arrays do not need to be managed as dynamic memory.

A functional style of programming is a lot about preserving your invariants with immutability, yet immutable vectors are not a solution (yet). Still, numerical code has its own invariants, and therefore some immutability is possible, and should be exploited to write robust code.  In my case, I have vectors and matrices; Their sizes are invariant, and so is the structure of the matrices, they have lots of zeros, they are sparse.  The approach is then to specialize a state monad to maintain size, and structure in an immutable manner, while allocating mutable vectors only within the core of the numerical engine.

Long ago I wrote an matrix assembly algorithm. You can see it as a sparse matrix structure inference code. It was really pretty simple: the final matrix is an assembly of many small sparse matrices, finding the final structure can be seen as solving a set of matrix structure equations, the algorithm defines an order of traversal for all this matrix points and then proceeds to "crawl" across all matrices at once, working out all the local structures. In a functional style of programming, I can use the same newVar like concept you would find in a constraint or inference solver, extend it with an equivalent newEq, to build a monadic environment to support my numerical code. The result is that "application code" just uses variable and defines equations, these are then assembled and only later does the code allocate just a few very large arrays to store the real data associated to the problem

This approach is "easy" because all invariants are built first and then used. Much modern applied math depends on hieararchical and adaptive grids. Hierarchy is probably easy to bring in to a functional style, adaptive grids maybe less obvious.