#include "mtt_Hybrd_Solver.hh"

// http://www.netlib.org/minpack/hybrd.f
// used by Octave's fsolve

MTT::Hybrd_Solver *MTT::Hybrd_Solver::static_ptr;

ColumnVector
MTT::Hybrd_Solver::f_hybrd (const ColumnVector &tryUi)
{
  MTT::Hybrd_Solver::static_ptr->_yz = MTT::Hybrd_Solver::static_ptr->eval(tryUi);
  return MTT::Hybrd_Solver::static_ptr->_yz;
}

void
MTT::Hybrd_Solver::Solve (void)
{    
  static std::fstream ferr ("MTT.Hybrd_messages", std::ios::out | std::ios::trunc | std::ios::app);
  int info;
  static int input_errors;
  static int user_errors;
  static int convergences;
  static int progress_errors;
  static int limit_errors;
  static int unknown_errors;
  
  NLFunc fcn(&Hybrd_Solver::f_hybrd);
  NLEqn	 eqn(Solver::_ui,fcn);
  eqn.set_tolerance(1.0e-20);
  Solver::_ui = eqn.solve(info);

  switch (info)
    {
    case 1:
      convergences++;
      break;
    case -2:
      input_errors++;
      break;
    case -1:
      user_errors++;
      break;
    case 3:
      progress_errors++;
      break;
    case 4:
      limit_errors++;
      break;
    default:
      unknown_errors++;
      break;
    }
  std::clog.setf (std::ios::scientific);
  if (1 != info)
    {
      std::clog
	<< "\r"
	<< " time " << _t << " \t"
	<< " converge (" << convergences << ") "
	<< " limit (" << limit_errors << ")"
	<< " progress (" << progress_errors << ")"
	<< " other (" << input_errors + user_errors + unknown_errors << ") "
	<< " (max error = " << std::abs (eval(_ui).max()) << ")"
	<< std::endl;
    }
  else
    {
      std::clog
	<< "\r time " << _t << "\t max error = " << std::abs (eval(_ui).max());
    }
  ferr << info << " ";
}