Imports MathNet.Numerics.LinearAlgebra.Double

Public Class NDSpline2

Private BasisCache(,,)() As DenseMatrix

Public Structure SplineInterpResult
    Public z() As Double
    Public dzdx(,) As Double
End Structure

....a bunch of private variables declared as dynamic arrays at class level,
then a "New" subroutine that appropriately dimensions all these arrays.
Perhaps most relevant is the call within the New subroutine to the below
routine which populates a cache matrices which hopefully will make the
interpolate function down below run more quickly.  The concept in all this is
that interpolate is called lots of times and needs to be as fast as possible.

Private Sub BasisCacher()

    ReDim BasisCache(2, _Dimensions - 1, 1)
    For k = 0 To 2
        For d = 0 To _Dimensions - 1
            ReDim BasisCache(k, d, 0)(_NumberOfKnots(d))
            ReDim BasisCache(k, d, 1)(_NumberOfKnots(d))
            For LIminus3 = 0 To _NumberOfKnots(d)
                BasisCache(k, d, 0)(LIminus3) = New DenseMatrix(k + 1, k + 2)
                BasisCache(k, d, 1)(LIminus3) = New DenseMatrix(k + 1, k + 2)
                For j = 0 To k
                    Index1Holder = LIminus3 + 4 + j
                    Index2Holder = LIminus3 + 3 + j - k
                    BasisCache(k, d, 0)(LIminus3)(j, j) =
                        _KnotVector(d) (Index1Holder) /
                        (_KnotVector(d)(Index1Holder) -
                        _KnotVector(d)(Index2Holder))
                    BasisCache(k, d, 0)(LIminus3)(j, j + 1) =
                        -_KnotVector(d)(Index2Holder) /
                        (_KnotVector(d)(Index1Holder) -
                        _KnotVector(d)(Index2Holder))
                    BasisCache(k, d, 1)(LIminus3)(j, j) = -1 /
                        (_KnotVector(d)(Index1Holder) -
                        _KnotVector(d)(Index2Holder))
                    BasisCache(k, d, 1)(LIminus3)(j, j + 1) = 1 /
                        (_KnotVector(d)(Index1Holder) -
                        _KnotVector(d)(Index2Holder))
                Next j
            Next LIminus3
        Next d
    Next k

End Sub

Function Interpolate_NonUniformNonRecurringKnots(x() As Double) As
SplineInterpResult

 'Build Basis Functions with cached values
    For d = 0 To _Dimensions - 1
        LowerIndex(d) = Array.BinarySearch(_KnotVectorStorage(d), x(d))
        If LowerIndex(d) < 0 Then LowerIndex(d) = -LowerIndex(d) - 2
        For k = 0 To 2 'degree=k+1
            BasisHolder(k) = BasisCache(k, d, 0)
            (LowerIndex(d)).Add(BasisCache(k, d, 1)
            (LowerIndex(d)).Multiply(x(d)))
        Next k
        BasisHolderAllDims(d) = (BasisHolder(0) * BasisHolder(1) *
                                BasisHolder(2)).Row(0)
    Next d

    ....stuff that I haven't built yet which will calculate z and dzdx.
    The lines above about BasisHolder(k) and BasisHolderAllDims(d) are
    exemplary of where the exception occurs.

    Return MySplineResult

End Function