(require "f2cl_macros")
(require "linmin")

(defun powell (p xi n np ftol iter fret &key (nmax 20) (itmax 200))
 (declare (type (simple-array double-float (*)) p))
 (declare (type (simple-array double-float (* *)) xi))
 (declare (type fixnum n)) (declare (type fixnum np))
 (declare (type double-float ftol)) (declare (type fixnum iter))
 (declare (type double-float fret)) (declare (type fixnum nmax))
 (declare (type fixnum itmax))
 (prog
  ((pt (make-array '(nmax) :element-type 'double-float))
   (ptt (make-array '(nmax) :element-type 'double-float))
   (xit (make-array '(nmax) :element-type 'double-float)) (t_ 0.0d0)
   (fptt 0.0d0) (i 0) (del 0.0d0) (ibig 0) (fp 0.0d0) (j 0)
  )
  (declare (type (simple-array double-float (*)) pt))
  (declare (type (simple-array double-float (*)) ptt))
  (declare (type (simple-array double-float (*)) xit))
  (declare (type double-float t_)) (declare (type double-float fptt))
  (declare (type fixnum i)) (declare (type double-float del))
  (declare (type fixnum ibig)) (declare (type double-float fp))
  (declare (type fixnum j)) (setf fret (func p))
  (fdo ((j 1 (+ j 1))) ((> j n) nil) (tagbody (fset (fref pt j) (fref p j))))
  (setf iter 0) label1 (setf iter (+ iter 1)) (setf fp fret) (setf ibig 0)
  (setf del 0.0)
  (fdo ((i 1 (+ i 1))) ((> i n) nil)
   (tagbody
    (fdo ((j 1 (+ j 1))) ((> j n) nil)
     (tagbody (fset (fref xit j) (fref xi j i)))
    )
    (multiple-value-setq (p xit n fret) (linmin p xit n fret))
    (cond
     ((> (abs (+ fp (- fret))) del) (setf del (abs (+ fp (- fret))))
      (setf ibig i)
  ))))
  (if (<= (* 2.0 (abs (+ fp (- fret)))) (* ftol (+ (abs fp) (abs fret))))
   (go end_label)
  )
  (if (= iter itmax) (error "Powell exceeding maximum iterations."))
  (fdo ((j 1 (+ j 1))) ((> j n) nil)
   (tagbody (fset (fref ptt j) (+ (* 2.0 (fref p j)) (- (fref pt j))))
    (fset (fref xit j) (+ (fref p j) (- (fref pt j))))
    (fset (fref pt j) (fref p j))
  ))
  (setf fptt (func ptt)) (if (>= fptt fp) (go label1))
  (setf t_
   (+
    (* (* 2.0 (+ (+ fp (* (* -1 2.0) fret)) fptt))
     (expt (+ (+ fp (- fret)) (- del)) 2)
    )
    (* (* -1 del) (expt (+ fp (- fptt)) 2))
  ))
  (if (>= t_ 0.0) (go label1))
  (multiple-value-setq (p xit n fret) (linmin p xit n fret))
  (fdo ((j 1 (+ j 1))) ((> j n) nil)
   (tagbody (fset (fref xi j ibig) (fref xit j)))
  )
       end_label
  (return (values p xi n np ftol iter fret))
))