stepit <- function(linemat, ips, dblwrd, MAXSTEP)
{
#STEPIT computes next step size in line search algorithm
#  Arguments:
#  LINEMAT:  Row 1 contains step values
#            Row 2 contains slope values
#            Row 3 contains function values
#  IPS:      If 1, previous slope was positive
#  DBLWRD:   Vector of length 2:  dblwrd[1] TRUE means step halved
#                                 dblwrd[2] TRUE means step doubled
#  MAXSTEP:  maximum size of step

#  Last modified 3 January 2008

test1 <- abs(linemat[2,5]) < abs(linemat[2,1])/10
test2 <- linemat[3,5] > linemat[3,1]
test3 <- linemat[2,5] > 0
if ((test1 || !test3) && test2) {
   #  ************************************************************
   #  function is worse and either slope is satisfory or negative
   ips <- 0        #  step is halved
   if (dblwrd[2]) {
      ind = -5
      return(list(linemat, ips, ind, dblwrd))
   }
   linemat[1,5] <- min(c(linemat[1,5]/2, MAXSTEP))
   linemat[,2] <- linemat[,1]
   linemat[,3] <- linemat[,1]
   dblwrd <- c(1, 0)
   ind <- 2
   return(list(linemat, ips, ind, dblwrd))
}
#  *********************************************************
if (test1) {
   #  test1 means successful convergence
   ind <- 0
   return(list(linemat, ips, ind, dblwrd))
}
#  **********************************************************
if (test3) {
   #  Current slope is positive
   ips <- 1
   linemat[,4] <- linemat[,5]
   deltaf <- linemat[3,3] - linemat[3,5]
   z <- (3/(linemat[1,5] - linemat[1,3])) * deltaf + linemat[2,3] + linemat[2,5]
   w <- z * z - linemat[2,3] * linemat[2,5]
   if (abs(linemat[2,3] + linemat[2,5] + 2 * z) >= 1e-05 & w > 0) {
     w <- sqrt(w)
     linemat[1,5] <- linemat[1,3] + (1 - ((linemat[2,5] + w - z)/(linemat[2,5] - linemat[2,3] +
         2 * w))) * (linemat[1,5] - linemat[1,3])
   } else {
           #  linear interpolation necessary
           aerror <- linemat[1,3]
           if (linemat[1,5] > linemat[1,3]) aerror <- linemat[1,5]
           linemat[1,5] <- linemat[1,3] - linemat[2,3] *((linemat[1,5] - linemat[1,3])/(linemat[2,5] - linemat[2,3]))
           if (linemat[1,5] > 2 * aerror) linemat[1,5] <- 2 * aerror
   }
   linemat[1,5] <- min(c(linemat[1,5], MAXSTEP))
   dblwrd <- c(0,0)
   ind <- 2
   return(list(linemat, ips, ind, dblwrd))
}
#  *************************************************************
#  Current slope is negative or zero
linemat[,2] <- linemat[,3]
linemat[,3] <- linemat[,5]
if (ips == 1) {
   #  *****************************************************
   #  previous slope is positive
   deltaf <- linemat[3,5] - linemat[3,4]
   z <- (3/(linemat[1,4] - linemat[1,5])) * deltaf + linemat[2,5] + linemat[2,4]
   w <- z * z - linemat[2,5] * linemat[2,4]
   if (abs(linemat[2,5] + linemat[2,4] + 2 * z) >= 1e-05 && w > 0) {
     w <- sqrt(w)
     linemat[1,5] <- linemat[1,5] + (1 - ((linemat[2,4] + w - z)/(linemat[2,4] - linemat[2,5] + 2 * w))) * (linemat[1,4] - linemat[1,5])
   } else {
           aerror <- linemat[1,5]
           if (linemat[1,4] > linemat[1,5]) aerror <- linemat[1,4]
           linemat[1,5] <- linemat[1,5] - linemat[2,5] *((linemat[1,4] - linemat[1,5])/(linemat[2,4] - linemat[2,5]))
           if (linemat[1,5] > 2 * aerror) linemat[1,5] <- 2 * aerror
   }
   linemat[1,5] <- min(c(linemat[1,5], MAXSTEP))
   dblwrd <- c(0,0)
   ind <- 2
   return(list(linemat, ips, ind, dblwrd))
}
#  ******************************************************
if ((linemat[2,3] - linemat[2,2]) * (linemat[1,3] - linemat[1,2]) > 0) {
   #  previous slope is negative
   z <- (3/(linemat[1,3] - linemat[1,2])) * (linemat[3,2] - linemat[3,3]) + linemat[2,2] + linemat[2,3]
   w <- z * z - linemat[2,2] * linemat[2,3]
   if (abs(linemat[2,2] + linemat[2,3] + 2 * z) >= 1e-05 & w > 0) {
           w <- sqrt(w)
           linemat[1,5] <- linemat[1,2] + (1 - ((linemat[2,3] + w - z)/(linemat[2,3] - linemat[2,2] + 2 * w))) * (linemat[1,3] - linemat[1,2])
   } else {
     linemat[1,5] <- linemat[1,2] - linemat[2,2] * ((linemat[1,3] - linemat[1,2])/(linemat[2,3] - linemat[2,2]))
   }
   linemat[1,5] <- min(c(linemat[1,5], MAXSTEP))
   dblwrd <- c(0,0)
   ind <- 2
   return(list(linemat, ips, ind, dblwrd))
} else {
   #  previous slope also negative but not as much
   if (dblwrd[1]) {
           ind <- 5
           return(list(linemat, ips, ind, dblwrd))
   } else {
           linemat[1,5] <- 2 * linemat[1,5]
           linemat[1,5] <- min(c(linemat[1,5], MAXSTEP))
           dblwrd <- c(0,1)
           ind <- 2
           return(list(linemat, ips, ind, dblwrd))
   }
}
ind <- 2
return(list(linemat, ips, ind, dblwrd))
}