LinSearch[f_,x0_,p0_,maxstep_,steptol_,fx0_,gx0_] :=
Block[{x = x0, p = p0, alphatol = 10^(-4), RETCODE = -1, steplength, 
       linegradient, lambdamin, lambda, lambdatmp,
       fx, fx2, temp1, temp2, lambdadiff, discr, a, b},
  steplength = Sqrt[Apply[Plus, p^2]];
  If[steplength > maxstep,
    p = p*maxstep/steplength; steplength = maxstep];
  linegradient = Apply[Plus, gx0*p];
  lambdamin = steptol/Max[Abs[p]/Max[Abs[x],1]];
  lambda = 1.0;
  While[RETCODE < 0,
    x = x0 + lambda*p; fx = Apply[f, x];
    If[fx <= (fx0 + alphatol * lambda * linegradient),
      (* sufficiently good solution found *)
      RETCODE = 0,
      If[lambda < lambdamin,  (* didn't find a minimum distinct from x0 *)
        x = x0; RETCODE = 1,
        If[lambda == 1.0,
          (* calculate new lambda for the first backtracking *)
          lambdatmp = -linegradient/(2*(fx - fx0 - linegradient)),
          (* calculate new lambda for later backtrackings *)
          temp1 = fx - fx0 - lambda*linegradient;
          temp2 = fx2 - fx0 - lambda2*linegradient;
          lambdadiff = 1/(lambda-lambda2);
          a = lambdadiff*(temp1/lambda^2-temp2/lambda2^2);
          b = lambdadiff*(lambda*temp2/lambda2^2-lambda2*temp1/lambda^2);
          discr = b^2 - 3.0*a*linegradient;
          If[a == 0.0,	
            lambdatmp = -linegradient/(2.0*b),
            lambdatmp = (Sqrt[discr] - b)/(3.0*a)];
          If[lambdatmp > 0.5*lambda, lambdatmp = 0.5*lambda]
        ] (* end of If *)
      ] (* end of If *)
    ]; (* end of If *)
    If[RETCODE < 0, fx2 = fx; lambda2 = lambda;
      lambda = Max[lambdatmp, 0.1*lambda]]
  ];  (* end of While *)
  {x, fx, RETCODE}  (* return values *)
] (* end of Block *)
