Index: mttroot/mtt/bin/mtt ================================================================== --- mttroot/mtt/bin/mtt +++ mttroot/mtt/bin/mtt @@ -14,10 +14,15 @@ ############################################################### ## Version control history ############################################################### ## $Header$ ## $Log$ +## Revision 1.309.2.1 2001/05/04 04:07:24 geraint +## Numerical solution of algebraic equations. +## sys_ae.cc written for unsolved inputs. +## Solution of equations using hybrd from MINPACK (as used by Octave fsolve). +## ## Revision 1.309 2001/04/28 03:15:03 geraint ## Fixed comment (interfered with "mtt help representations"). ## ## Revision 1.308 2001/04/25 22:17:45 geraint ## Fixed icd.txt2 dependency. @@ -1088,10 +1093,13 @@ verbose=' -s' # Default integration method integration_method=implicit; +# Default algebraic equation solver +algebraic_solver=Reduce_Solver; + # Default no info info_switch='' # Default use m, not oct files m='m'; @@ -1152,10 +1160,33 @@ ;; *) echo $1 is an unknown integration method - use euler, rk4 or implicit; exit;; esac;; + -ae ) + mtt_switches="$mtt_switches $1"; + case $2 in + fsolve | hybrd) + mtt_switches="$mtt_switches $2"; + algebraic_solver=Hybrd_Solver; + shift; + ;; + hooke) + mtt_switches="$mtt_switches $2"; + algebraic_solver=HJ_Solver; + shift; + ;; + reduce) + mtt_switches="$mtt_switches $2"; + algebraic_solver=Reduce_Solver; + Solve='-A'; + shift; + ;; + *) + echo $1 is an unknown solver - use hybrd, hooke or reduce; + exit;; + esac;; -s ) mtt_switches="$mtt_switches $1"; sensitivity=sensitivity ;; -ss ) mtt_switches="$mtt_switches $1"; @@ -1288,10 +1319,11 @@ echo ' -c c-code generation' echo ' -d use directory ' echo ' -dc Maximise derivative (not integral) causality' echo ' -dc Maximise derivative (not integral) causality' echo ' -i Use implicit, euler or rk4 integration' + echo ' -ae Solve algebraic equations with Reduce, hybrd (fsolve) or "Hooke and Jeeves"' echo ' -o ode is same as dae' echo ' -oct use oct files in place of m files where appropriate' echo ' -opt optimise code generation' echo ' -p print environment variables' echo ' -partition partition hierachical system' @@ -1925,10 +1957,13 @@ mtt_m2cc.sh $1 \$* cat mtt_%.cc: ${MTT_LIB}/cc/mtt_%.cc cp ${MTT_LIB}/cc/mtt_\$*.cc mtt_\$*.cc +mtt_%.hh: ${MTT_LIB}/cc/mtt_%.hh + cp ${MTT_LIB}/cc/mtt_\$*.hh mtt_\$*.hh + ## .o files .PRECIOUS: $1_%.o $1_%.o: $1_%.cc $1_def.h $1_sympar.h $1_cr.h echo Compiling $1_\$*.cc ${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $< -DSTANDALONE @@ -2227,10 +2262,15 @@ @echo "Creating $1_ode2odes_implicit.o" ar -cr \$@ \$^ $1_ode2odes_implicit.oct: $1_cseo.oct $1_csex.oct $1_smxa.oct $1_smxax.oct mtt_implicit.oct @echo > /dev/null +mtt_Solver.cc: mtt_Solver.hh +$1_ode2odes_${algebraic_solver}.cc: mtt_Solver.cc mtt_${algebraic_solver}.hh mtt_${algebraic_solver}.cc +$1_ode2odes_${algebraic_solver}.o: mtt_Solver.o mtt_${algebraic_solver}.o + @echo "Creating $1_ode2odes_${algebraic_solver}.o" + ar -cr \$@ \$^ #SUMMARY numpar numerical parameter declaration (m) $1_numpar.m: $1_numpar.txt $1_sympars.txt mtt_txt2m $1 numpar @@ -2529,11 +2569,11 @@ $1_csex.m $1_cseo.m $1_logic.m ifeq ($using_oct,yes) touch $1_ode2odes.m # Create a dummy which wont' be used mtt $mtt_switches -q -u $1 ode2odes oct else - make_ode2odes $1 m $integration_method + make_ode2odes $1 m $integration_method $algebraic_solver endif endif ifneq ($integration_method,implicit) $1_ode2odes.m : $1_def.r $1_sympars.txt\ $1_simpar.m $1_numpar.m $1_state.m $1_input.m \ @@ -2541,27 +2581,27 @@ ifeq ($using_oct,yes) echo "*** Warning: Shouldn't be here! Creating dummy $1_ode2odes.m" touch $1_ode2odes.m # Create a dummy which wont' be used mtt $mtt_switches -q -u $1 ode2odes oct else - make_ode2odes $1 m $integration_method + make_ode2odes $1 m $integration_method $algebraic_solver endif endif #SUMMARY ode2odes Simulation function (m) #SUMMARY ode2odes Simulation function (cc) #SUMMARY ode2odes Simulation function (oct) #SUMMARY ode2odes Simulation function (exe) $1_ode2odes.exe: $1_def.h $1_sympar.h\ - $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o + $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o echo Creating $1_ode2odes.exe ${MTT_CXX} ${MTT_CXXFLAGS} -o $1_ode2odes.exe\ - $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o\ + $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o\ ${MTT_LDFLAGS} ${MTT_CXXLIBS} -$1_ode2odes.o: $1_ode2odes.cc $1_ode2odes_common.o $1_ode2odes_${integration_method}.o +$1_ode2odes.o: $1_ode2odes.cc $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o echo Creating $1_ode2odes.o ${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $1_ode2odes.cc -DSTANDALONE $1_ode2odes.oct: $1_ode2odes.cc $1_ode2odes_common.oct $1_ode2odes_${integration_method}.oct touch $1_ode2odes.m @@ -2568,13 +2608,13 @@ echo Creating $1_ode2odes.oct $MKOCTFILE $1_ode2odes.cc $1_ode2odes.cc: $1_def.r $1_sympars.txt\ $1_ode2odes_common.m $1_ode2odes_common.cc\ - $1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc + $1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc mtt_Solver.cc mtt_${algebraic_solver}.cc mtt_${algebraic_solver}.hh touch $1_ode2odes.m - make_ode2odes $1 cc $integration_method + make_ode2odes $1 cc $integration_method $algebraic_solver #Conversion of m to p to c #SUMMARY ode ordinary differential equations (c) #SUMMARY ode ordinary differential equations (p) #SUMMARY state state declaration (c) Index: mttroot/mtt/bin/trans/make_ode2odes ================================================================== --- mttroot/mtt/bin/trans/make_ode2odes +++ mttroot/mtt/bin/trans/make_ode2odes @@ -7,10 +7,15 @@ ############################################################### ## Version control history ############################################################### ## $Id$ ## $Log$ +## Revision 1.57.2.1 2001/05/04 04:07:24 geraint +## Numerical solution of algebraic equations. +## sys_ae.cc written for unsolved inputs. +## Solution of equations using hybrd from MINPACK (as used by Octave fsolve). +## ## Revision 1.57 2001/04/01 03:38:54 geraint ## Reset row to zero after write to file, ready for subsequent runs. ## Eliminates SIGSEGV in Octave when _ode2odes called multiple times. ## ## Revision 1.56 2001/03/30 15:13:58 gawthrop @@ -223,10 +228,16 @@ if [ -n "$3" ]; then method=$3 else method=implicit fi + +if [ -n "$4" ]; then + algebraic_solver=$4 +else + algebraic_solver="Reduce_Solver" +fi echo Creating $filename with $method integration method # Find system constants Nx=`mtt_getsize $sys x` # States @@ -348,11 +359,10 @@ cat < $filename #include #include #include -#include #include #include #ifndef STANDALONE #include @@ -362,10 +372,12 @@ #include "${sys}_sympar.h" #ifdef STANDALONE #include #include + +#include "mtt_${algebraic_solver}.hh" extern ColumnVector F${sys}_ae ( ColumnVector &x, ColumnVector &u, const double &t, @@ -444,108 +456,10 @@ void set_signal_handlers (void); #endif // STANDALONE -class Solver -{ -public: - Solver (void) { - solve_errors.input = 0; - solve_errors.user = 0; - solve_errors.converge = 0; - solve_errors.progress = 0; - solve_errors.limit = 0; - solve_errors.unknown = 0; - } - ColumnVector solve (const ColumnVector &x, - const ColumnVector &u, - const double &t, - const ColumnVector &par) - { - initialise(x,u,t,par); - NLFunc fcn(&Solver::evaluate); - NLEqn eqn(Ui,fcn); - Ui = eqn.solve(info); - switch (info) - { - case -2: - solve_errors.input++; - write_stats(); - break; - case -1: - solve_errors.user++; - write_stats(); - break; - case 1: - solve_errors.converge++; - break; - case 3: - solve_errors.progress++; - write_stats(); - break; - case 4: - solve_errors.limit++; - write_stats(); - break; - default: - solve_errors.unknown++; - write_stats(); - break;; - } - U.insert (Ui,MTTNU); - return U; - } - void write_stats (void) - { - cerr << "input error (" << solve_errors.input << ") " - << ", user (" << solve_errors.user << ") " - << ", converge (" << solve_errors.converge << ") " - << ", progress (" << solve_errors.progress << ") " - << ", limit (" << solve_errors.limit << ") " - << ", unknown (" << solve_errors.unknown << ") " - << endl; - } -private: - static ColumnVector evaluate (const ColumnVector &tryUi) - { - U.insert(tryUi,MTTNU); - return F${sys}_ae(X,U,T,P); - } - void initialise (const ColumnVector &x, - const ColumnVector &u, - const double &t, - const ColumnVector &par) - { - X = x; - T = t; - P = par; - U.insert(u,0); - Ui = ColumnVector(MTTNYZ,0.0); - } -private: - struct { - long int input; - long int user; - long int converge; - long int progress; - long int limit; - long int unknown; - } solve_errors; - static double T; - static ColumnVector X; - static ColumnVector P; - static ColumnVector U; - static ColumnVector Ui; - int info; -}; -double Solver::T = 0.0; -ColumnVector Solver::U(MTTNU+MTTNYZ,0.0); -ColumnVector Solver::X(MTTNX,0.0); -ColumnVector Solver::P(MTTNPAR,0.0); -ColumnVector Solver::Ui(MTTNYZ,0.0); - inline ColumnVector mtt_input (ColumnVector &x, ColumnVector &y, const double &t, ColumnVector &par) @@ -552,11 +466,12 @@ { #ifdef STANDALONE static ColumnVector u (MTTNU); static ColumnVector ui (MTTNYZ); static ColumnVector U (MTTNU+MTTNYZ); - static Solver ae; + + static ${algebraic_solver} ae(F${sys}_ae,MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ); u = F${sys}_input (x, y, t, par); if (MTTNYZ == 0) { return u; Index: mttroot/mtt/cc/include/useful-functions.hh ================================================================== --- mttroot/mtt/cc/include/useful-functions.hh +++ mttroot/mtt/cc/include/useful-functions.hh @@ -4,36 +4,33 @@ #define HAVE_USEFUL_FUNCTIONS_HH #ifndef __cplusplus #define inline /* strip */ -typedef double doubleref_t; -#else -typedef double &doubleref_t; #endif // ! __cplusplus static inline double -max (const doubleref_t x1, const doubleref_t x2) +max (const double x1, const double x2) { return static_cast((x1 >= x2) ? x1 : (x1 < x2) ? x2 : 0); } static inline double -min (const doubleref_t x1, const doubleref_t x2) +min (const double x1, const double x2) { return static_cast((x1 <= x2) ? x1 : (x1 > x2) ? x2 : 0); } static inline double -nonsingular (const doubleref_t x) +nonsingular (const double x) { return static_cast((x == 0) ? 1.0e-30 : x); } static inline double -sign (const doubleref_t x) +sign (const double x) { return static_cast((x > 0) ? +1 : (x < 0) ? -1 : 0); } ADDED mttroot/mtt/lib/cc/mtt_HJ_Solver.cc Index: mttroot/mtt/lib/cc/mtt_HJ_Solver.cc ================================================================== --- /dev/null +++ mttroot/mtt/lib/cc/mtt_HJ_Solver.cc @@ -0,0 +1,279 @@ + +#include "mtt_HJ_Solver.hh" + +// http://www.netlib.org/opt/hooke.c +// Hooke and Jeeves solution + +HJ_Solver *HJ_Solver::static_ptr; + +double +HJ_Solver::f (double tryUi[], int nyz) +{ + for (int i = 0; i < nyz; i++) + HJ_Solver::static_ptr->_ui(i) = tryUi[i]; + + HJ_Solver::static_ptr->_yz = HJ_Solver::static_ptr->eval(HJ_Solver::static_ptr->_ui); + + double ans = 0.0; + for (int i = 0; i < nyz; i++) + ans += HJ_Solver::static_ptr->_yz(i)*HJ_Solver::static_ptr->_yz(i); + return ans; +} + +void +HJ_Solver::Solve (void) +{ + double startpt[_nyz]; + double endpt[_nyz]; + for (int i = 0; i < _nyz; i++) + startpt[i] = _ui(i); + hooke(_nyz,startpt,endpt); + for (int i = 0; i < _nyz; i++) + _ui(i) = endpt[i]; +} + + // adapted from http://www.netlib.org/opt/hooke.c: + +/* Nonlinear Optimization using the algorithm of Hooke and Jeeves */ +/* 12 February 1994 author: Mark G. Johnson */ + + +/* Find a point X where the nonlinear function f(X) has a local */ +/* minimum. X is an n-vector and f(X) is a scalar. In mathe- */ +/* matical notation f: R^n -> R^1. The objective function f() */ +/* is not required to be continuous. Nor does f() need to be */ +/* differentiable. The program does not use or require */ +/* derivatives of f(). */ + +/* The software user supplies three things: a subroutine that */ +/* computes f(X), an initial "starting guess" of the minimum point */ +/* X, and values for the algorithm convergence parameters. Then */ +/* the program searches for a local minimum, beginning from the */ +/* starting guess, using the Direct Search algorithm of Hooke and */ +/* Jeeves. */ + +/* This C program is adapted from the Algol pseudocode found in */ +/* "Algorithm 178: Direct Search" by Arthur F. Kaupe Jr., Commun- */ +/* ications of the ACM, Vol 6. p.313 (June 1963). It includes the */ +/* improvements suggested by Bell and Pike (CACM v.9, p. 684, Sept */ +/* 1966) and those of Tomlin and Smith, "Remark on Algorithm 178" */ +/* (CACM v.12). The original paper, which I don't recommend as */ +/* highly as the one by A. Kaupe, is: R. Hooke and T. A. Jeeves, */ +/* "Direct Search Solution of Numerical and Statistical Problems", */ +/* Journal of the ACM, Vol. 8, April 1961, pp. 212-229. */ + +/* Calling sequence: */ +/* int hooke(nvars, startpt, endpt, rho, epsilon, itermax) */ +/* */ +/* nvars {an integer} This is the number of dimensions */ +/* in the domain of f(). It is the number of */ +/* coordinates of the starting point (and the */ +/* minimum point.) */ +/* startpt {an array of doubles} This is the user- */ +/* supplied guess at the minimum. */ +/* endpt {an array of doubles} This is the location of */ +/* the local minimum, calculated by the program */ +/* rho {a double} This is a user-supplied convergence */ +/* parameter (more detail below), which should be */ +/* set to a value between 0.0 and 1.0. Larger */ +/* values of rho give greater probability of */ +/* convergence on highly nonlinear functions, at a */ +/* cost of more function evaluations. Smaller */ +/* values of rho reduces the number of evaluations */ +/* (and the program running time), but increases */ +/* the risk of nonconvergence. See below. */ +/* epsilon {a double} This is the criterion for halting */ +/* the search for a minimum. When the algorithm */ +/* begins to make less and less progress on each */ +/* iteration, it checks the halting criterion: if */ +/* the stepsize is below epsilon, terminate the */ +/* iteration and return the current best estimate */ +/* of the minimum. Larger values of epsilon (such */ +/* as 1.0e-4) give quicker running time, but a */ +/* less accurate estimate of the minimum. Smaller */ +/* values of epsilon (such as 1.0e-7) give longer */ +/* running time, but a more accurate estimate of */ +/* the minimum. */ +/* itermax {an integer} A second, rarely used, halting */ +/* criterion. If the algorithm uses >= itermax */ +/* iterations, halt. */ + + +/* The user-supplied objective function f(x,n) should return a C */ +/* "double". Its arguments are x -- an array of doubles, and */ +/* n -- an integer. x is the point at which f(x) should be */ +/* evaluated, and n is the number of coordinates of x. That is, */ +/* n is the number of coefficients being fitted. */ + +/* rho, the algorithm convergence control */ +/* The algorithm works by taking "steps" from one estimate of */ +/* a minimum, to another (hopefully better) estimate. Taking */ +/* big steps gets to the minimum more quickly, at the risk of */ +/* "stepping right over" an excellent point. The stepsize is */ +/* controlled by a user supplied parameter called rho. At each */ +/* iteration, the stepsize is multiplied by rho (0 < rho < 1), */ +/* so the stepsize is successively reduced. */ +/* Small values of rho correspond to big stepsize changes, */ +/* which make the algorithm run more quickly. However, there */ +/* is a chance (especially with highly nonlinear functions) */ +/* that these big changes will accidentally overlook a */ +/* promising search vector, leading to nonconvergence. */ +/* Large values of rho correspond to small stepsize changes, */ +/* which force the algorithm to carefully examine nearby points */ +/* instead of optimistically forging ahead. This improves the */ +/* probability of convergence. */ +/* The stepsize is reduced until it is equal to (or smaller */ +/* than) epsilon. So the number of iterations performed by */ +/* Hooke-Jeeves is determined by rho and epsilon: */ +/* rho**(number_of_iterations) = epsilon */ +/* In general it is a good idea to set rho to an aggressively */ +/* small value like 0.5 (hoping for fast convergence). Then, */ + + +/* if the user suspects that the reported minimum is incorrect */ +/* (or perhaps not accurate enough), the program can be run */ +/* again with a larger value of rho such as 0.85, using the */ +/* result of the first minimization as the starting guess to */ +/* begin the second minimization. */ + +/* Normal use: (1) Code your function f() in the C language */ +/* (2) Install your starting guess {or read it in} */ +/* (3) Run the program */ +/* (4) {for the skeptical}: Use the computed minimum */ +/* as the starting point for another run */ + +/* Data Fitting: */ +/* Code your function f() to be the sum of the squares of the */ +/* errors (differences) between the computed values and the */ +/* measured values. Then minimize f() using Hooke-Jeeves. */ +/* EXAMPLE: you have 20 datapoints (ti, yi) and you want to */ +/* find A,B,C such that (A*t*t) + (B*exp(t)) + (C*tan(t)) */ +/* fits the data as closely as possible. Then f() is just */ +/* f(x) = SUM (measured_y[i] - ((A*t[i]*t[i]) + (B*exp(t[i])) */ +/* + (C*tan(t[i]))))^2 */ +/* where x[] is a 3-vector consisting of {A, B, C}. */ + +/* */ +/* The author of this software is M.G. Johnson. */ +/* Permission to use, copy, modify, and distribute this software */ +/* for any purpose without fee is hereby granted, provided that */ +/* this entire notice is included in all copies of any software */ +/* which is or includes a copy or modification of this software */ +/* and in all copies of the supporting documentation for such */ +/* software. THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT */ +/* ANY EXPRESS OR IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE */ +/* AUTHOR NOR AT&T MAKE ANY REPRESENTATION OR WARRANTY OF ANY */ +/* KIND CONCERNING THE MERCHANTABILITY OF THIS SOFTWARE OR ITS */ +/* FITNESS FOR ANY PARTICULAR PURPOSE. */ +/* */ + +double +HJ_Solver::best_nearby(double *delta, + double *point, + double prevbest, + int nvars) +{ + double z[VARS]; + double minf, ftmp; + int i; + minf = prevbest; + for (i = 0; i < nvars; i++) + z[i] = point[i]; + for (i = 0; i < nvars; i++) { + z[i] = point[i] + delta[i]; + ftmp = f(z, nvars); + if (ftmp < minf) + minf = ftmp; + else { + delta[i] = 0.0 - delta[i]; + z[i] = point[i] + delta[i]; + ftmp = f(z, nvars); + if (ftmp < minf) + minf = ftmp; + else + z[i] = point[i]; + } + } + for (i = 0; i < nvars; i++) + point[i] = z[i]; + return (minf); +} + + +int +HJ_Solver::hooke (int nvars, // MTTNYZ + double startpt[], // user's initial guess + double endpt[], // result + double rho, // geometric shrink factor + double epsilon, // end value stepsize + int itermax) // max # iterations +{ + const int VARS = _nyz; + double delta[VARS]; + double newf, fbefore, steplength, tmp; + double xbefore[VARS], newx[VARS]; + int i, j, keep; + int iters, iadj; + for (i = 0; i < nvars; i++) { + newx[i] = xbefore[i] = startpt[i]; + delta[i] = fabs(startpt[i] * rho); + if (delta[i] == 0.0) + delta[i] = rho; + } + iadj = 0; + steplength = rho; + iters = 0; + fbefore = f(newx, nvars); + newf = fbefore; + while ((iters < itermax) && (steplength > epsilon)) { + iters++; + iadj++; + /* find best new point, one coord at a time */ + for (i = 0; i < nvars; i++) { + newx[i] = xbefore[i]; + } + newf = best_nearby(delta, newx, fbefore, nvars); + /* if we made some improvements, pursue that direction */ + keep = 1; + while ((newf < fbefore) && (keep == 1)) { + iadj = 0; + for (i = 0; i < nvars; i++) { + /* firstly, arrange the sign of delta[] */ + if (newx[i] <= xbefore[i]) + delta[i] = 0.0 - fabs(delta[i]); + else + delta[i] = fabs(delta[i]); /* now, move further in this direction */ + tmp = xbefore[i]; + xbefore[i] = newx[i]; + newx[i] = newx[i] + newx[i] - tmp; + } + fbefore = newf; + newf = best_nearby(delta, newx, fbefore, nvars); + /* if the further (optimistic) move was bad.... */ + if (newf >= fbefore) + break; + /* make sure that the differences between the new */ + /* and the old points are due to actual */ + /* displacements; beware of roundoff errors that */ + /* might cause newf < fbefore */ + keep = 0; + for (i = 0; i < nvars; i++) { + keep = 1; + if (fabs(newx[i] - xbefore[i]) > + (0.5 * fabs(delta[i]))) + break; + else + keep = 0; + } + } + if ((steplength >= epsilon) && (newf >= fbefore)) { + steplength = steplength * rho; + for (i = 0; i < nvars; i++) { + delta[i] *= rho; + } + } +} +for (i = 0; i < nvars; i++) + endpt[i] = xbefore[i]; +return (iters); +} ADDED mttroot/mtt/lib/cc/mtt_HJ_Solver.hh Index: mttroot/mtt/lib/cc/mtt_HJ_Solver.hh ================================================================== --- /dev/null +++ mttroot/mtt/lib/cc/mtt_HJ_Solver.hh @@ -0,0 +1,51 @@ + +#include "mtt_Solver.hh" + +class HJ_Solver : public Solver { + + // http://www.netlib.org/opt/hooke.c + // Hooke and Jeeves solution + +public: + + HJ_Solver (sys_ae ae, + const int npar, + const int nu, + const int nx, + const int ny, + const int nyz) + : Solver (ae,npar,nu,nx,ny,nyz) + { static_ptr = this; VARS = nyz; }; + + static double + f (double tryUi[], int nyz); + +protected: + + void + Solve (void); + + double + best_nearby (double delta[], + double point[], + double prevbest, + int nvars); + + int + hooke (int nvars, // MTTNYZ + double startpt[], // user's initial guess + double endpt[], // result + double rho = 0.05, // geometric shrink factor + double epsilon = 1e-3, // end value stepsize + int itermax = 5000); // max # iterations + +private: + + int VARS; + +public: + + static HJ_Solver *static_ptr; + +}; + ADDED mttroot/mtt/lib/cc/mtt_Hybrd_Solver.cc Index: mttroot/mtt/lib/cc/mtt_Hybrd_Solver.cc ================================================================== --- /dev/null +++ mttroot/mtt/lib/cc/mtt_Hybrd_Solver.cc @@ -0,0 +1,67 @@ + +#include "mtt_Hybrd_Solver.hh" + +// http://www.netlib.org/minpack/hybrd.f +// used by Octave's fsolve + +Hybrd_Solver *Hybrd_Solver::static_ptr; + +ColumnVector +Hybrd_Solver::f_hybrd (const ColumnVector &tryUi) +{ + Hybrd_Solver::static_ptr->_yz = Hybrd_Solver::static_ptr->eval(tryUi); + return Hybrd_Solver::static_ptr->_yz; +} + +void +Hybrd_Solver::Solve (void) +{ + 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(0.01); + 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: + // if (abs(eval(_ui).max()) > 1.0e-6) + limit_errors++; + // else + // convergences++; + break; + default: + unknown_errors++; + break; + } + if (1 != info) + { + cerr << "input (" << input_errors << ") " + << " user (" << user_errors << ") " + << " converge (" << convergences << ") " + << " progress (" << progress_errors << ") " + << " limit (" << limit_errors << ")" + << " unknown (" << unknown_errors << ")" + << " (max error = " << abs(eval(_ui).max()) << ")" << endl; + } +} + ADDED mttroot/mtt/lib/cc/mtt_Hybrd_Solver.hh Index: mttroot/mtt/lib/cc/mtt_Hybrd_Solver.hh ================================================================== --- /dev/null +++ mttroot/mtt/lib/cc/mtt_Hybrd_Solver.hh @@ -0,0 +1,36 @@ + +#include "mtt_Solver.hh" +#include + +class Hybrd_Solver : public Solver { + + // http://www.netlib.org/minpack/hybrd.f + // used by Octave's fsolve + +public: + + Hybrd_Solver (sys_ae ae, + const int npar, + const int nu, + const int nx, + const int ny, + const int nyz) + : Solver (ae,npar,nu,nx,ny,nyz) + { + static_ptr = this; + } + + static ColumnVector + f_hybrd (const ColumnVector &tryUi); + +protected: + + void + Solve (void); + +public: + + static Hybrd_Solver *static_ptr; + +}; + ADDED mttroot/mtt/lib/cc/mtt_Reduce_Solver.cc Index: mttroot/mtt/lib/cc/mtt_Reduce_Solver.cc ================================================================== --- /dev/null +++ mttroot/mtt/lib/cc/mtt_Reduce_Solver.cc @@ -0,0 +1,11 @@ + +#include "mtt_Reduce_Solver.hh" + + +void +Reduce_Solver::Solve (void) +{ + cerr << "Error:" + << " Symbolic solution of equations failed during model build" << endl + << " Try using one of the other algebraic solution methods" << endl; +} ADDED mttroot/mtt/lib/cc/mtt_Reduce_Solver.hh Index: mttroot/mtt/lib/cc/mtt_Reduce_Solver.hh ================================================================== --- /dev/null +++ mttroot/mtt/lib/cc/mtt_Reduce_Solver.hh @@ -0,0 +1,25 @@ + +#include "mtt_Solver.hh" + +class Reduce_Solver : public Solver { + + // Dummy class + // This will not be used unless the Reduce solver has failed earlier + // in the model build process + +public: + + Reduce_Solver (sys_ae ae, + const int npar, + const int nu, + const int nx, + const int ny, + const int nyz) + : Solver (ae,npar,nu,nx,ny,nyz) + { ; }; + + void + Solve (void); + +}; + ADDED mttroot/mtt/lib/cc/mtt_Solver.cc Index: mttroot/mtt/lib/cc/mtt_Solver.cc ================================================================== --- /dev/null +++ mttroot/mtt/lib/cc/mtt_Solver.cc @@ -0,0 +1,42 @@ + +#include "mtt_Solver.hh" + +Solver::Solver (sys_ae ae, + const int npar, + const int nu, + const int nx, + const int ny, + const int nyz) +{ + _ae = ae; + _np = npar; + _nu = nu; + _nx = nx; + _ny = ny; + _nyz = nyz; + + _uui = ColumnVector (_nu+_nyz); +}; + +ColumnVector +Solver::solve (const ColumnVector &x, + const ColumnVector &u, + const double &t, + const ColumnVector &par) +{ + _x = x; + _uui.insert(u,0); + _t = t; + _par = par; + _ui = ColumnVector(_nyz,0.2); + Solve (); + _uui.insert(_ui,_nu); + return _uui; +} + +ColumnVector +Solver::eval (const ColumnVector &ui) +{ + _uui.insert(ui,_nu); + return _ae (_x, _uui, _t, _par); +} ADDED mttroot/mtt/lib/cc/mtt_Solver.hh Index: mttroot/mtt/lib/cc/mtt_Solver.hh ================================================================== --- /dev/null +++ mttroot/mtt/lib/cc/mtt_Solver.hh @@ -0,0 +1,60 @@ + +#ifndef MTT_SOLVER +#define MTT_SOLVER + +#include +#include +#include + +#include + +class Solver { + + typedef ColumnVector (*sys_ae) // pointer to F${sys}_ae function + (ColumnVector &,ColumnVector &,const double &t,ColumnVector &); + +public: + + Solver (sys_ae ae, + const int npar, + const int nu, + const int nx, + const int ny, + const int nyz); + + ColumnVector + solve (const ColumnVector &x, + const ColumnVector &u, + const double &t, + const ColumnVector &par); + + ColumnVector + eval (const ColumnVector &ui); + +protected: + + virtual void + Solve (void) = 0; + +protected: + + ColumnVector _x; + ColumnVector _uui; + double _t; + ColumnVector _par; + + ColumnVector _ui; + ColumnVector _yz; + + int _nu; + int _np; + int _nx; + int _ny; + int _nyz; + + sys_ae _ae; + +}; + +#endif // MTT_SOLVER +