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
`vectorize`d 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.

11/12/1997