home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / liboctave / QPSOL.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  6KB  |  280 lines
  1. // QPSOL.cc                                                -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. #ifdef HAVE_CONFIG_H
  25. #include "config.h"
  26. #endif
  27.  
  28. #include <math.h>
  29. #include <float.h>
  30.  
  31. #ifndef QPSOL_MISSING
  32.  
  33. #include "QPSOL.h"
  34. #include "f77-uscore.h"
  35.  
  36. extern "C"
  37. {
  38.   int F77_FCN (qpsol) (int*, int*, int*, int*, int*, int*, int*, int*,
  39.                double*, double*, double*, double*, double*,
  40.                double*, double*,
  41.                int (*)(int*, int*, int*, int*, double*,
  42.                    double*, double*),
  43.                int*, int*, int*, int*, double*, int*, int*,
  44.                double*, double*, int*, int*, double*, int*);
  45.  
  46.   int F77_FCN (dgemv) (const char*, const int*, const int*,
  47.                const double*, const double*, const int*,
  48.                const double*, const int*, const double*,
  49.                double*, const int*, long);
  50. }
  51.  
  52. int
  53. qphess (int *pn, int *pnrowh, int *ncolh, int *pcol, double *hess,
  54.     double *x, double *hx)
  55. {
  56.   int n = *pn;
  57.   int nrowh = *pnrowh;
  58.   int jthcol = *pcol;
  59.  
  60.   if (jthcol > 0)
  61.     {
  62.       int hp = (jthcol - 1) * nrowh;
  63.       for (int i = 0; i < n; i++)
  64.     hx[i] = hess[hp+i];
  65.     }
  66.   else
  67.     {
  68.       char trans = 'N';
  69.       double alpha = 1.0;
  70.       double beta  = 0.0;
  71.       int i_one = 1;
  72.  
  73.       F77_FCN (dgemv) (&trans, pn, pn, &alpha, hess, pn, x, &i_one,
  74.                &beta, hx, &i_one, 1L);
  75.     }
  76.  
  77.   return 0;
  78. }
  79.  
  80. Vector
  81. QPSOL::minimize (double& objf, int& inform, Vector& lambda)
  82. {
  83.   int i;
  84.  
  85.   int n = x.capacity ();
  86.  
  87.   int itmax = (iteration_limit () < 0) ? 50 * n : iteration_limit ();
  88.   int msglvl = print_level ();
  89.   int nclin = lc.size ();
  90.   int nctotl = nclin + n;
  91.  
  92.   double bigbnd = infinite_bound ();
  93.  
  94.   double dummy;
  95.   double *pa = &dummy;
  96.   Matrix clin;
  97.   if (nclin > 0)
  98.     {
  99.       clin = lc.constraint_matrix ();
  100.       pa = clin.fortran_vec ();
  101.     }
  102.  
  103.   double *pbl = new double [nctotl];
  104.   double *pbu = new double [nctotl];
  105.  
  106.   if (bnds.size () > 0)
  107.     {
  108.       for (i = 0; i < n; i++)
  109.     {
  110.       pbl[i] = bnds.lower_bound (i);
  111.       pbu[i] = bnds.upper_bound (i);
  112.     }
  113.     }
  114.   else
  115.     {
  116.       for (i = 0; i < n; i++)
  117.     {
  118.       pbl[i] = -bigbnd;
  119.       pbu[i] = bigbnd;
  120.     }
  121.     }
  122.  
  123.   for (i = 0; i < nclin; i++)
  124.     {
  125.       pbl[i+n] = lc.lower_bound (i);
  126.       pbu[i+n] = lc.upper_bound (i);
  127.     }
  128.  
  129.   double *pc = c.fortran_vec ();
  130.  
  131.   double *featol = new double [nctotl];
  132.   double tmp = feasibility_tolerance ();
  133.   for (i = 0; i < nctotl; i++)
  134.     featol[i] = tmp;
  135.  
  136.   double *ph = H.fortran_vec ();
  137.  
  138.   int cold = 1;
  139.   int lp = 0;
  140.   int orthog = 1;
  141.  
  142.   int *istate = new int [nctotl];
  143.  
  144.   double *px = x.fortran_vec ();
  145.  
  146.   int iter = 0;
  147.   lambda.resize (nctotl);
  148.   double *pclambda = lambda.fortran_vec ();
  149.  
  150.   int leniw = 2 * n;
  151.  
  152.   int lenw;
  153.   int ncon = nclin > 1 ? nclin : 1;
  154.   if (lp == 0 || nclin >= n)
  155.     lenw = 2*n*(n + 2) + nclin + 2*ncon;
  156.   else
  157.     lenw = 2*ncon*(1 + ncon) + 4*n + nclin;
  158.  
  159.   int *iw = new int [leniw];
  160.   double *w = new double [lenw];
  161.  
  162.   F77_FCN (qpsol) (&itmax, &msglvl, &n, &nclin, &nctotl, &ncon, &n,
  163.            &n, &bigbnd, pa, pbl, pbu, pc, featol, ph, qphess,
  164.            &cold, &lp, &orthog, istate, px, &inform, &iter,
  165.            &objf, pclambda, iw, &leniw, w, &lenw);
  166.  
  167.   delete [] pbl;
  168.   delete [] pbu;
  169.   delete [] featol;
  170.   delete [] istate;
  171.   delete [] iw;
  172.   delete [] w;
  173.  
  174.   return x;
  175. }
  176.  
  177. QPSOL_options::QPSOL_options (void)
  178. {
  179.   init ();
  180. }
  181.  
  182. QPSOL_options::QPSOL_options (const QPSOL_options& opt)
  183. {
  184.   copy (opt);
  185. }
  186.  
  187. QPSOL_options&
  188. QPSOL_options::operator = (const QPSOL_options& opt)
  189. {
  190.   if (this != &opt)
  191.     copy (opt);
  192.  
  193.   return *this;
  194. }
  195.  
  196. QPSOL_options::~QPSOL_options (void)
  197. {
  198. }
  199.  
  200. void
  201. QPSOL_options::init (void)
  202. {
  203.   x_feasibility_tolerance = sqrt (DBL_EPSILON);
  204.   x_infinite_bound = 1.0e+30;
  205.   x_iteration_limit = -1;
  206.   x_print_level = 0;
  207. }
  208.  
  209. void
  210. QPSOL_options::copy (const QPSOL_options& opt)
  211. {
  212.   x_feasibility_tolerance = opt.x_feasibility_tolerance;
  213.   x_infinite_bound = opt.x_infinite_bound;
  214.   x_iteration_limit = opt.x_iteration_limit;
  215.   x_print_level = opt.x_print_level;
  216. }
  217.  
  218. void
  219. QPSOL_options::set_default_options (void)
  220. {
  221.   init ();
  222. }
  223.  
  224. void
  225. QPSOL_options::set_feasibility_tolerance (double val)
  226. {
  227.   x_feasibility_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
  228. }
  229.  
  230. void
  231. QPSOL_options::set_infinite_bound (double val)
  232. {
  233.   x_infinite_bound = (val > 0.0) ? val : 1.0e+30;
  234. }
  235.  
  236. void
  237. QPSOL_options::set_iteration_limit (int val)
  238. {
  239.   x_iteration_limit = (val > 0) ? val : -1;
  240. }
  241.  
  242. void
  243. QPSOL_options::set_print_level (int val)
  244. {
  245.   x_print_level = (val >= 0) ? val : 0;
  246. }
  247.  
  248. double
  249. QPSOL_options::feasibility_tolerance (void)
  250. {
  251.   return x_feasibility_tolerance;
  252. }
  253.  
  254. double
  255. QPSOL_options::infinite_bound (void)
  256. {
  257.   return x_infinite_bound;
  258. }
  259.  
  260. int
  261. QPSOL_options::iteration_limit (void)
  262. {
  263.   return x_iteration_limit;
  264. }
  265.  
  266. int
  267. QPSOL_options::print_level (void)
  268. {
  269.   return x_print_level;
  270. }
  271.  
  272. #endif /* QPSOL_MISSING */
  273.  
  274. /*
  275. ;;; Local Variables: ***
  276. ;;; mode: C++ ***
  277. ;;; page-delimiter: "^/\\*" ***
  278. ;;; End: ***
  279. */
  280.