Mercurial > hg > octave-terminal
changeset 1546:a272c4056bab
[project @ 1995-10-07 03:38:40 by jwe]
author | jwe |
---|---|
date | Sat, 07 Oct 1995 03:39:26 +0000 |
parents | 8bdfa6fe5d69 |
children | 75e822649946 |
files | liboctave/DASSL.cc liboctave/LinConst.cc liboctave/Quad.cc liboctave/Range.cc |
diffstat | 4 files changed, 103 insertions(+), 78 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/DASSL.cc +++ b/liboctave/DASSL.cc @@ -102,7 +102,7 @@ info [i] = 0; } -DAE::DAE (const Vector& state, double time, DAEFunc& f) +DAE::DAE (const ColumnVector& state, double time, DAEFunc& f) { n = state.capacity (); t = time; @@ -129,7 +129,8 @@ info [i] = 0; } -DAE::DAE (const Vector& state, const Vector& deriv, double time, DAEFunc& f) +DAE::DAE (const ColumnVector& state, const ColumnVector& deriv, + double time, DAEFunc& f) { if (deriv.capacity () != state.capacity ()) { @@ -169,14 +170,14 @@ delete iwork; } -Vector +ColumnVector DAE::deriv (void) { return xdot; } void -DAE::initialize (const Vector& state, double time) +DAE::initialize (const ColumnVector& state, double time) { integration_error = 0; restart = 1; @@ -187,7 +188,8 @@ } void -DAE::initialize (const Vector& state, const Vector& deriv, double time) +DAE::initialize (const ColumnVector& state, + const ColumnVector& deriv, double time) { integration_error = 0; restart = 1; @@ -200,9 +202,9 @@ ddassl_f (const double& time, double *state, double *deriv, double *delta, int& ires, double *, int *) { - Vector tmp_deriv (nn); - Vector tmp_state (nn); - Vector tmp_delta (nn); + ColumnVector tmp_deriv (nn); + ColumnVector tmp_state (nn); + ColumnVector tmp_delta (nn); for (int i = 0; i < nn; i++) { @@ -227,8 +229,8 @@ ddassl_j (const double& time, double *, double *, double *pd, const double& cj, double *, int *) { - Vector tmp_state (nn); - Vector tmp_deriv (nn); + ColumnVector tmp_state (nn); + ColumnVector tmp_deriv (nn); // XXX FIXME XXX @@ -252,7 +254,7 @@ return 0; } -Vector +ColumnVector DAE::integrate (double tout) { integration_error = 0; @@ -361,7 +363,7 @@ } Matrix -DAE::integrate (const Vector& tout, Matrix& xdot_out) +DAE::integrate (const ColumnVector& tout, Matrix& xdot_out) { Matrix retval; int n_out = tout.capacity (); @@ -396,7 +398,8 @@ } Matrix -DAE::integrate (const Vector& tout, Matrix& xdot_out, const Vector& tcrit) +DAE::integrate (const ColumnVector& tout, Matrix& xdot_out, + const ColumnVector& tcrit) { Matrix retval; int n_out = tout.capacity ();
--- a/liboctave/LinConst.cc +++ b/liboctave/LinConst.cc @@ -42,8 +42,8 @@ (*current_liboctave_error_handler) ("fatal LinConst error: %s", msg); } -LinConst::LinConst (const Matrix& a_eq, const Vector& b_eq, - const Matrix& a_ineq, const Vector& b_ineq) +LinConst::LinConst (const Matrix& a_eq, const ColumnVector& b_eq, + const Matrix& a_ineq, const ColumnVector& b_ineq) { // Need some checks here. @@ -123,10 +123,10 @@ return retval; } -Vector +ColumnVector LinConst::eq_constraint_vector (void) const { - Vector retval (nb); + ColumnVector retval (nb); int count = 0; for (int i = 0; i < nb; i++) { @@ -140,10 +140,10 @@ return retval; } -Vector +ColumnVector LinConst::ineq_constraint_vector (void) const { - Vector retval (2*nb); + ColumnVector retval (2*nb); int count = 0; for (int i = 0; i < nb; i++) {
--- a/liboctave/Quad.cc +++ b/liboctave/Quad.cc @@ -29,9 +29,6 @@ #include <config.h> #endif -#include <cfloat> -#include <cmath> - #include "Quad.h" #include "f77-uscore.h" #include "sun-utils.h"
--- a/liboctave/Range.cc +++ b/liboctave/Range.cc @@ -29,7 +29,9 @@ #include <config.h> #endif +#include <cfloat> #include <climits> +#include <cmath> #include <iostream.h> @@ -136,6 +138,81 @@ return is; } +// C See Knuth, Art Of Computer Programming, Vol. 1, Problem 1.2.4-5. +// C +// C===Tolerant FLOOR function. +// C +// C X - is given as a Double Precision argument to be operated on. +// C It is assumed that X is represented with M mantissa bits. +// C CT - is given as a Comparison Tolerance such that +// C 0.LT.CT.LE.3-SQRT(5)/2. If the relative difference between +// C X and A whole number is less than CT, then TFLOOR is +// C returned as this whole number. By treating the +// C floating-point numbers as a finite ordered set note that +// C the heuristic EPS=2.**(-(M-1)) and CT=3*EPS causes +// C arguments of TFLOOR/TCEIL to be treated as whole numbers +// C if they are exactly whole numbers or are immediately +// C adjacent to whole number representations. Since EPS, the +// C "distance" between floating-point numbers on the unit +// C interval, and M, the number of bits in X'S mantissa, exist +// C on every floating-point computer, TFLOOR/TCEIL are +// C consistently definable on every floating-point computer. +// C +// C For more information see the following references: +// C (1) P. E. Hagerty, "More On Fuzzy Floor And Ceiling," APL QUOTE +// C QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5. +// C (2) L. M. Breed, "Definitions For Fuzzy Floor And Ceiling", APL +// C QUOTE QUAD 8(3):16-23, March 1978. This paper cites FL1 through +// C FL5, the history of five years of evolutionary development of +// C FL5 - the seven lines of code below - by open collaboration +// C and corroboration of the mathematical-computing community. +// C +// C Penn State University Center for Academic Computing +// C H. D. Knoble - August, 1978. + +static inline double +tfloor (double x, double ct) +{ +// C---------FLOOR(X) is the largest integer algebraically less than +// C or equal to X; that is, the unfuzzy FLOOR function. + +// DINT (X) = X - DMOD (X, 1.0); +// FLOOR (X) = DINT (X) - DMOD (2.0 + DSIGN (1.0, X), 3.0); + +// C---------Hagerty's FL5 function follows... + + double q = 1.0; + + if (x < 0.0) + q = 1.0 - ct; + + double rmax = q / (2.0 - ct); + + double t1 = 1.0 + floor (x); + t1 = (ct / q) * (t1 < 0.0 ? -t1 : t1); + t1 = rmax < t1 ? rmax : t1; + t1 = ct > t1 ? ct : t1; + t1 = floor (x + t1); + + if (x <= 0.0 || t1 - x < rmax) + return t1; + else + return t1 - 1.0; +} + +static inline double +tceil (double x, double ct) +{ + return -tfloor (-x, ct); +} + +static inline double +round (double x, double ct) +{ + return tfloor (x+0.5, ct); +} + + // Find an approximate number of intervals, then do the best we can to // find the number of intervals that we would get if we had done // something like @@ -152,65 +229,13 @@ int Range::nelem_internal (void) const { - // We can't have more than INT_MAX elements in the range. - - double d_n_intervals = (rng_limit - rng_base) / rng_inc; - int max_intervals = INT_MAX - 1; - double d_max_val = (double) max_intervals; - - if (d_n_intervals > d_max_val) - return -1; + double ct = 3.0 * DBL_EPSILON; - int n_intervals = (d_n_intervals > 0) - ? ((int) (d_n_intervals + 0.5)) - : ((int) (d_n_intervals - 0.5)); - - if (rng_limit > rng_base && rng_inc > 0) - { - // Our approximation may have been too big. - - while (rng_base + n_intervals * rng_inc > rng_limit && n_intervals > 0) - n_intervals--; - - // Now that we are close, get the actual number. Try to avoid - // problems with extended precision registers. + double tmp = round ((rng_limit - rng_base + rng_inc) / rng_inc, ct); - for (;;) - { - volatile double tmp_inc = (n_intervals + 1) * rng_inc; - volatile double tmp_val = rng_base + tmp_inc; - if (tmp_val <= rng_limit && n_intervals < max_intervals) - n_intervals++; - else - break; - } - } - else if (rng_limit < rng_base && rng_inc < 0) - { - // Our approximation may have been too big. - - while (rng_base + n_intervals * rng_inc < rng_limit && n_intervals > 0) - n_intervals--; + int n_intervals = (int) (tmp > 0.0 ? tmp : 0); - // Now that we are close, get the actual number. Try to avoid - // problems with extended precision registers. - - for (;;) - { - volatile double tmp_inc = (n_intervals + 1) * rng_inc; - volatile double tmp_val = rng_base + tmp_inc; - if (tmp_val >= rng_limit && n_intervals < max_intervals) - n_intervals++; - else - break; - } - } - else if (rng_limit == rng_base) - n_intervals = 0; - else - n_intervals = -1; - - return (n_intervals >= max_intervals) ? -1 : n_intervals + 1; + return (n_intervals >= INT_MAX - 1) ? -1 : n_intervals; } /*