> Gradmat
function(parvec, infcn, eps = 1e-06)
{
# Function to calculate the difference-quotient approx gradient
#   (matrix) of an arbitrary input (vector) function    infcn
# Now recoded to use central differences !
        dd = length(parvec)
        aa = length(infcn(parvec))
        epsmat = (diag(dd) * eps)/2
        gmat = array(0, dim = c(aa, dd))
        for(i in 1:dd)
                gmat[, i] <- (infcn(parvec + epsmat[, i]) - 
                    infcn(parvec - epsmat[, i]))/eps
        if(aa > 1)  gmat   else c(gmat)
}



NRroot
function(inipar, infcn, nmax = 25, stoptol = 1e-05, 
        eps = 1e-06, gradfunc = NULL)
{
    if(is.null(gradfunc))
        gradfunc = function(x) Gradmat(x, infcn, eps)
    ctr = 0
    newpar = inipar
    oldpar = inipar - 1
    while(ctr < nmax & sqrt(sum((newpar - oldpar)^2)) > stoptol) {
        oldpar <- newpar
        newpar <- oldpar - solve(gradfunc(oldpar), c(infcn(oldpar)))
        ctr <- ctr + 1
    }
    list(nstep = ctr, initial = inipar, final = newpar, 
        funcval = infcn(newpar))
}




> GradSrch
function(inipar, infcn, step, nmax = 25, stoptol = 1e-05,
        unitfac = F, eps = 1e-06, gradfunc = NULL)
{
# Function to implement Steepest-descent. The   unitfac
#  condition indicates whether or not the supplied step-length
#  factor(s) multiply the negative gradient itself, or
#  the unit vector in the same direction.
   if(is.null(gradfunc))
         gradfunc = function(x) Gradmat(x, infcn, eps)
   steps = if(length(step) > 1) step else rep(step, nmax)
   newpar = inipar
   oldpar = newpar - 1
   ctr = 0
   while(ctr < nmax & sqrt(sum((newpar - oldpar)^2)) > stoptol) {
         ctr = ctr + 1
         oldpar = newpar
         newstep = gradfunc(oldpar)
         newstep = if(unitfac) newstep/sum(newstep^2) else newstep
         newpar = oldpar - steps[ctr] * newstep
   }
   list(nstep = ctr, initial = inipar, final = newpar,
         funcval = infcn(newpar))
}