1 files changed, 82 insertions, 0 deletions
diff --git a/methods/linmin/mnbrak.hpp b/methods/linmin/mnbrak.hpp
new file mode 100644
index 0000000..9871473
--- /dev/null
+++ b/methods/linmin/mnbrak.hpp
@@ -0,0 +1,82 @@
+#ifndef MNBRAK_HPP
+#define MNBRAK_HPP 
+#define OPT_HEADER
+//#include "optimizer.hpp"
+#include "bas_util.hpp"
+namespace opt_utilities
+{
+
+  
+  template <typename T>
+  void mnbrak(T& ax,T& bx,T& cx,T& fa,T& fb,T& fc,func_obj<T,T>& func)
+  {
+    const T GOLD=1.618034;
+    const T GLIMIT=100;
+    const T TINY=std::numeric_limits<T>::epsilon();
+    T ulim,u,r,q,fu;
+    fa=func.eval(ax);
+    fb=func.eval(bx);
+    
+    if(fb>fa)
+      {
+	//shft(dum,ax,bx,dum);
+	//shft(dum,fb,fa,dum);
+	std::swap(ax,bx);
+	std::swap(fa,fb);
+      }
+
+    cx=bx+GOLD*(bx-ax);
+    fc=func.eval(cx);
+    while(fb>fc)
+      {
+	r=(bx-ax)*(fb-fc);
+	q=(bx-cx)*(fb-fa);
+	u=bx-T((bx-cx)*q-(bx-ax)*r)/
+	  T(T(2.)*sign(T(max(T(tabs(T(q-r))),T(TINY))),T(q-r)));
+	ulim=bx+GLIMIT*(cx-bx);
+	if((bx-u)*(u-cx)>0.)
+	  {
+	    fu=func.eval(u);
+	    if(fu<fc)
+	      {
+		ax=bx;
+		bx=u;
+		fa=fb;
+		fb=fu;
+		return;
+	      }
+	    else if(fu>fb)
+	      {
+		cx=u;
+		fc=fu;
+		return;
+	      }
+	    u=cx+GOLD*(cx-bx);
+	    fu=func.eval(u);
+	  }
+	else if((cx-u)*(u-ulim)>0.)
+	  {
+	    fu=func.eval(u);
+	    if(fu<fc)
+	      {
+		shft3(bx,cx,u,T(cx+GOLD*(cx-bx)));
+		shft3(fb,fc,fu,func.eval(u));
+	      }
+	  }
+	else if((u-ulim)*(ulim-cx)>=0)
+	  {
+	    u=ulim;
+	    fu=func.eval(u);
+	  }
+	else
+	  {
+	    u=cx+GOLD*(cx-bx);
+	    fu=func.eval(u);
+	  }
+	shft3(ax,bx,cx,u);
+	shft3(fa,fb,fc,fu);
+      }
+  }
+}
+
+#endif