// QPSOL.cc                                                -*- C++ -*-
/*

Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton

This file is part of Octave.

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

*/

#if defined (__GNUG__)
#pragma implementation
#endif

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <math.h>
#include <float.h>

#ifndef QPSOL_MISSING

#include "QPSOL.h"
#include "f77-uscore.h"

extern "C"
{
  int F77_FCN (qpsol, QPSOL) (int&, int&, int&, int&, int&, int&,
			      int&, int&, double&, double*, double*,
			      double*, double*, double*, double*,
			      int (*)(const int&, const int&,
				      const int&, const int&, double*,
				      double*, double*),
			      int&, int&, int&, int*, double*, int&,
			      int&, double&, double*, int*, int&,
			      double*, int&);

  int F77_FCN (dgemv, DGEMV) (const char*, const int&, const int&,
			      const double&, const double*,
			      const int&, const double*, const int&,
			      const double&, double*, const int&,
			      long);
}

int
qphess (const int& n, const int& nrowh, const int& ncolh,
	const int& jthcol, double *hess, double *x, double *hx)
{
  if (jthcol > 0)
    {
      int hp = (jthcol - 1) * nrowh;
      for (int i = 0; i < n; i++)
	hx[i] = hess[hp+i];
    }
  else
    {
      F77_FCN (dgemv, DGEMV) ("N", n, n, 1.0, hess, n, x, 1, 0.0, hx,
			      1, 1L);
    }

  return 0;
}

Vector
QPSOL::minimize (double& objf, int& inform, Vector& lambda)
{
  int i;

  int n = x.capacity ();
 
  int itmax = (iteration_limit () < 0) ? 50 * n : iteration_limit ();
  int msglvl = print_level ();
  int nclin = lc.size ();
  int nctotl = nclin + n;

  double bigbnd = infinite_bound ();

  double dummy;
  double *pa = &dummy;
  Matrix clin;
  if (nclin > 0)
    {
      clin = lc.constraint_matrix ();
      pa = clin.fortran_vec ();
    }

  double *pbl = new double [nctotl];
  double *pbu = new double [nctotl];

  if (bnds.size () > 0)
    {
      for (i = 0; i < n; i++)
	{
	  pbl[i] = bnds.lower_bound (i);
	  pbu[i] = bnds.upper_bound (i);
	}
    }
  else
    {
      for (i = 0; i < n; i++)
	{
	  pbl[i] = -bigbnd;
	  pbu[i] = bigbnd;
	}
    }

  for (i = 0; i < nclin; i++)
    {
      pbl[i+n] = lc.lower_bound (i);
      pbu[i+n] = lc.upper_bound (i);
    }

  double *pc = c.fortran_vec ();

  double *featol = new double [nctotl];
  double tmp = feasibility_tolerance ();
  for (i = 0; i < nctotl; i++)
    featol[i] = tmp;

  double *ph = H.fortran_vec ();

  int cold = 1;
  int lp = 0;
  int orthog = 1;

  int *istate = new int [nctotl];

  double *px = x.fortran_vec ();

  int iter = 0;
  lambda.resize (nctotl);
  double *pclambda = lambda.fortran_vec ();

  int leniw = 2 * n;

  int lenw;
  int ncon = nclin > 1 ? nclin : 1;
  if (lp == 0 || nclin >= n)
    lenw = 2*n*(n + 2) + nclin + 2*ncon;
  else
    lenw = 2*ncon*(1 + ncon) + 4*n + nclin;

  int *iw = new int [leniw];
  double *w = new double [lenw];

  F77_FCN (qpsol, QPSOL) (itmax, msglvl, n, nclin, nctotl, ncon, n,
			  n, bigbnd, pa, pbl, pbu, pc, featol, ph,
			  qphess, cold, lp, orthog, istate, px,
			  inform, iter, objf, pclambda, iw, leniw, w,
			  lenw);

  delete [] pbl;
  delete [] pbu;
  delete [] featol;
  delete [] istate;
  delete [] iw;
  delete [] w;

  return x;
}

QPSOL_options::QPSOL_options (void)
{
  init ();
}

QPSOL_options::QPSOL_options (const QPSOL_options& opt)
{
  copy (opt);
}

QPSOL_options&
QPSOL_options::operator = (const QPSOL_options& opt)
{
  if (this != &opt)
    copy (opt);

  return *this;
}

QPSOL_options::~QPSOL_options (void)
{
}

void
QPSOL_options::init (void)
{
  x_feasibility_tolerance = sqrt (DBL_EPSILON);
  x_infinite_bound = 1.0e+30;
  x_iteration_limit = -1;
  x_print_level = 0;
}

void
QPSOL_options::copy (const QPSOL_options& opt)
{
  x_feasibility_tolerance = opt.x_feasibility_tolerance;
  x_infinite_bound = opt.x_infinite_bound;
  x_iteration_limit = opt.x_iteration_limit;
  x_print_level = opt.x_print_level;
}

void
QPSOL_options::set_default_options (void)
{
  init ();
}

void
QPSOL_options::set_feasibility_tolerance (double val)
{
  x_feasibility_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
}

void
QPSOL_options::set_infinite_bound (double val)
{
  x_infinite_bound = (val > 0.0) ? val : 1.0e+30;
}

void
QPSOL_options::set_iteration_limit (int val)
{
  x_iteration_limit = (val > 0) ? val : -1;
}

void
QPSOL_options::set_print_level (int val)
{
  x_print_level = (val >= 0) ? val : 0;
}

double
QPSOL_options::feasibility_tolerance (void)
{
  return x_feasibility_tolerance;
}

double
QPSOL_options::infinite_bound (void)
{
  return x_infinite_bound;
}

int
QPSOL_options::iteration_limit (void)
{
  return x_iteration_limit;
}

int
QPSOL_options::print_level (void)
{
  return x_print_level;
}

#endif /* QPSOL_MISSING */

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; page-delimiter: "^/\\*" ***
;;; End: ***
*/