LinSearch[f_,x0_,p0_,maxstep_,steptol_,fx0_,gx0_] :=
Block[{x = x0, p = p0, alphatol = 10^(-4), RCODE = -1, 
       steplen, linegrad, fx, lambdamin, lambda, lmb2, 
       lmbtmp, fx2, temp1, temp2, lmbdiff, discr, a, b},
  steplen = Sqrt[Apply[Plus, p^2]];
  If[steplen > maxstep,
    p = p*maxstep/steplen; steplen = maxstep];
  linegrad = Apply[Plus, gx0*p];
  lambdamin = steptol/Max[Abs[p]/Max[Abs[x],1]];
  lambda = 1.0;
  While[RCODE < 0,
    x = x0 + lambda*p; fx = Apply[f, x];
    If[fx <= (fx0 + alphatol * lambda * linegrad),
      (* tarpeeksi hyvä ratkaisu löytyi *)
      RCODE = 0,
      If[lambda < lambdamin,  
        (* ei löytynyt x0:sta poikkeavaa ratkaisua *)
        x = x0; RCODE = 1,
        If[lambda == 1.0,
          (* laske uusi lambda ensimmäisessä vaiheessa *)
          lmbtmp = -linegrad/
                      (2*(fx - fx0 - linegrad)),
          (* laske uusi lambda myöh. käyttöä varten *)
          temp1 = fx - fx0 - lambda*linegrad;
          temp2 = fx2 - fx0 - lmb2*linegrad;
          lmbdiff = 1/(lambda-lmb2);
          a = lmbdiff*(temp1/lambda^2-temp2/lmb2^2);
          b = lmbdiff*(lambda*temp2/lmb2^2-
              lmb2*temp1/lambda^2);
          discr = b^2 - 3.0*a*linegrad;
          If[a == 0.0,	
            lmbtmp = -linegrad/(2.0*b),
            lmbtmp = (Sqrt[discr] - b)/(3.0*a)];
          If[lmbtmp > 0.5*lambda, lmbtmp = 0.5*lambda]
        ]
      ] 
    ]; 
    If[RCODE < 0, fx2 = fx; lmb2 = lambda;
      lambda = Max[lmbtmp, 0.1*lambda]]
  ];  
  {x, fx, RCODE}  (* palauta arvot *)
] 
