Decoupling Numerical Algorithms

Related Sites

I have a client who is asking about using OO design and C++ for an algorithmic intensive numerical application. The application implements various solving algorithms and formulas, and then they assemble them into scripts that solve complex calculations like air flow over wings and things like that.

They have an existing application written in FORTRAN. It has problems where changes propagate throughout the code and make the changes very expensive and error prone. They like the formulas they've written in FORTRAN because they are aggressively optimized by their FORTRAN compiler. [...]

It seems like the way they use the formulas as building blocks for the scripts lends itself to a component based approach for organizing the formulas. [...] But I am concerned about the performance issues. [...]

I have a similar client, who has a very similar problem. And they asked me a very similar question.

Their problem is twofold. First, they are breaking new ground with these kinds of algorithms, and are constantly having to revise them. Unfortunately the algorithmic code is highly coupled, depending upon huge arrays stuffed with lots of intermediate values. Often a minor change in the algorithm requires a major change in the software. Secondly, the algorithms sometimes fail in predictable ways, (e.g. floating point resolution problems, or infinities). If an algorithm fails in one way there is an alternative algorithm which might work. That alternative uses some of the intermediate values. But starting the new algorithm when the old one fails is tricky. The intermediate values no longer exist. So they have to start over from scratch and recreate them.

We decided to model the algorithm as a graph, (i.e. a set of nodes and edges). Every node represents a transformation of the data, every edge represents an intermediate state. Nodes can have several edges leaving them, some that represent failure states.

We create classes that represent the nodes. They have DO and UNDO functions. We create an algorithm manager that walks the graph of edges and nodes, calling the DO functions and watching for errors. If an error occurs, the managers calls UNDO walking backwards until the data is at a state where a different algorithm can be applied.

This structure does two things. It decouples the algorithmic code to some extent, and it also allows a strategy to manipulate the state of the intermediate values.

Numerical methods in fields such as computational fluid dynamics attempt to solve a system of differential equation, e.g. the Euler equations. The methods used to solve the system of equations vary based upon computational grid type (tetrahedral, prism, etc.), computational grid dimensions (2D, 3D), Node/Cell centered discretization schemes, etc. One of the goals in developing our system was to abstract the underlining solution method from the application. In other words, we wanted to plug and play solvers.

The approach I took was to abstract the solver from the application by using a sequencer. The sequencer provided an interface to the application for running a generic integration scheme. Each numerical method would require a unique sequencer. The sequencer would be responsible for "sequencing" the integration of the system of differential equations based upon the underlining method. This is analogous Robert's algorithm manager in that it "directs" the execution of events required to perform the integration.

Todd Veldhuizen, The Object-Oriented Numerics Page

John J. Barton and Lee R. Nackman, Scientific and Engineering C++