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 / NLEqn.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  5KB  |  326 lines
  1. // NLEqn.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 <float.h>
  29. #include <math.h>
  30.  
  31. #include "NLEqn.h"
  32. #include "dMatrix.h"
  33. #include "f77-uscore.h"
  34. #include "lo-error.h"
  35.  
  36. extern "C"
  37. {
  38.   int F77_FCN (hybrd1) (int (*)(int*, double*, double*, int*),
  39.             const int*, double*, double*, const double*,
  40.             int*, double*, const int*);
  41.  
  42.   int F77_FCN (hybrj1) (int (*)(int*, double*, double*, double*, int*, int*),
  43.             const int*, double*, double*, double*, const int*,
  44.             const double*, int*, double*, const int*);
  45. }
  46.  
  47. static nonlinear_fcn user_fun;
  48. static jacobian_fcn user_jac;
  49.  
  50. // error handling
  51.  
  52. void
  53. NLEqn::error (const char* msg)
  54. {
  55.   (*current_liboctave_error_handler) ("fatal NLEqn error: %s", msg);
  56. }
  57.  
  58. // Constructors
  59.  
  60. NLEqn::NLEqn (void) : NLFunc (), n (0), x () { }
  61.  
  62. NLEqn::NLEqn (const Vector& xvec, const NLFunc f) 
  63.   : NLFunc (f), n (xvec.capacity ()), x (xvec) { }
  64.  
  65. NLEqn::NLEqn (const NLEqn& a) : NLFunc (a.fun, a.jac), n (a.n), x (a.x) { }
  66.  
  67. void
  68. NLEqn::resize (int nn)
  69. {
  70.   if (n != nn)
  71.     {
  72.       n = nn;
  73.       x.resize (n);
  74.     }
  75. }
  76.  
  77. int
  78. NLEqn::size (void) const
  79. {
  80.   return n;
  81. }
  82.  
  83. // Assignment
  84.  
  85. NLEqn&
  86. NLEqn::operator = (const NLEqn& a)
  87. {
  88.   fun = a.fun;
  89.   jac = a.jac;
  90.   x = a.n;
  91.  
  92.   return *this;
  93. }
  94.  
  95. Vector
  96. NLEqn::states (void) const
  97. {
  98.   return x;
  99. }
  100.  
  101. void
  102. NLEqn::set_states (const Vector& xvec)
  103. {
  104.   if (xvec.capacity () != n)
  105.     {
  106.       error ("dimension error");
  107.       return;
  108.     }
  109.  
  110.   x = xvec;
  111. }
  112.  
  113. // Other operations
  114.  
  115. Vector
  116. NLEqn::solve (const Vector& xvec)
  117. {
  118.   set_states (xvec);
  119.   int info;
  120.   return solve (info);
  121. }
  122.  
  123. Vector
  124. NLEqn::solve (const Vector& xvec, int& info)
  125. {
  126.   set_states (xvec);
  127.   return solve (info);
  128. }
  129.  
  130. Vector
  131. NLEqn::solve (void)
  132. {
  133.   int info;
  134.   return solve (info);
  135. }
  136.  
  137. int
  138. hybrd1_fcn (int *n, double *x, double *fvec, int *iflag)
  139. {
  140.   int nn = *n;
  141.   Vector tmp_f (nn);
  142.   Vector tmp_x (nn);
  143.  
  144.   for (int i = 0; i < nn; i++)
  145.     tmp_x.elem (i) = x[i];
  146.  
  147.   tmp_f = (*user_fun) (tmp_x);
  148.  
  149.   if (tmp_f.length () == 0)
  150.     *iflag = -1;
  151.   else
  152.     {
  153.       for (i = 0; i < nn; i++)
  154.     fvec[i] = tmp_f.elem (i);
  155.     }
  156.  
  157.   return 0;
  158. }
  159.  
  160. int
  161. hybrj1_fcn (int *n, double *x, double *fvec, double *fjac,
  162.         int *ldfjac, int *iflag)
  163. {
  164.   int nn = *n;
  165.   Vector tmp_x (nn);
  166.  
  167.   for (int i = 0; i < nn; i++)
  168.     tmp_x.elem (i) = x[i];
  169.  
  170.   int flag = *iflag;
  171.   if (flag == 1)
  172.     {
  173.       Vector tmp_f (nn);
  174.  
  175.       tmp_f = (*user_fun) (tmp_x);
  176.  
  177.       if (tmp_f.length () == 0)
  178.     *iflag = -1;
  179.       else
  180.     {
  181.       for (i = 0; i < nn; i++)
  182.         fvec[i] = tmp_f.elem (i);
  183.     }
  184.     }
  185.   else
  186.     {
  187.       Matrix tmp_fj (nn, nn);
  188.  
  189.       tmp_fj = (*user_jac) (tmp_x);
  190.  
  191.       if (tmp_fj.rows () == 0 || tmp_fj.columns () == 0)
  192.     *iflag = -1;
  193.       else
  194.     {
  195.       int ld = *ldfjac;
  196.       for (int j = 0; j < nn; j++)
  197.         for (i = 0; i < nn; i++)
  198.           fjac[j*ld+i] = tmp_fj.elem (i, j);
  199.     }
  200.     }
  201.  
  202.   return 0;
  203. }
  204.  
  205. Vector
  206. NLEqn::solve (int& info)
  207. {
  208.   int tmp_info = 0;
  209.  
  210.   if (n == 0)
  211.     {
  212.       error ("equation set not initialized");
  213.       return Vector ();
  214.     }
  215.  
  216.   double tol = tolerance ();
  217.  
  218.   double *fvec = new double [n];
  219.   double *px = new double [n];
  220.   for (int i = 0; i < n; i++)
  221.     px[i] = x.elem (i);
  222.  
  223.   user_fun = fun;
  224.   user_jac = jac;
  225.  
  226.   if (jac)
  227.     {
  228.       int lwa = (n*(n+13))/2;
  229.       double *wa = new double [lwa];
  230.       double *fjac = new double [n*n];
  231.  
  232.       F77_FCN (hybrj1) (hybrj1_fcn, &n, px, fvec, fjac, &n, &tol,
  233.             &tmp_info, wa, &lwa);
  234.  
  235.       delete [] wa;
  236.       delete [] fjac;
  237.     }
  238.   else
  239.     {
  240.       int lwa = (n*(3*n+13))/2;
  241.       double *wa = new double [lwa];
  242.  
  243.       F77_FCN (hybrd1) (hybrd1_fcn, &n, px, fvec, &tol, &tmp_info, wa, &lwa);
  244.  
  245.       delete [] wa;
  246.     }
  247.  
  248.   Vector retval;
  249.  
  250.   info = tmp_info;
  251.  
  252.   if (info >= 0)
  253.     {
  254.       retval.resize (n);
  255.  
  256.       for (i = 0; i < n; i++)
  257.     retval.elem (i) = px[i];
  258.     }
  259.  
  260.   delete [] fvec;
  261.   delete [] px;
  262.  
  263.   return retval;
  264. }
  265.  
  266. NLEqn_options::NLEqn_options (void)
  267. {
  268.   init ();
  269. }
  270.  
  271. NLEqn_options::NLEqn_options (const NLEqn_options& opt)
  272. {
  273.   copy (opt);
  274. }
  275.  
  276. NLEqn_options&
  277. NLEqn_options::operator = (const NLEqn_options& opt)
  278. {
  279.   if (this != &opt)
  280.     copy (opt);
  281.  
  282.   return *this;
  283. }
  284.  
  285. NLEqn_options::~NLEqn_options (void)
  286. {
  287. }
  288.  
  289. void
  290. NLEqn_options::init (void)
  291. {
  292.   double sqrt_eps = sqrt (DBL_EPSILON);
  293.   x_tolerance = sqrt_eps;
  294. }
  295.  
  296. void
  297. NLEqn_options::copy (const NLEqn_options& opt)
  298. {
  299.   x_tolerance = opt.x_tolerance;
  300. }
  301.  
  302. void
  303. NLEqn_options::set_default_options (void)
  304. {
  305.   init ();
  306. }
  307.  
  308. void
  309. NLEqn_options::set_tolerance (double val)
  310. {
  311.   x_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
  312. }
  313.  
  314. double
  315. NLEqn_options::tolerance (void)
  316. {
  317.   return x_tolerance;
  318. }
  319.  
  320. /*
  321. ;;; Local Variables: ***
  322. ;;; mode: C++ ***
  323. ;;; page-delimiter: "^/\\*" ***
  324. ;;; End: ***
  325. */
  326.