Nearest-neighbor interpolation

Deformations begin from the task of selecting a value val from an array vec(ix), ix=1,nx. The points of the array are at locations x = x0+dx*(ix-1). Given the location x of the desired value we backsolve for ix. In Fortran, conversion of a real value to an integer is done by truncating the fractional part of the real value. To get rounding up as well as down, we add a half before conversion to an integer, namely ix=int(1.5+(x-x0)/dx). This gives the nearest neighbor. The adjoint to extracting a value from a vector is putting it back. A convenient subroutine for nearest-neighbor interpolation is spot0().  

# Nearest neighbor interpolation, essentially:   val = vec( 1.5 + (t-t0)/dt)
#
subroutine spot0( adj, add, nt,t0,dt, t, val,    vec    )
integer it,       adj, add, nt
real                           t0,dt, t, val,    vec( nt)
call adjnull(     adj, add,              val, 1, vec, nt)
it = 1.5 + (t-t0) / dt
if( 0 < it && it <= nt) 
        if( adj == 0 )                          # add value onto vector
                vec( it) = vec( it) + val
        else                                    # take value from vector
                val      = val      + vec( it)
return; end

Recall subroutine advance() . For jump==0 its matrix equivalent is an identity matrix. For other values of jump, the identity matrix has its diagonal shifted up or down. Now examine subroutine spot0() and think about its matrix equivalent. Since its input is a single value and its output is a vector, that means its matrix is a column vector so the adjoint operator is a row vector. The vector is all zeros except for somewhere where there is a ``1''.