Mercurial > hg > octave-terminal
changeset 1945:8c4bce5e773e
[project @ 1996-02-14 01:03:03 by jwe]
author | jwe |
---|---|
date | Wed, 14 Feb 1996 01:03:40 +0000 |
parents | 8cb4d3008c76 |
children | ec07d85b4152 |
files | liboctave/DASSL.cc liboctave/DASSL.h liboctave/LSODE.cc liboctave/LSODE.h liboctave/NPSOL.cc liboctave/QPSOL.cc |
diffstat | 6 files changed, 184 insertions(+), 190 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/DASSL.cc +++ b/liboctave/DASSL.cc @@ -63,12 +63,10 @@ liw = 0; lrw = 0; - info = new int [15]; - iwork = (int *) 0; - rwork = (double *) 0; + info.resize (15); for (int i = 0; i < 15; i++) - info [i] = 0; + info.elem (i) = 0; } DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f) @@ -82,12 +80,10 @@ liw = 20 + n; lrw = 40 + 9*n + n*n; - info = new int [15]; - iwork = new int [liw]; - rwork = new double [lrw]; + info.resize (15); for (int i = 0; i < 15; i++) - info [i] = 0; + info.elem (i) = 0; } DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv, @@ -105,19 +101,10 @@ liw = 20 + n; lrw = 40 + 9*n + n*n; - info = new int [15]; - iwork = new int [liw]; - rwork = new double [lrw]; + info.resize (15); for (int i = 0; i < 15; i++) - info [i] = 0; -} - -DASSL::~DASSL (void) -{ - delete [] info; - delete [] rwork; - delete [] iwork; + info.elem (i) = 0; } void @@ -199,12 +186,26 @@ ColumnVector DASSL::do_integrate (double tout) { + ColumnVector retval; + + if (restart) + { + restart = 0; + info.elem (0) = 0; + } + + if (iwork.length () != liw) + iwork.resize (liw); + + if (rwork.length () != lrw) + rwork.resize (lrw); + integration_error = 0; if (DAEFunc::jacobian_function ()) - iwork [4] = 1; + iwork.elem (4) = 1; else - iwork [4] = 0; + iwork.elem (4) = 0; double *px = x.fortran_vec (); double *pxdot = xdot.fortran_vec (); @@ -215,93 +216,91 @@ if (stop_time_set) { - info [3] = 1; - rwork [0] = stop_time; + info.elem (3) = 1; + rwork.elem (0) = stop_time; } else - info [3] = 0; + info.elem (3) = 0; double abs_tol = absolute_tolerance (); double rel_tol = relative_tolerance (); if (initial_step_size () >= 0.0) { - rwork[2] = initial_step_size (); - info[7] = 1; + rwork.elem (2) = initial_step_size (); + info.elem (7) = 1; } else - info[7] = 0; + info.elem (7) = 0; if (maximum_step_size () >= 0.0) { - rwork[2] = maximum_step_size (); - info[6] = 1; + rwork.elem (2) = maximum_step_size (); + info.elem (6) = 1; } else - info[6] = 0; + info.elem (6) = 0; double *dummy = 0; int *idummy = 0; - if (restart) - { - restart = 0; - info[0] = 0; - } + int *pinfo = info.fortran_vec (); + int *piwork = iwork.fortran_vec (); + double *prwork = rwork.fortran_vec (); // again: - F77_FCN (ddassl, DDASSL) (ddassl_f, n, t, px, pxdot, tout, info, - rel_tol, abs_tol, idid, rwork, lrw, iwork, - liw, dummy, idummy, ddassl_j); + F77_XFCN (ddassl, DDASSL, (ddassl_f, n, t, px, pxdot, tout, pinfo, + rel_tol, abs_tol, idid, prwork, lrw, + piwork, liw, dummy, idummy, ddassl_j)); - switch (idid) + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dassl"); + else { - case 1: // A step was successfully taken in the - // intermediate-output mode. The code has not yet reached - // TOUT. - break; + switch (idid) + { + case 1: // A step was successfully taken in intermediate-output + // mode. The code has not yet reached TOUT. + case 2: // The integration to TSTOP was successfully completed + // (T=TSTOP) by stepping exactly to TSTOP. + case 3: // The integration to TOUT was successfully completed + // (T=TOUT) by stepping past TOUT. Y(*) is obtained by + // interpolation. YPRIME(*) is obtained by interpolation. - case 2: // The integration to TSTOP was successfully completed - // (T=TSTOP) by stepping exactly to TSTOP. - break; - - case 3: // The integration to TOUT was successfully completed - // (T=TOUT) by stepping past TOUT. Y(*) is obtained by - // interpolation. YPRIME(*) is obtained by interpolation. - break; + retval = x; + t = tout; + break; - case -1: // A large amount of work has been expended. (About 500 steps). - case -2: // The error tolerances are too stringent. - case -3: // The local error test cannot be satisfied because you - // specified a zero component in ATOL and the - // corresponding computed solution component is zero. - // Thus, a pure relative error test is impossible for - // this component. - case -6: // DDASSL had repeated error test failures on the last - // attempted step. - case -7: // The corrector could not converge. - case -8: // The matrix of partial derivatives is singular. - case -9: // The corrector could not converge. There were repeated - // error test failures in this step. - case -10: // The corrector could not converge because IRES was - // equal to minus one. - case -11: // IRES equal to -2 was encountered and control is being - // returned to the calling program. - case -12: // DDASSL failed to compute the initial YPRIME. - case -33: // The code has encountered trouble from which it cannot - // recover. A message is printed explaining the trouble - // and control is returned to the calling program. For - // example, this occurs when invalid input is detected. - - default: - integration_error = 1; - break; + case -1: // A large amount of work has been expended. (~500 steps). + case -2: // The error tolerances are too stringent. + case -3: // The local error test cannot be satisfied because you + // specified a zero component in ATOL and the + // corresponding computed solution component is zero. + // Thus, a pure relative error test is impossible for + // this component. + case -6: // DDASSL had repeated error test failures on the last + // attempted step. + case -7: // The corrector could not converge. + case -8: // The matrix of partial derivatives is singular. + case -9: // The corrector could not converge. There were repeated + // error test failures in this step. + case -10: // The corrector could not converge because IRES was + // equal to minus one. + case -11: // IRES equal to -2 was encountered and control is being + // returned to the calling program. + case -12: // DDASSL failed to compute the initial YPRIME. + case -33: // The code has encountered trouble from which it cannot + // recover. A message is printed explaining the trouble + // and control is returned to the calling program. For + // example, this occurs when invalid input is detected. + default: + integration_error = 1; + break; + } } - t = tout; - - return x; + return retval; } Matrix
--- a/liboctave/DASSL.h +++ b/liboctave/DASSL.h @@ -119,7 +119,7 @@ DASSL (const ColumnVector& x, const ColumnVector& xdot, double time, DAEFunc& f); - ~DASSL (void); + ~DASSL (void) { } void force_restart (void); @@ -146,9 +146,9 @@ int liw; int lrw; int idid; - int *info; - int *iwork; - double *rwork; + Array<int> info; + Array<int> iwork; + Array<double> rwork; friend int ddassl_j (double *time, double *state, double *deriv, double *pd, double *cj, double *rpar, int *ipar);
--- a/liboctave/LSODE.cc +++ b/liboctave/LSODE.cc @@ -72,14 +72,6 @@ liw = 20 + n; lrw = 22 + n * (9 + n); - - iwork = new int [liw]; - rwork = new double [lrw]; - for (int i = 4; i < 9; i++) - { - iwork[i] = 0; - rwork[i] = 0.0; - } } LSODE::LSODE (const ColumnVector& state, double time, const ODEFunc& f) @@ -100,20 +92,6 @@ liw = 20 + n; lrw = 22 + n * (9 + n); - - iwork = new int [liw]; - rwork = new double [lrw]; - for (int i = 4; i < 9; i++) - { - iwork[i] = 0; - rwork[i] = 0.0; - } -} - -LSODE::~LSODE (void) -{ - delete [] rwork; - delete [] iwork; } void @@ -180,6 +158,30 @@ ColumnVector LSODE::do_integrate (double tout) { + ColumnVector retval; + + if (restart) + { + restart = 0; + istate = 1; + } + + if (iwork.length () != liw) + { + iwork.resize (liw); + + for (int i = 4; i < 9; i++) + iwork.elem (i) = 0; + } + + if (rwork.length () != lrw) + { + rwork.resize (lrw); + + for (int i = 4; i < 9; i++) + rwork.elem (i) = 0; + } + if (jac) method_flag = 21; else @@ -199,13 +201,12 @@ // Try 5000 steps before giving up. - iwork[5] = 5000; - int working_too_hard = 0; + iwork.elem (5) = 5000; if (stop_time_set) { itask = 4; - rwork [0] = stop_time; + rwork.elem (0) = stop_time; } else { @@ -215,64 +216,67 @@ double abs_tol = absolute_tolerance (); double rel_tol = relative_tolerance (); - rwork[4] = (initial_step_size () >= 0.0) ? initial_step_size () : 0.0; - rwork[5] = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0; - rwork[6] = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0; + rwork.elem (4) = (initial_step_size () >= 0.0) ? initial_step_size () : 0.0; + rwork.elem (5) = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0; + rwork.elem (6) = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0; - if (restart) - { - restart = 0; - istate = 1; - } + int *piwork = iwork.fortran_vec (); + double *prwork = rwork.fortran_vec (); + + working_too_hard = 0; again: - F77_FCN (lsode, LSODE) (lsode_f, n, xp, t, tout, itol, rel_tol, - abs_tol, itask, istate, iopt, rwork, lrw, - iwork, liw, lsode_j, method_flag); + F77_XFCN (lsode, LSODE, (lsode_f, n, xp, t, tout, itol, rel_tol, + abs_tol, itask, istate, iopt, prwork, lrw, + piwork, liw, lsode_j, method_flag)); - switch (istate) + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in lsode"); + else { - case -13: // Return requested in user-supplied function. - case -6: // error weight became zero during problem. (solution - // component i vanished, and atol or atol(i) = 0.) - case -5: // repeated convergence failures (perhaps bad jacobian - // supplied or wrong choice of mf or tolerances). - case -4: // repeated error test failures (check all inputs). - case -3: // illegal input detected (see printed message). - case -2: // excess accuracy requested (tolerances too small). - integration_error = 1; - return ColumnVector (); - break; + switch (istate) + { + case -13: // Return requested in user-supplied function. + case -6: // error weight became zero during problem. (solution + // component i vanished, and atol or atol(i) = 0.) + case -5: // repeated convergence failures (perhaps bad jacobian + // supplied or wrong choice of mf or tolerances). + case -4: // repeated error test failures (check all inputs). + case -3: // illegal input detected (see printed message). + case -2: // excess accuracy requested (tolerances too small). + integration_error = 1; + break; - case -1: // excess work done on this call (perhaps wrong mf). - working_too_hard++; - if (working_too_hard > 20) - { + case -1: // excess work done on this call (perhaps wrong mf). + working_too_hard++; + if (working_too_hard > 20) + { + (*current_liboctave_error_handler) + ("giving up after more than %d steps attempted in lsode", + iwork.elem (5) * 20); + integration_error = 1; + } + else + { + istate = 2; + goto again; + } + break; + + case 2: // lsode was successful + retval = x; + t = tout; + break; + + default: // Error? (*current_liboctave_error_handler) - ("giving up after more than %d steps attempted in lsode", - iwork[5] * 20); - integration_error = 1; - return ColumnVector (); + ("unrecognized value of istate returned from lsode"); + break; } - else - { - istate = 2; - goto again; - } - break; - - case 2: // lsode was successful - break; - - default: - // Error? - break; } - t = tout; - - return x; + return retval; } #if 0
--- a/liboctave/LSODE.h +++ b/liboctave/LSODE.h @@ -128,7 +128,7 @@ LSODE (const ColumnVector& state, double time, const ODEFunc& f); - ~LSODE (void); + ~LSODE (void) { } void force_restart (void); @@ -157,14 +157,15 @@ int integration_error; int restart; int method_flag; - int *iwork; - double *rwork; + Array<int> iwork; + Array<double> rwork; int istate; int itol; int itask; int iopt; int liw; int lrw; + int working_too_hard; friend int lsode_f (int *neq, double *t, double *y, double *ydot);
--- a/liboctave/NPSOL.cc +++ b/liboctave/NPSOL.cc @@ -183,45 +183,40 @@ // Constraint stuff. + double bigbnd = infinite_bound (); + Matrix clin = lc.constraint_matrix (); double *pclin = clin.fortran_vec (); - Array<double> aclow (n+nclin+ncnln); - double *clow = aclow.fortran_vec (); - - Array<double> acup (n+nclin+ncnln); - double *cup = acup.fortran_vec (); + ColumnVector aclow (n+nclin+ncnln); + ColumnVector acup (n+nclin+ncnln); if (bnds.size () > 0) { - for (int i = 0; i < n; i++) - { - clow[i] = bnds.lower_bound (i); - cup[i] = bnds.upper_bound (i); - } + aclow.insert (bnds.lower_bounds (), 0); + acup.insert (bnds.upper_bounds (), 0); } else { - double huge = 1.0e30; - for (int i = 0; i < n; i++) - { - clow[i] = -huge; - cup[i] = huge; - } + aclow.fill (-bigbnd, 0, n-1); + acup.fill (bigbnd, 0, n-1); } - for (int i = 0; i < nclin; i++) + if (nclin > 0) { - clow[i+n] = lc.lower_bound (i); - cup[i+n] = lc.upper_bound (i); + aclow.insert (lc.lower_bounds (), n); + acup.insert (lc.upper_bounds (), n); } - for (int i = 0; i < ncnln; i++) + if (ncnln > 0) { - clow[i+n+nclin] = nlc.lower_bound (i); - cup[i+n+nclin] = nlc.upper_bound (i); + aclow.insert (nlc.lower_bounds (), n+nclin); + acup.insert (nlc.upper_bounds (), n+nclin); } + double *clow = aclow.fortran_vec (); + double *cup = acup.fortran_vec (); + Array<double> ac (ncnln); double *c = ac.fortran_vec ();
--- a/liboctave/QPSOL.cc +++ b/liboctave/QPSOL.cc @@ -87,13 +87,8 @@ double bigbnd = infinite_bound (); - double *pa = 0; - Matrix clin; - if (nclin > 0) - { - clin = lc.constraint_matrix (); - pa = clin.fortran_vec (); - } + Matrix clin = lc.constraint_matrix (); + double *pa = clin.fortran_vec (); ColumnVector bl (n+nclin); ColumnVector bu (n+nclin); @@ -105,8 +100,8 @@ } else { - bl.fill (-bigbnd); - bu.fill (bigbnd); + bl.fill (-bigbnd, 0, n-1); + bu.fill (bigbnd, 0, n-1); } if (nclin > 0)