This week, the third post in a series on the new finite-difference framework. It started here, in case you need to catch up.

Mandatory plug: there are still places availables for my Introduction to QuantLib Development course. For details, click on the link in the banner above.

Follow me on Twitter if you want to be notified of new posts, or add me to your Google+ circles, or subscribe via RSS: the buttons for that are in the footer. Also, make sure to check my Training page.

Operators, iterators and layouts

The interface of the FdmLinearOp class is shown in the listing below, and its main feature is an apply method that returns the result of applying the operator to a given array; you might remember a similar applyTo method defined by the old-style operators. (I’ll gloss over the toMatrix method, which returns some kind of representation of the operator as a matrix and seems to be only used for inspection.)

class FdmLinearOp {
    typedef Array array_type;
    virtual ~FdmLinearOp() { }
    virtual Disposable<array_type> apply(
                           const array_type& r) const = 0;
    virtual Disposable<SparseMatrix> toMatrix() const = 0;

However, you probably noticed that the array_type used in the interface is defined as a typedef to Array, which implements a one-dimensional array. What of my claim that we can also cover multi-dimensional operators, then? Well, the framework maps the elements of a multi-dimensional array to the elements of a one-dimensional one by flattening the former and turning it into the latter.

As an example, let’s say that we have a three-dimensional array \( a_{i,j,k} \) with \( 0 \le i < M \), \( 0 \le j < N \), and \( 0 \le k < P \). We’ll map it into a one-dimensional sequence via a simple lexicographic ordering of the elements. We start from \( a_{0,0,0} \) as the first element, followed by \( a_{1,0,0} \) up to \( a_{M-1,0,0} \). Then we increase the second index and reset the first, obtaining \( a_{0,1,0} \) and following with \( a_{1,1,0} \) to \( a_{M-1,1,0} \). The next element is \( a_{0,2,0} \), and this repeats until we get to \( a_{M-1,N-1,0} \), at which point we follow with \( a_{0,0,1} \). You got the gist; we go through the whole thing again to \( a_{M-1,N-1,1} \), then \( a_{0,0,2} \) and so on until we finish with \( a_{M-1,N-1,P-1} \).

The FdmLinearOpIterator and FdmLinearOpLayout classes, whose interfaces are shown below, cooperate to implement the logic above. And yes, even though both classes are named after the FdmLinearOp class, they actually work on the arrays to which operators are applied.

class FdmLinearOpIterator {
    explicit FdmLinearOpIterator(const std::vector<Size>& dim);
    explicit FdmLinearOpIterator(Size index = 0);
    FdmLinearOpIterator(const std::vector<Size>& dim,
                        const std::vector<Size>& coordinates,
                        Size index)

    void operator++();
    bool operator!=(const FdmLinearOpIterator&);

class FdmLinearOpLayout {
    FdmLinearOpLayout(const std::vector<Size>& dim);

    FdmLinearOpIterator begin() const;
    FdmLinearOpIterator end() const;

    Size index(const std::vector<Size>& coordinates) const;

    Size neighbourhood(const FdmLinearOpIterator& iterator,
                       Size i, Integer offset) const;
    Size neighbourhood(const FdmLinearOpIterator& iterator,
                       Size i1, Integer offset1,
                       Size i2, Integer offset2) const;

Instances of the FdmLinearOpIterator keep track of a tuple of indices into the several dimensions of an array (like the \( (i,j,k) \) tuples above) and of the corresponding single index into the flattened array; to do that, they must also store the sizes of the array along each dimension (the ones that I called \( M \), \( N \) and \( P \) in the example).

Depending on the constructor you’ll call, you can pass explicitly all the state or just part of it. The constructor taking only a vector of dimensions stores the sizes of the array and implicitly initializes all the indices to 0; that is, it returns the iterator pointing to the beginning of the array. The constructor taking a single index builds a kind of sentinel value; it cannot be incremented, because it doesn’t know the size of the array, but when it is passed the total number of elements of the array it returns an end iterator to which others can be compared to end a loop, as in:

std::vector<Size> dims(...); // { M, N, P }
FdmLinearOpIterator begin(dims);
FdmLinearOpIterator end(M*N*P);
for (FdmLinearOpIterator i=begin; i!=end; ++i)
    ... // do something

Finally, the third constructor takes the dimensions and the whole set of indices, and thus allows you to create an iterator pointing to any element of the array. For performance reasons (I guess) it leaves to the caller to check that the parameters are consistent, that is, that the set of indices into the array actually corresponds to the single index into the flattened array.

A FdmLinearOpIterator instance is not really an iterator in the C++ sense, since it doesn’t implement the whole required interface; in particular, there’s no dereference operator returning the value it points to. What it does implement is operator++, which increments the iterator, and operator!=, which compares two iterators; the two methods are shown in the next listing. (Providing only operator!= and not operator== is enough to support the i != end idiom used in for loops, but might be a tad extreme.)

void FdmLinearOpIterator::operator++() {
    for (Size i=0; i < dim_.size(); ++i) {
        if (++coordinates_[i] == dim_[i])
            coordinates_[i] = 0;

bool FdmLinearOpIterator::operator!=(
                    const FdmLinearOpIterator& iterator) {
    return index_ != iterator.index_;

To round up their interface, iterators also declare the coordinates inspector, which returns the set of indices into the array, and index, returning the corresponding single index.

The increment operator literally implements the iterating strategy I described in the previous page. It starts by incrementing the first index; if doing so causes it to reach the size of the first dimension, it resets the index to 0 and repeats both the increment and check on the next dimension; otherwise, it breaks out of the loop. Thus, for instance, calling operator++ with the indices set to \( (M-1,N-1,2) \) will increment the first and get to \( (M,N-1,2) \), then reset it and increment the second to \( (0,N,2) \), then reset the second too and increment the third to \( (0,0,3) \), then finally exit. Calling it with the indices set to \( (2,0,0) \), instead, will just increment the first to \( (3,0,0) \) and exit. In any case, the single index into the flattened array is also incremented so that it stays in sync with the other indices.

Finally, operator!= just compares two iterators by comparing their two indices into the flattened array (also known as the indices for which I should probably get a shorter name at this point). This allows iterators to work as sentinels even if they don’t have a full set of indices and can’t be incremented.

The FdmLinearOpLayout class (whose interface we have seen in a previous listing, and whose implementation is shown in the next listing) provides some utilities for working with iterators based on a given array layout, or set of dimensions.

                                const std::vector<Size>& dim)
: dim_(dim), spacing_(dim.size()) {
    spacing_[0] = 1;
    std::partial_sum(dim.begin(), dim.end()-1,
        spacing_.begin()+1, std::multiplies<Size>());
    size_ = spacing_.back()*dim.back();

FdmLinearOpIterator FdmLinearOpLayout::begin() const {
    return FdmLinearOpIterator(dim_);

FdmLinearOpIterator FdmLinearOpLayout::end() const {
    return FdmLinearOpIterator(size_);

Size FdmLinearOpLayout::index(
                  const std::vector<Size>& coordinates) const {
    return std::inner_product(coordinates.begin(),
                              spacing_.begin(), Size(0));

Most of its calculations are based on what it calls spacing between the elements of the array along the various directions. I’ll use the same example as before, and consider a three-dimensional array with sizes \( (M,N,P) \) along the three dimensions. If we take two elements \( (i,j,k) \) and \( (i+1,j,k) \), they occupy consecutive places on the one-dimensional flattened array (remember that we increment the first index first while enumerating the elements) so we say that the spacing between them is \( 1 \). If we take two elements \( (i,j,k) \) and \( (i,j+1,k) \), their positions in the flattened array are separated by a whole turn of the first index: enumerating the elements between them will cause you to go through \( (i+1,j,k) \), \( (i+2,j,k) \) and so on up to \( (M-1,j,k) \), then to \( (0,j+1,k) \), and then again through \( (1,j+1,k) \), \( (2,j+1,k) \) up to \( (i,j+1,k) \). A bit of thinking will convince you that the spacing between the two elements is \( M \), and a similar argument will find a spacing of \( M \times N \) between elements \( (i,j,k) \) and \( (i,j,k+1) \).

The constructor of a FdmLinearOpLayout instance takes a general set of sizes \( (n_1, \ldots, n_N) \) and calculates the spacings \( (1, n_1, n_1 \times n_2, \ldots, \prod_1^{N-1} n_i) \) as well as the total size \( \prod_1^N n_i \); the calculations is written in STL speak, but it shouldn’t be difficult to map it to the formulas.

Given the sizes and the spacings, the class can implement a number of methods that let you (and library code) work with iterators. The begin and end methods return the corresponding iterators, so that once you have a layout you no longer need to care about sentinels and iterator constructors: loops can be written as:

std::vector<Size> dims(...);
FdmLinearOpLayout l(dims);
for (FdmLinearOpIterator i=l.begin(); i!=l.end(); ++i)
    ... // do something

The index method provides a conversion between a set of indices and the single index into the flattened array. Other methods, such as neighborhood, work at a slightly higher level; they return the index (or the iterator) of the element you get by starting from a given element, whose position is described by an iterator, and moving a given number of steps in a given dimension. An overload that moves in two directions is also available for convenience, even though it’s equivalent to two consecutive calls to the first version.

At this point, we can finally make sense of the interface of the FdmMesher class shown in the previous post. The methods dplus, dminus and location all take as first argument an iterator to specify a point in the mesh; the second argument specifies the direction along which we’re moving. (The points in the mesh are in one-to-one correspondence with the elements of the array containing the values of the instrument being priced on the mesh, thus an iterator can denote either one.) Finally, the layout method returns an instance of FdmLinearOpLayout that can be used to create iterators.

In next post: a few actual operators.