At this point we have only implemented one solver, SEP's traditional workhorse, the conjugate gradient solverClaerbout (1994).
Creating an effective interface for this subprogram presented a challenge because model space and data space may each be of any dimensionality, such as a model which is sep_1d and a residual which is a sep_2d; and both spaces are often made up of several parts, which may themselves be of different dimensionality. A problem may have, for instance, a data residual which is a vector and a model residual which is a plane. This presents far too many possibilities to consider an overloaded solver. Such a thing is at any rate needless, since the conjugate gradient algorithm needs only to take dot products; the Fortran77 incarnation of the solver assumes that whatever matrices it is given are vectors.
To handle diversely-dimensioned, compound spaces, we decided on a similar tactic. Anything being passed to cg should be in the form of a 2D data set, a vector of pointers to traces. This is accomplished by collapsing larger dimensionality data sets via the function vectorize. Several distinct spaces (a model residual and a data residual) may be vectorized and then joined together via the function merge. A call to subroutine cg would then look like this:
call cg(iter, niter, merge(vectorize(model),vectorize(filter)), & merge(vectorize(dmodel),vectorize(dfilter)), & merge(vectorize(modelresidual),vectorize(dataresidual)), & merge(vectorize(dmodelresidual),vectorize(ddataresidual)))
In fact, a linear inverse problem could be entirely contained within the call to the solver, since procedures may be passed as arguments.