#include "libscl.h"

using namespace std;
using namespace scl;

class nleqns : public nleqns_base {
public:
  bool get_f(const realmat& x, realmat& f) 
  {
    if (x.nrow() != 2) error("Error, nleqns, wrong dimension for x");
    if (f.nrow() != 2)f.resize(2,1);
    f[1] = 2.0 - exp(1.0*x[1] + 1.0*x[2]);
    f[2] = 2.0 - exp(2.0*x[2]);
    return true;
  }
  bool get_F(const realmat& x, realmat& f, realmat& F)
  {
    if (! this->get_f(x,f) ) return false;
    if ((F.nrow() !=2)||(F.ncol()!=2)) F.resize(2,2);
    F(1,1) = f[1]-2.0; F(1,2) = f[2]-2.0;
    F(2,1) = 0.0;      F(2,2) = 2.0*(f[2]-2.0);
    return true;
  }
};

class nleqns_df : public nleqns_base {
public:
  bool get_f(const realmat& x, realmat& f) 
  {
    if (x.nrow() != 2) error("Error, nleqns, wrong dimension for x");
    if (f.nrow() != 2) f.resize(2,1);
    f[1] = 2.0 - exp(1.0*x[1] + 1.0*x[2]);
    f[2] = 2.0 - exp(2.0*x[2]);
    return true;
  }
  bool get_F(const realmat& x, realmat& f, realmat& F)
  {
    /*
    if (! this->get_f(x,f) ) return false;
    return nleqns_base::df(x,F);
    */

    if (this->get_f(x,f)) {
      return nleqns_base::df(x,F);
    }
    else {
      return false;
    }

  }
};

class nleqns_trivial : public nleqns_base {
  bool get_f(const realmat&x, realmat&f) 
  {
    f = x;
    return true;
  }
  bool get_F(const realmat& x, realmat& f, realmat& F) 
  {
    f = x;
    INTEGER d = x.size();
    F.resize(d,d,0.0);
    for (INTEGER i=1; i<=d; ++i) F(i,i) = 1.0;
    return true;
  }
};    

int main(int argc, char* argp[], char* envp[])
{
  cerr << '\n';
  cerr << "Warning: " << *argp << " writes to tstnls.out, the usage\n" 
       << "         " << *argp << " > tstnls.out is incorrect\n";

  cerr << '\n';
  cerr << "first test's error messages:" << '\n';

  ofstream fout("tstnls.out");

  fout << starbox("/Using analytic derivatives//");

  nleqns exp_eqns;
  nlsolve solver(exp_eqns);

  solver.set_output(true, &fout);
  solver.set_check_derivatives(true);
  solver.set_warning_messages(true);

  realmat x_start(2,1,0.0); realmat x_stop;
  if (solver.solve(x_start, x_stop)) {
    fout << starbox("/The Answer!//") << x_stop; 
  }
  else {
   fout << starbox("/Failure!//"); 
  }

  //#if defined GNU_GPP_COMPILER

  x_start[1] = 0.0;
  x_start[2] = -3.0;
  if (solver.solve(x_start, x_stop)) {
    fout << starbox("/The Answer!//") << x_stop; 
  }
  else {
   fout << starbox("/Failure!//"); 
  }

  //#endif

  cerr << "second test's error messages:" << '\n';

  fout << starbox("/Using numerical differentiator df//");
  
  nleqns_df exp_eqns_df;
  nlsolve solver_df(exp_eqns_df);

  solver_df.set_output(true, &fout);
  solver_df.set_check_derivatives(true);
  solver_df.set_warning_messages(true);

  realmat x_start_df(2,1,0.0); realmat x_stop_df;
  if (solver_df.solve(x_start_df, x_stop_df)) {
    fout << starbox("/The Answer!//") << x_stop_df; 
  }
  else {
   fout << starbox("/Failure!//"); 
  }

  fout << solver_df.get_norm_f() << '\n';
  fout << solver_df.get_rank_F() << '\n';
  fout << solver_df.get_iter_count() << '\n';

  cerr << "third test's error messages:" << '\n';

  fout << starbox("/Trivial problem//");

  nleqns_trivial trivial_eqns;
  nlsolve trivial_solver(trivial_eqns);

  trivial_solver.set_output(true, &fout);
  trivial_solver.set_check_derivatives(true);
  trivial_solver.set_warning_messages(true);

  realmat x_start_trivial(10,1,0.0);
  realmat x_stop_trivial;
  if (trivial_solver.solve(x_start_trivial, x_stop_trivial)) {
    fout << starbox("/The Answer!//") << x_stop_trivial; 
  }
  else {
   fout << starbox("/Failure!//"); 
  }

  fout << trivial_solver.get_norm_f() << '\n';
  fout << trivial_solver.get_rank_F() << '\n';
  fout << trivial_solver.get_iter_count() << '\n';

  x_start_trivial.resize(20,1,1.1);
  if (trivial_solver.solve(x_start_trivial, x_stop_trivial)) {
    fout << starbox("/The Answer!//") << x_stop_trivial; 
  }
  else {
   fout << starbox("/Failure!//"); 
  }

  fout << trivial_solver.get_norm_f() << '\n';
  fout << trivial_solver.get_rank_F() << '\n';
  fout << trivial_solver.get_iter_count() << '\n';

  return 0;
}
