#include <iostream>
/**************************************************************/
/* This program demonstrates how to avoid global variables    */
/* in C++ programs when using Fortran subroutines. An example */
/* is from minimization of a smooth function                  */
/**************************************************************/
// External Fortran subroutine
extern "C"
void uncmin_(const int* N, float* X0, void (*FCN)(const int* N, const float* X, float* F),
	     float* X, float* F, int* INFO, float* W, int* LW);

// First class 
class Parabola{
  float u0, u1;
  static Parabola* pt; // static pointer which allows static member function FCN to access class data
public:
  Parabola(float u0_, float u1_);
  double Value(double x0, double x1)
    { return (x0-u0)*(x0-u0)+(x1-u1)*(x1-u1);}
  static void FCN(const int* N, const float* X, float* F);// Static member function to be used with Fortran sub.
};
Parabola* Parabola::pt;

Parabola::Parabola(float u0_, float u1_) : u0(u0_), u1(u1_)
{Parabola::pt = this;}

void Parabola::FCN(const int* N, const float* X, float* F)
{
  if (*N!=2) std::cerr<<"Not right number of variables!"<<std::endl;
  *F = pt->Value(X[0],X[1]);
}

// Second class 
class Parabola1{
  float u0, u1, u2;
  static Parabola1* pt; // static pointer which allows static member function FCN to access class data
public:
  Parabola1(float u0_, float u1_, float u2_);
  double Value(double x0, double x1, double x2)
    { return (x0-u0)*(x0-u0)+(x1-u1)*(x1-u1)+(x2-u2)*(x2-u2);}
  static void FCN(const int* N, const float* X, float* F);// Static member function to be used with Fortran sub.
};
Parabola1* Parabola1::pt;

Parabola1::Parabola1(float u0_, float u1_, float u2_) : u0(u0_), u1(u1_), u2(u2_)
{Parabola1::pt = this;}

void Parabola1::FCN(const int* N, const float* X, float* F)
{
  if (*N!=3) std::cerr<<"Not right number of variables!"<<std::endl;
  *F = pt->Value(X[0],X[1],X[2]);
}

// C++ wrapper for the Fortran subroutine. It works with any function
// or class, therefore we used template functor.
template <class functor>
float Minimize(int N, float x[], functor& Fun)
{
  float* x0 = new float[N];
  int lw = N*(N+10);
  float* w = new float[lw];
  for (int i=0; i<N; i++) x0[i]=x[i];
  float f; int info;
  uncmin_(&N,x0,Fun,x,&f,&info,w,&lw);
  delete[] w;
  delete[] x0;
  return f;
}

int main()
{
  float X0[2]={1,1};
  float F;
  
  Parabola parab(4.,5.);
  
  F = Minimize(2, X0, Parabola::FCN);
  
  std::cout<<X0[0]<<" "<<X0[1]<<" "<<F<<std::endl;

  Parabola1 param1(4.,5.,6.);
  
  F = Minimize(3, X0, Parabola1::FCN);
  
  std::cout<<X0[0]<<" "<<X0[1]<<" "<<X0[2]<<" "<<F<<std::endl;
  return 0;
}