# What are the benefits and drawbacks inherent to using classes to encapsulate numerical algorithms?

Many algorithms used in scientific computing have a different inherent structure than algorithms commonly considered in less math-intensive forms of software engineering. In particular, individual mathematical algorithms tend to be highly complex, often involving hundreds or thousands of lines of code, yet nonetheless involve no state (i.e. are not acting upon a complex data structure) and can often be boiled down -- in terms of programmatic interface -- to a single function acting on an array (or two).

This suggests that a function, and not a class, is the natural interface to most algorithms encountered in scientific computing. Yet this argument offers little insight regarding how the implementation of complex, multi-part algorithms should be handled.

While the traditional approach has been to simply have one function that calls a number of other functions, passing the relevant arguments along the way, OOP offers a different approach, wherein algorithms can be encapsulated as classes. For clarity, by encapsulating an algorithm in a class, I mean creating a class wherein the algorithm inputs are entered into the class constructor, and then a public method is called to actually invoke the algorithm. Such an implementation of multigrid in C++ psuedocode might look like:

class multigrid {
private:
x_, b_
[grid structure]

restrict(...)
interpolate(...)
relax(...)
public:
multigrid(x,b) : x_(x), b_(b) { }
run()
}

multigrid::run() {
[call restrict, interpolate, relax, etc.]
}


My question is then as follows: what are the benefits and drawbacks of this kind of practice as compared to a more traditional approach without classes? Are there issues of extensibility or maintainability? To be clear, I am not intending to solicit opinion, but rather to better understand the downstream effects (i.e. those that might not arise until a codebase becomes quite large) of adopting such a coding practice.

• It's always a bad sign when your class name is an adjective rather than a noun. – David Ketcheson Jul 12 '12 at 6:26
• A class could serve as a stateless namespace for organizing functions in order to manage complexity, but there are other ways to manage complexity in languages that provide classes. (Namespaces in C++ and modules in Python come to mind.) – Geoff Oxberry Jul 12 '12 at 7:02
• @GeoffOxberry I can't speak to whether this is good or bad usage -- which is why I'm asking in the first place -- but classes, unlike namespaces or modules, can also manage "temporary state", e.g. the grid hierarchy in multigrid, that is discarded upon completion of the algorithm. – Ben Jul 12 '12 at 7:20

Having done numerical software for 15 years, I can unambiguously state the following:

• Encapsulation is important. You do not want to pass around pointers to data (as you suggest) since it exposes the data storage scheme. If you expose the storage scheme, you can never change it again because you will access the data all over the entire program. The only way to avoid this is to encapsulate the data into private member variables of a class and let only member functions act on it. If I read your question, you think of a function that computes the eigenvalues of a matrix as stateless, taking a pointer to the matrix entries as argument and returning the eigenvalues in some way. I think this is the wrong way to think about it. In my view, this function should be a "const" memberfunction of a class -- not because it changes the matrix, but because it is one that operates with the data.

• Most OO programming languages allow you to have private member functions. This is your way to break apart one large algorithm into a smaller one. For example, the various helper functions you need for the eigenvalue computation still operate on the matrix, and so would naturally be private member functions of a matrix class.

• Compared to many other software systems, it may be true that class hierarchies are often less important than, say, in graphical user interfaces. There are certainly places in numerical software where they are prominent -- Jed outlines one in another answer to this thread, namely the many ways one can represent a matrix (or, more generally, a linear operator on a finite dimensional vector space). PETSc does this very consistently, with virtual functions for all operations that act on matrices (they don't call it "virtual functions", but that's what it is). There are other areas in typical finite element codes where one uses this design principle of OO software. The ones that come to mind are the many kinds of quadrature formulas and the many kinds of finite elements, all of which are naturally represented as one interface/many implementations. Material law descriptions would also fall into this group. But it may be true that that's about it and that the rest of a finite element code doesn't use inheritance as pervasively as one may use it, say, in GUIs.

From only these three points, it should be clear that object oriented programming is most definitely applicable to numerical codes as well, and that it would be foolish to ignore the many benefits of this style. It may be true that BLAS/LAPACK don't use this paradigm (and that the usual interface exposed by MATLAB doesn't either) but I would venture the guess that every successful numerical software written in the past 10 years is, in fact, object oriented.

Encapsulation and data hiding are extremely important for extensible libraries in scientific computing. Consider matrices and linear solvers as two examples. A user just needs to know that an operator is linear, but it may have internal structure such as sparsity, a kernel, a hierarchical representation, a tensor product, or a Schur complement. In all cases, Krylov methods do not depend on the details of the operator, they only depend on the action of the MatMult function (and perhaps its adjoint). Similarly, the user of a linear solver interface (e.g. a nonlinear solver) only cares that the linear problem is solved, and should not need or want to specify the algorithm that is used. Indeed, specifying such things would impede the capability of the nonlinear solver (or other outer interface).

Interfaces are good. Depending on an implementation is bad. Whether you achieve this using C++ classes, C objects, Haskell typeclasses, or some other language feature is inconsequential. The capability, robustness, and extensibility of an interface is what matters in scientific libraries.

Classes should be used only if the structure of the code is hierarchical. Since you are mentioning Algorithms, their natural structure is a flow chart, not a hierarchy of objects.

In case of OpenFOAM, the algorithmic part is implemented in terms of generic operators (div, grad, curl, etc) that are basically abstract functions operating on different types of tensors, using different types of numerical schemes. This part of the code is basically built from a lot of generic algorithms operating on classes. This allows the client to write something like:

solve(ddt(U) + div(phi, U)  == rho*g + ...);


Hierarchies such as transport models, turbulence models, differencing schemes, gradient schemes, boundary conditions, etc are implemented in terms of C++ classes (again, generic on the tensor quantities).

I've noticed a similar structure in the CGAL library, where the various algorithms are packed together as groups of function objects bundled with geometrical information to form geometrical Kernels (classes), but this is again done to separate operations from geometry (point removal from a face, from a point data type).

Hierarchical structure ==> classes

Procedural, flow chart ==> algorithms

Even if this is an old question, I think that it is worth mentioning Julia's particular solution. What this language does is "class-less OOP": the main constructs are types, i.e., composite data objects akin to structs in C, on which an inheritance relation is defined. The types do not have "member functions", but each function has a type signature and accepts subtypes. For instance, you could have an abstract Matrix type and subtypes DenseMatrix, SparseMatrix, and have a generic method do_something(a::Matrix, b::Matrix) with specializationdo_something(a::SparseMatrix, b::SparseMatrix). Multiple dispatching is used to select the most appropriate version to call.

This approach is more powerful than class-based OOP, which is equivalent on dispatching based on inheritance on the first argument only, if you adopt the convention that "a method is a function with this as its first parameter" (common e.g. in Python). Some form of multiple dispatch can be emulated in, say, C++, but with considerable contortions.

The main distinction is that methods do not belong to classes, but they exist as separate entities and inheritance can happen on all parameters.

Some references:

http://docs.julialang.org/en/release-0.4/manual/methods/

http://assoc.tumblr.com/post/71454527084/cool-things-you-can-do-in-julia

https://thenewphalls.wordpress.com/2014/03/06/understanding-object-oriented-programming-in-julia-inheritance-part-2/

Two advantages of the OO approach could be:

• a long calculation from which you might or might not want different results. For example, if $\beta$ is the final output, but it depends on an intermediate result $\alpha$, you could have a calculate_alpha() method that caches the $\alpha$ result inside the instance. Then when you call calculate_beta(), it also calls calculate_alpha() if no $\alpha$ result has been cached yet.

• a calculation with several inputs, where if one input changes then the whole calculation doesn't necessarily have to be done over. For example, the calculate_f() method returns $f(x,y,z)$. If you then decide to redo the calculation for another value of $z$, you can call set_z() and the $z$ parameter is marked internally as 'dirty', so that when you call calculate_f() again, only the portion of the calculation that depends on $z$ is redone.