Mercurial > hg > octave-thorsten

changeset 3912:f56cd411adb4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-04-28 03:12:27 by jwe]
author jwe
date Sun, 28 Apr 2002 03:12:28 +0000
parents 8389e78e67d4
children f11cfbfc84b5
files kpathsea/c-auto.in libcruft/ChangeLog libcruft/slatec-err/ixsav.f libcruft/slatec-err/xerrwd.f liboctave/ChangeLog liboctave/DASPK.cc liboctave/DASPK.h liboctave/Makefile.in src/ChangeLog src/DLD-FUNCTIONS/daspk.cc src/Makefile.in
diffstat 11 files changed, 1359 insertions(+), 43 deletions(-) [+]
line wrap: on
line diff
--- a/kpathsea/c-auto.in
+++ b/kpathsea/c-auto.in
@@ -1,4 +1,4 @@
-/* c-auto.in.  Generated automatically from configure.in by autoheader.  */
+/* c-auto.in.  Generated from configure.in by autoheader.  */
 /* acconfig.h -- used by autoheader when generating c-auto.in.
 
    If you're thinking of editing acconfig.h to fix a configuration
@@ -22,7 +22,6 @@
 #undef HAVE_STRSTR
 
 
-
 /* Define if your compiler understands prototypes.  */
 #undef HAVE_PROTOTYPES
 
@@ -90,88 +89,106 @@
    of features).  */
 #undef NOTOOL
 
-/* Define if the `closedir' function returns void instead of `int'. */
+/* Define to 1 if the `closedir' function returns void instead of `int'. */
 #undef CLOSEDIR_VOID
 
-/* Define if you have the <assert.h> header file. */
+/* Define to 1 if you have the <assert.h> header file. */
 #undef HAVE_ASSERT_H
 
-/* Define if you have the `basename' function. */
+/* Define to 1 if you have the `basename' function. */
 #undef HAVE_BASENAME
 
-/* Define if you have the `bcopy' function. */
+/* Define to 1 if you have the `bcopy' function. */
 #undef HAVE_BCOPY
 
-/* Define if you have the <dirent.h> header file, and it defines `DIR'. */
+/* Define to 1 if you have the <dirent.h> header file, and it defines `DIR'.
+   */
 #undef HAVE_DIRENT_H
 
-/* Define if you have the <float.h> header file. */
+/* Define to 1 if you have the <float.h> header file. */
 #undef HAVE_FLOAT_H
 
-/* Define if you have the `getcwd' function. */
+/* Define to 1 if you have the `getcwd' function. */
 #undef HAVE_GETCWD
 
-/* Define if you have the `getwd' function. */
+/* Define to 1 if you have the `getwd' function. */
 #undef HAVE_GETWD
 
-/* Define if you have the <inttypes.h> header file. */
+/* Define to 1 if you have the <inttypes.h> header file. */
 #undef HAVE_INTTYPES_H
 
-/* Define if you have the <limits.h> header file. */
+/* Define to 1 if you have the <limits.h> header file. */
 #undef HAVE_LIMITS_H
 
-/* Define if you have the <memory.h> header file. */
+/* Define to 1 if you have the <memory.h> header file. */
 #undef HAVE_MEMORY_H
 
-/* Define if you have the <ndir.h> header file, and it defines `DIR'. */
+/* Define to 1 if you have the <ndir.h> header file, and it defines `DIR'. */
 #undef HAVE_NDIR_H
 
-/* Define if you have the `putenv' function. */
+/* Define to 1 if you have the `putenv' function. */
 #undef HAVE_PUTENV
 
-/* Define if you have the <pwd.h> header file. */
+/* Define to 1 if you have the <pwd.h> header file. */
 #undef HAVE_PWD_H
 
-/* Define if you have the <stdint.h> header file. */
+/* Define to 1 if you have the <stdint.h> header file. */
 #undef HAVE_STDINT_H
 
-/* Define if you have the <stdlib.h> header file. */
+/* Define to 1 if you have the <stdlib.h> header file. */
 #undef HAVE_STDLIB_H
 
-/* Define if you have the `strcasecmp' function. */
+/* Define to 1 if you have the `strcasecmp' function. */
 #undef HAVE_STRCASECMP
 
-/* Define if you have the <strings.h> header file. */
+/* Define to 1 if you have the <strings.h> header file. */
 #undef HAVE_STRINGS_H
 
-/* Define if you have the <string.h> header file. */
+/* Define to 1 if you have the <string.h> header file. */
 #undef HAVE_STRING_H
 
-/* Define if you have the `strstr' function. */
+/* Define to 1 if you have the `strstr' function. */
 #undef HAVE_STRSTR
 
-/* Define if you have the `strtol' function. */
+/* Define to 1 if you have the `strtol' function. */
 #undef HAVE_STRTOL
 
-/* Define if you have the <sys/dir.h> header file, and it defines `DIR'. */
+/* Define to 1 if you have the <sys/dir.h> header file, and it defines `DIR'.
+   */
 #undef HAVE_SYS_DIR_H
 
-/* Define if you have the <sys/ndir.h> header file, and it defines `DIR'. */
+/* Define to 1 if you have the <sys/ndir.h> header file, and it defines `DIR'.
+   */
 #undef HAVE_SYS_NDIR_H
 
-/* Define if you have the <sys/param.h> header file. */
+/* Define to 1 if you have the <sys/param.h> header file. */
 #undef HAVE_SYS_PARAM_H
 
-/* Define if you have the <sys/stat.h> header file. */
+/* Define to 1 if you have the <sys/stat.h> header file. */
 #undef HAVE_SYS_STAT_H
 
-/* Define if you have the <sys/types.h> header file. */
+/* Define to 1 if you have the <sys/types.h> header file. */
 #undef HAVE_SYS_TYPES_H
 
-/* Define if you have the <unistd.h> header file. */
+/* Define to 1 if you have the <unistd.h> header file. */
 #undef HAVE_UNISTD_H
 
-/* Define if you have the ANSI C header files. */
+/* Define to the address where bug reports for this package should be sent. */
+#undef PACKAGE_BUGREPORT
+
+/* Define to the full name of this package. */
+#undef PACKAGE_NAME
+
+/* Define to the full name and version of this package. */
+#undef PACKAGE_STRING
+
+/* Define to the one symbol short name of this package. */
+#undef PACKAGE_TARNAME
+
+/* Define to the version of this package. */
+#undef PACKAGE_VERSION
+
+/* Define to 1 if you have the ANSI C header files. */
 #undef STDC_HEADERS
 
 /* Define to empty if `const' does not conform to ANSI C. */
--- a/libcruft/ChangeLog
+++ b/libcruft/ChangeLog
@@ -1,5 +1,7 @@
 2002-04-27  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* slatec-err/ixsav.f, slatec-err/xerrwd.f: New files.
+
 	* daspk: New directory.
 	* Makefile.in (CRUFT_DIRS): Add it to the list
 
new file mode 100644
--- /dev/null
+++ b/libcruft/slatec-err/ixsav.f
@@ -0,0 +1,70 @@
+*DECK IXSAV
+      INTEGER FUNCTION IXSAV (IPAR, IVALUE, ISET)
+C***BEGIN PROLOGUE  IXSAV
+C***SUBSIDIARY
+C***PURPOSE  Save and recall error message control parameters.
+C***LIBRARY   MATHLIB
+C***CATEGORY  R3C
+C***TYPE      ALL (IXSAV-A)
+C***AUTHOR  Hindmarsh, Alan C., (LLNL)
+C***DESCRIPTION
+C
+C  IXSAV saves and recalls one of two error message parameters:
+C    LUNIT, the logical unit number to which messages are printed, and
+C    MESFLG, the message print flag.
+C  This is a modification of the SLATEC library routine J4SAVE.
+C
+C  Saved local variables..
+C   LUNIT  = Logical unit number for messages.
+C   LUNDEF = Default logical unit number, data-loaded to 6 below
+C            (may be machine-dependent).
+C   MESFLG = Print control flag..
+C            1 means print all messages (the default).
+C            0 means no printing.
+C
+C  On input..
+C    IPAR   = Parameter indicator (1 for LUNIT, 2 for MESFLG).
+C    IVALUE = The value to be set for the parameter, if ISET = .TRUE.
+C    ISET   = Logical flag to indicate whether to read or write.
+C             If ISET = .TRUE., the parameter will be given
+C             the value IVALUE.  If ISET = .FALSE., the parameter
+C             will be unchanged, and IVALUE is a dummy argument.
+C
+C  On return..
+C    IXSAV = The (old) value of the parameter.
+C
+C***SEE ALSO  XERMSG, XERRWD, XERRWV
+C***ROUTINES CALLED  NONE
+C***REVISION HISTORY  (YYMMDD)
+C   921118  DATE WRITTEN
+C   930329  Modified prologue to SLATEC format. (FNF)
+C   941025  Minor modification re default unit number. (ACH)
+C***END PROLOGUE  IXSAV
+C
+C**End
+      LOGICAL ISET
+      INTEGER IPAR, IVALUE
+C-----------------------------------------------------------------------
+      INTEGER LUNIT, LUNDEF, MESFLG
+C-----------------------------------------------------------------------
+C The following Fortran-77 declaration is to cause the values of the
+C listed (local) variables to be saved between calls to this routine.
+C-----------------------------------------------------------------------
+      SAVE LUNIT, LUNDEF, MESFLG
+      DATA LUNIT/-1/, LUNDEF/6/, MESFLG/1/
+C
+C***FIRST EXECUTABLE STATEMENT  IXSAV
+      IF (IPAR .EQ. 1) THEN
+        IF (LUNIT .EQ. -1) LUNIT = LUNDEF
+        IXSAV = LUNIT
+        IF (ISET) LUNIT = IVALUE
+        ENDIF
+C
+      IF (IPAR .EQ. 2) THEN
+        IXSAV = MESFLG
+        IF (ISET) MESFLG = IVALUE
+        ENDIF
+C
+      RETURN
+C----------------------- End of Function IXSAV -------------------------
+      END
new file mode 100644
--- /dev/null
+++ b/libcruft/slatec-err/xerrwd.f
@@ -0,0 +1,97 @@
+
+*DECK XERRWD
+      SUBROUTINE XERRWD (MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
+C***BEGIN PROLOGUE  XERRWD
+C***SUBSIDIARY
+C***PURPOSE  Write error message with values.
+C***LIBRARY   MATHLIB
+C***CATEGORY  R3C
+C***TYPE      DOUBLE PRECISION (XERRWV-S, XERRWD-D)
+C***AUTHOR  Hindmarsh, Alan C., (LLNL)
+C***DESCRIPTION
+C
+C  Subroutines XERRWD, XSETF, XSETUN, and the function routine IXSAV,
+C  as given here, constitute a simplified version of the SLATEC error
+C  handling package.
+C
+C  All arguments are input arguments.
+C
+C  MSG    = The message (character array).
+C  NMES   = The length of MSG (number of characters).
+C  NERR   = The error number (not used).
+C  LEVEL  = The error level..
+C           0 or 1 means recoverable (control returns to caller).
+C           2 means fatal (run is aborted--see note below).
+C  NI     = Number of integers (0, 1, or 2) to be printed with message.
+C  I1,I2  = Integers to be printed, depending on NI.
+C  NR     = Number of reals (0, 1, or 2) to be printed with message.
+C  R1,R2  = Reals to be printed, depending on NR.
+C
+C  Note..  this routine is machine-dependent and specialized for use
+C  in limited context, in the following ways..
+C  1. The argument MSG is assumed to be of type CHARACTER, and
+C     the message is printed with a format of (1X,A).
+C  2. The message is assumed to take only one line.
+C     Multi-line messages are generated by repeated calls.
+C  3. If LEVEL = 2, control passes to the statement   STOP
+C     to abort the run.  This statement may be machine-dependent.
+C  4. R1 and R2 are assumed to be in double precision and are printed
+C     in D21.13 format.
+C
+C***ROUTINES CALLED  IXSAV
+C***REVISION HISTORY  (YYMMDD)
+C   920831  DATE WRITTEN
+C   921118  Replaced MFLGSV/LUNSAV by IXSAV. (ACH)
+C   930329  Modified prologue to SLATEC format. (FNF)
+C   930407  Changed MSG from CHARACTER*1 array to variable. (FNF)
+C   930922  Minor cosmetic change. (FNF)
+C***END PROLOGUE  XERRWD
+C
+C*Internal Notes:
+C
+C For a different default logical unit number, IXSAV (or a subsidiary
+C routine that it calls) will need to be modified.
+C For a different run-abort command, change the statement following
+C statement 100 at the end.
+C-----------------------------------------------------------------------
+C Subroutines called by XERRWD.. None
+C Function routine called by XERRWD.. IXSAV
+C-----------------------------------------------------------------------
+C**End
+C
+C  Declare arguments.
+C
+      DOUBLE PRECISION R1, R2
+      INTEGER NMES, NERR, LEVEL, NI, I1, I2, NR
+      CHARACTER*(*) MSG
+C
+C  Declare local variables.
+C
+      INTEGER LUNIT, IXSAV, MESFLG
+C
+C  Get logical unit number and message print flag.
+C
+C***FIRST EXECUTABLE STATEMENT  XERRWD
+      LUNIT = IXSAV (1, 0, .FALSE.)
+      MESFLG = IXSAV (2, 0, .FALSE.)
+      IF (MESFLG .EQ. 0) GO TO 100
+C
+C  Write the message.
+C
+      WRITE (LUNIT,10)  MSG
+ 10   FORMAT(1X,A)
+      IF (NI .EQ. 1) WRITE (LUNIT, 20) I1
+ 20   FORMAT(6X,'In above message,  I1 =',I10)
+      IF (NI .EQ. 2) WRITE (LUNIT, 30) I1,I2
+ 30   FORMAT(6X,'In above message,  I1 =',I10,3X,'I2 =',I10)
+      IF (NR .EQ. 1) WRITE (LUNIT, 40) R1
+ 40   FORMAT(6X,'In above message,  R1 =',D21.13)
+      IF (NR .EQ. 2) WRITE (LUNIT, 50) R1,R2
+ 50   FORMAT(6X,'In above,  R1 =',D21.13,3X,'R2 =',D21.13)
+C
+C  Abort the run if LEVEL = 2.
+C
+ 100  IF (LEVEL .NE. 2) RETURN
+      STOP
+C----------------------- End of Subroutine XERRWD ----------------------
+      END
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,8 @@
+2002-04-27  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* DASPK.h, DASPK.cc: New files.
+	* Makefile.in: Add them to the appropriate lists.
+
 2002-04-23  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* Array2-idx.h (Array2<T>::index (idx_vector&, idx_vector&) const):
new file mode 100644
--- /dev/null
+++ b/liboctave/DASPK.cc
@@ -0,0 +1,515 @@
+/*
+
+Copyright (C) 1996, 1997, 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if defined (__GNUG__)
+#pragma implementation
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cfloat>
+#include <cmath>
+
+#include "DASPK.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+
+typedef int (*daspk_fcn_ptr) (const double&, const double*,
+			      const double*, const double&,
+			      double*, int&, double*, int*);
+
+typedef int (*daspk_jac_ptr) (const double&, const double*,
+			      const double*, double*,
+			      const double&, double*, int*);
+
+typedef int (*daspk_psol_ptr) (const int&, const double&,
+			       const double*, const double*,
+			       const double*, const double&,
+			       const double*, double*, int*,
+			       double*, const double&, int&,
+			       double*, int*);
+
+extern "C"
+int F77_FUNC (ddaspk, DDASPK) (daspk_fcn_ptr, const int&, double&,
+			      double*, double*, double&, const int*,
+			      const double&, const double&, int&,
+			      double*, const int&, int*, const int&,
+			      const double*, const int*,
+			      daspk_jac_ptr, daspk_psol_ptr);
+
+static DAEFunc::DAERHSFunc user_fun;
+static DAEFunc::DAEJacFunc user_jac;
+static int nn;
+
+DASPK::DASPK (void) : DAE ()
+{
+  stop_time_set = 0;
+  stop_time = 0.0;
+
+  sanity_checked = 0;
+
+  info.resize (15);
+
+  for (int i = 0; i < 15; i++)
+    info.elem (i) = 0;
+}
+
+DASPK::DASPK (const ColumnVector& state, double time, DAEFunc& f)
+  : DAE (state, time, f)
+{
+  n = size ();
+
+  stop_time_set = 0;
+  stop_time = 0.0;
+
+  sanity_checked = 0;
+
+  info.resize (20);
+
+  for (int i = 0; i < 20; i++)
+    info.elem (i) = 0;
+}
+
+DASPK::DASPK (const ColumnVector& state, const ColumnVector& deriv,
+	  double time, DAEFunc& f)
+  : DAE (state, deriv, time, f)
+{
+  n = size ();
+
+  stop_time_set = 0;
+  stop_time = 0.0;
+
+  DAEFunc::set_function (f.function ());
+  DAEFunc::set_jacobian_function (f.jacobian_function ());
+
+  sanity_checked = 0;
+
+  info.resize (20);
+
+  for (int i = 0; i < 20; i++)
+    info.elem (i) = 0;
+}
+
+void
+DASPK::force_restart (void)
+{
+  restart = 1;
+  integration_error = 0;
+}
+
+void
+DASPK::set_stop_time (double tt)
+{
+  stop_time_set = 1;
+  stop_time = tt;
+}
+
+void
+DASPK::clear_stop_time (void)
+{
+  stop_time_set = 0;
+}
+
+int
+ddaspk_f (const double& time, const double *state, const double *deriv,
+	  const double&, double *delta, int& ires, double *, int *)
+{
+  ColumnVector tmp_deriv (nn);
+  ColumnVector tmp_state (nn);
+  ColumnVector tmp_delta (nn);
+
+  for (int i = 0; i < nn; i++)
+    {
+      tmp_deriv.elem (i) = deriv [i];
+      tmp_state.elem (i) = state [i];
+    }
+
+  tmp_delta = user_fun (tmp_state, tmp_deriv, time, ires);
+
+  if (ires >= 0)
+    {
+      if (tmp_delta.length () == 0)
+	ires = -2;
+      else
+	{
+	  for (int i = 0; i < nn; i++)
+	    delta [i] = tmp_delta.elem (i);
+	}
+    }
+
+  return 0;
+}
+
+//NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT,
+//C                          WP, IWP, B, EPLIN, IER, RPAR, IPAR)
+
+int
+ddaspk_psol (const int& neq, const double& time, const double *state,
+	     const double *deriv, const double *savr,
+	     const double& cj, const double *wght, double *wp,
+	     int *iwp, double *b, const double& eplin, int& ier,
+	     double *, int*)
+{
+  abort ();
+}
+
+int
+ddaspk_j (const double& time, const double *, const double *, double *pd,
+	  const double& cj, double *, int *)
+{
+  ColumnVector tmp_state (nn);
+  ColumnVector tmp_deriv (nn);
+
+  // XXX FIXME XXX
+
+  Matrix tmp_dfdxdot (nn, nn);
+  Matrix tmp_dfdx (nn, nn);
+
+  DAEFunc::DAEJac tmp_jac;
+  tmp_jac.dfdxdot = &tmp_dfdxdot;
+  tmp_jac.dfdx    = &tmp_dfdx;
+
+  tmp_jac = user_jac (tmp_state, tmp_deriv, time);
+
+  // Fix up the matrix of partial derivatives for daspk.
+
+  tmp_dfdx = tmp_dfdx +  cj * tmp_dfdxdot;
+
+  for (int j = 0; j < nn; j++)
+    for (int i = 0; i < nn; i++)
+      pd [nn * j + i] = tmp_dfdx.elem (i, j);
+
+  return 0;
+}
+
+ColumnVector
+DASPK::do_integrate (double tout)
+{
+  ColumnVector retval;
+
+  if (restart)
+    {
+      restart = 0;
+      info.elem (0) = 0;
+    }
+
+  liw = 40 + n;
+  if (info(9) == 1 || info(9) == 3)
+    liw += n;
+  if (info (10) == 1 || info(15) == 1)
+    liw += n;
+
+  lrw = 50 + 9*n;
+  if (info(5) == 0)
+    lrw += n*n;
+  if (info(15) == 1)
+    lrw += n;
+
+  if (iwork.length () != liw)
+    iwork.resize (liw);
+
+  if (rwork.length () != lrw)
+    rwork.resize (lrw);
+
+  integration_error = 0;
+
+  if (DAEFunc::jacobian_function ())
+    info.elem (4) = 1;
+  else
+    info.elem (4) = 0;
+
+  double *px    = x.fortran_vec ();
+  double *pxdot = xdot.fortran_vec ();
+
+  nn = n;
+  user_fun = DAEFunc::fun;
+  user_jac = DAEFunc::jac;
+
+  if (! sanity_checked)
+    {
+      int ires = 0;
+
+      ColumnVector res = (*user_fun) (x, xdot, t, ires);
+
+      if (res.length () != x.length ())
+	{
+	  (*current_liboctave_error_handler)
+	    ("daspk: inconsistent sizes for state and residual vectors");
+
+	  integration_error = 1;
+	  return retval;
+	}
+
+      sanity_checked = 1;
+    }
+  
+  if (stop_time_set)
+    {
+      rwork.elem (0) = stop_time;
+      info.elem (3) = 1;
+    }
+  else
+    info.elem (3) = 0;
+
+  double abs_tol = absolute_tolerance ();
+  double rel_tol = relative_tolerance ();
+
+  if (initial_step_size () >= 0.0)
+    {
+      rwork.elem (2) = initial_step_size ();
+      info.elem (7) = 1;
+    }
+  else
+    info.elem (7) = 0;
+
+  if (maximum_step_size () >= 0.0)
+    {
+      rwork.elem (1) = maximum_step_size ();
+      info.elem (6) = 1;
+    }
+  else
+    info.elem (6) = 0;
+
+  double *dummy = 0;
+  int *idummy = 0;
+
+  int *pinfo = info.fortran_vec ();
+  int *piwork = iwork.fortran_vec ();
+  double *prwork = rwork.fortran_vec ();
+
+// again:
+
+  F77_XFCN (ddaspk, DDASPK, (ddaspk_f, n, t, px, pxdot, tout, pinfo,
+			     rel_tol, abs_tol, idid, prwork, lrw,
+			     piwork, liw, dummy, idummy, ddaspk_j,
+			     ddaspk_psol));
+
+  if (f77_exception_encountered)
+    {
+      integration_error = 1;
+      (*current_liboctave_error_handler) ("unrecoverable error in daspk");
+    }
+  else
+    {
+      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.
+
+	  retval = x;
+	  t = tout;
+	  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: // DDASPK 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: // DDASPK 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;
+	}
+    }
+
+  return retval;
+}
+
+Matrix
+DASPK::do_integrate (const ColumnVector& tout)
+{
+  Matrix dummy;
+  return integrate (tout, dummy);
+}
+
+Matrix
+DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out)
+{
+  Matrix retval;
+  int n_out = tout.capacity ();
+
+  if (n_out > 0 && n > 0)
+    {
+      retval.resize (n_out, n);
+      xdot_out.resize (n_out, n);
+
+      for (int i = 0; i < n; i++)
+	{
+	  retval.elem (0, i) = x.elem (i);
+	  xdot_out.elem (0, i) = xdot.elem (i);
+	}
+
+      for (int j = 1; j < n_out; j++)
+	{
+	  ColumnVector x_next = do_integrate (tout.elem (j));
+
+	  if (integration_error)
+	    return retval;
+
+	  for (int i = 0; i < n; i++)
+	    {
+	      retval.elem (j, i) = x_next.elem (i);
+	      xdot_out.elem (j, i) = xdot.elem (i);
+	    }
+	}
+    }
+
+  return retval;
+}
+
+Matrix
+DASPK::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
+{
+  Matrix dummy;
+  return integrate (tout, dummy, tcrit);
+}
+
+Matrix
+DASPK::integrate (const ColumnVector& tout, Matrix& xdot_out,
+		  const ColumnVector& tcrit) 
+{
+  Matrix retval;
+  int n_out = tout.capacity ();
+
+  if (n_out > 0 && n > 0)
+    {
+      retval.resize (n_out, n);
+      xdot_out.resize (n_out, n);
+
+      for (int i = 0; i < n; i++)
+	{
+	  retval.elem (0, i) = x.elem (i);
+	  xdot_out.elem (0, i) = xdot.elem (i);
+	}
+
+      int n_crit = tcrit.capacity ();
+
+      if (n_crit > 0)
+	{
+	  int i_crit = 0;
+	  int i_out = 1;
+	  double next_crit = tcrit.elem (0);
+	  double next_out;
+	  while (i_out < n_out)
+	    {
+	      bool do_restart = false;
+
+	      next_out = tout.elem (i_out);
+	      if (i_crit < n_crit)
+		next_crit = tcrit.elem (i_crit);
+
+	      bool save_output;
+	      double t_out;
+
+	      if (next_crit == next_out)
+		{
+		  set_stop_time (next_crit);
+		  t_out = next_out;
+		  save_output = true;
+		  i_out++;
+		  i_crit++;
+		  do_restart = true;
+		}
+	      else if (next_crit < next_out)
+		{
+		  if (i_crit < n_crit)
+		    {
+		      set_stop_time (next_crit);
+		      t_out = next_crit;
+		      save_output = false;
+		      i_crit++;
+		      do_restart = true;
+		    }
+		  else
+		    {
+		      clear_stop_time ();
+		      t_out = next_out;
+		      save_output = true;
+		      i_out++;
+		    }
+		}
+	      else
+		{
+		  set_stop_time (next_crit);
+		  t_out = next_out;
+		  save_output = true;
+		  i_out++;
+		}
+
+	      ColumnVector x_next = do_integrate (t_out);
+
+	      if (integration_error)
+		return retval;
+
+	      if (save_output)
+		{
+		  for (int i = 0; i < n; i++)
+		    {
+		      retval.elem (i_out-1, i) = x_next.elem (i);
+		      xdot_out.elem (i_out-1, i) = xdot.elem (i);
+		    }
+		}
+
+	      if (do_restart)
+		force_restart ();
+	    }
+	}
+      else
+	{
+	  retval = integrate (tout, xdot_out);
+
+	  if (integration_error)
+	    return retval;
+	}
+    }
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
new file mode 100644
--- /dev/null
+++ b/liboctave/DASPK.h
@@ -0,0 +1,168 @@
+/*
+
+Copyright (C) 1996, 1997, 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#if !defined (octave_DASPK_h)
+#define octave_DASPK_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+#include <cfloat>
+#include <cmath>
+
+#include "DAE.h"
+
+class
+DASPK_options
+{
+public:
+
+  DASPK_options (void) { init (); }
+
+  DASPK_options (const DASPK_options& opt) { copy (opt); }
+
+  DASPK_options& operator = (const DASPK_options& opt)
+    {
+      if (this != &opt)
+	copy (opt);
+
+      return *this;
+    }
+
+  ~DASPK_options (void) { }
+
+  void init (void)
+    {
+      x_absolute_tolerance = DBL_EPSILON * DBL_EPSILON;
+      x_initial_step_size = -1.0;
+      x_maximum_step_size = -1.0;
+      x_minimum_step_size = 0.0;
+      x_relative_tolerance = sqrt (DBL_EPSILON);
+    }
+
+  void copy (const DASPK_options& opt)
+    {
+      x_absolute_tolerance = opt.x_absolute_tolerance;
+      x_initial_step_size = opt.x_initial_step_size;
+      x_maximum_step_size = opt.x_maximum_step_size;
+      x_minimum_step_size = opt.x_minimum_step_size;
+      x_relative_tolerance = opt.x_relative_tolerance;
+    }
+
+  void set_default_options (void) { init (); }
+
+  void set_absolute_tolerance (double val)
+    { x_absolute_tolerance = (val > 0.0) ? val : DBL_EPSILON * DBL_EPSILON; }
+
+  void set_initial_step_size (double val)
+    { x_initial_step_size = (val >= 0.0) ? val : -1.0; }
+
+  void set_maximum_step_size (double val)
+    { x_maximum_step_size = (val >= 0.0) ? val : -1.0; }
+
+  void set_minimum_step_size (double val)
+    { x_minimum_step_size = (val >= 0.0) ? val : 0.0; }
+
+  void set_relative_tolerance (double val)
+    { x_relative_tolerance = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); }
+
+  double absolute_tolerance (void) { return x_absolute_tolerance; }
+
+  double initial_step_size (void) { return x_initial_step_size; }
+
+  double maximum_step_size (void) { return x_maximum_step_size; }
+
+  double minimum_step_size (void) { return x_minimum_step_size; }
+
+  double relative_tolerance (void) { return x_relative_tolerance; }
+
+private:
+
+  double x_absolute_tolerance;
+  double x_initial_step_size;
+  double x_maximum_step_size;
+  double x_minimum_step_size;
+  double x_relative_tolerance;
+};
+
+class
+DASPK : public DAE, public DASPK_options
+{
+public:
+
+  DASPK (void);
+
+  DASPK (const ColumnVector& x, double time, DAEFunc& f);
+
+  DASPK (const ColumnVector& x, const ColumnVector& xdot,
+	 double time, DAEFunc& f);
+
+  ~DASPK (void) { }
+
+  void force_restart (void);
+
+  void set_stop_time (double t);
+  void clear_stop_time (void);
+
+  ColumnVector do_integrate (double t);
+
+  Matrix do_integrate (const ColumnVector& tout);
+
+  Matrix do_integrate (const ColumnVector& tout, const ColumnVector& tcrit); 
+
+  Matrix integrate (const ColumnVector& tout, Matrix& xdot_out);
+
+  Matrix integrate (const ColumnVector& tout, Matrix& xdot_out,
+		    const ColumnVector& tcrit); 
+
+private:
+
+  double stop_time;
+  int stop_time_set;
+
+  int n;
+  int integration_error;
+  int restart;
+  int liw;  
+  int lrw;
+  int idid;
+  int sanity_checked;
+  Array<int> info;
+  Array<int> iwork;
+  Array<double> rwork;
+
+  friend int ddaspk_j (double *time, double *state, double *deriv,
+		       double *pd, double *cj, double *rpar, int *ipar);
+
+  friend int ddaspk_f (double *time, double *state, double *deriv,
+		       double *delta, int *ires, double *rpar, int *ipar);
+
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- a/liboctave/Makefile.in
+++ b/liboctave/Makefile.in
@@ -44,7 +44,7 @@
 	vx-rv-cs.h vx-s-ccv.h vx-s-crv.h \
 	vx-rv-crv.h vx-cv-ccv.h vx-crv-rv.h vx-ccv-cv.h
 
-INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.h DASSL.h FEGrid.h \
+INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.h DASPK.h DASSL.h FEGrid.h \
 	LinConst.h LP.h LPsolve.h LSODE.h NLConst.h NLEqn.h NLFunc.h \
 	NLP.h ODE.h ODEFunc.h Objective.h QP.h Quad.h Range.h base-de.h \
 	base-min.h byte-swap.h cmd-edit.h cmd-hist.h data-conv.h \
@@ -84,16 +84,14 @@
 	vx-rv-cs.cc vx-s-ccv.cc vx-s-crv.cc \
 	vx-rv-crv.cc vx-cv-ccv.cc vx-crv-rv.cc vx-ccv-cv.cc
 
-LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc DAE.cc DASSL.cc FEGrid.cc \
-	LinConst.cc LPsolve.cc LSODE.cc NLEqn.cc Quad.cc Range.cc \
-	data-conv.cc dir-ops.cc file-ops.cc \
-	file-stat.cc glob-match.cc \
-	idx-vector.cc lo-ieee.cc lo-mappers.cc lo-specfun.cc \
-	lo-sysdep.cc lo-utils.cc mach-info.cc oct-alloc.cc \
-	oct-env.cc oct-fftw.cc oct-group.cc \
-	oct-passwd.cc oct-shlib.cc \
-	oct-syscalls.cc oct-time.cc prog-args.cc \
-	str-vec.cc \
+LIBOCTAVE_CXX_SOURCES := Bounds.cc CollocWt.cc DAE.cc DASPK.cc \
+	DASSL.cc FEGrid.cc LinConst.cc LPsolve.cc LSODE.cc \
+	NLEqn.cc Quad.cc Range.cc data-conv.cc dir-ops.cc \
+	file-ops.cc file-stat.cc glob-match.cc idx-vector.cc \
+	lo-ieee.cc lo-mappers.cc lo-specfun.cc lo-sysdep.cc \
+	lo-utils.cc mach-info.cc oct-alloc.cc oct-env.cc \
+	oct-fftw.cc oct-group.cc oct-passwd.cc oct-shlib.cc \
+	oct-syscalls.cc oct-time.cc prog-args.cc str-vec.cc \
 	$(TEMPLATE_SRC) \
 	$(TI_SRC) \
 	$(MATRIX_SRC) \
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@
+2002-04-27  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* DLD-FUNCTIONS/daspk.cc: New file.
+	* Makefile.in (DLD_XSRC): Add it to the list.
+
 2002-04-25  Paul Kienzle <pkienzle@users.sf.net>
 
 	* DLD-FUNCTIONS/kron.cc: New file.
new file mode 100644
--- /dev/null
+++ b/src/DLD-FUNCTIONS/daspk.cc
@@ -0,0 +1,439 @@
+/*
+
+Copyright (C) 1996, 1997, 2002 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <string>
+
+#include <iomanip>
+#include <iostream>
+
+#include "DASPK.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "ov-fcn.h"
+#include "pager.h"
+#include "unwind-prot.h"
+#include "utils.h"
+#include "variables.h"
+
+// Global pointer for user defined function required by daspk.
+static octave_function *daspk_fcn;
+
+static DASPK_options daspk_opts;
+
+// Is this a recursive call?
+static int call_depth = 0;
+
+ColumnVector
+daspk_user_function (const ColumnVector& x, const ColumnVector& xdot,
+		     double t, int& ires)
+{
+  ColumnVector retval;
+
+  int nstates = x.capacity ();
+
+  assert (nstates == xdot.capacity ());
+
+  octave_value_list args;
+  args(2) = t;
+
+  if (nstates > 1)
+    {
+      Matrix m1 (nstates, 1);
+      Matrix m2 (nstates, 1);
+      for (int i = 0; i < nstates; i++)
+	{
+	  m1 (i, 0) = x (i);
+	  m2 (i, 0) = xdot (i);
+	}
+      octave_value state (m1);
+      octave_value deriv (m2);
+      args(1) = deriv;
+      args(0) = state;
+    }
+  else
+    {
+      double d1 = x (0);
+      double d2 = xdot (0);
+      octave_value state (d1);
+      octave_value deriv (d2);
+      args(1) = deriv;
+      args(0) = state;
+    }
+
+  if (daspk_fcn)
+    {
+      octave_value_list tmp = daspk_fcn->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  gripe_user_supplied_eval ("daspk");
+	  return retval;
+	}
+
+      int tlen = tmp.length ();
+      if (tlen > 0 && tmp(0).is_defined ())
+	{
+	  retval = ColumnVector (tmp(0).vector_value ());
+
+	  if (tlen > 1)
+	    ires = tmp(1).int_value ();
+
+	  if (error_state || retval.length () == 0)
+	    gripe_user_supplied_eval ("daspk");
+	}
+      else
+	gripe_user_supplied_eval ("daspk");
+    }
+
+  return retval;
+}
+
+#define DASPK_ABORT() \
+  do \
+    { \
+      unwind_protect::run_frame ("Fdaspk"); \
+      return retval; \
+    } \
+  while (0)
+
+#define DASPK_ABORT1(msg) \
+  do \
+    { \
+      ::error ("daspk: " msg); \
+      DASPK_ABORT (); \
+    } \
+  while (0)
+
+#define DASPK_ABORT2(fmt, arg) \
+  do \
+    { \
+      ::error ("daspk: " fmt, arg); \
+      DASPK_ABORT (); \
+    } \
+  while (0)
+
+DEFUN_DLD (daspk, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{x}, @var{xdot}] =} daspk (@var{fcn}, @var{x0}, @var{xdot0}, @var{t}, @var{t_crit})\n\
+Return a matrix of states and their first derivatives with respect to\n\
+@var{t}.  Each row in the result matrices correspond to one of the\n\
+elements in the vector @var{t}.  The first element of @var{t}\n\
+corresponds to the initial state @var{x0} and derivative @var{xdot0}, so\n\
+that the first row of the output @var{x} is @var{x0} and the first row\n\
+of the output @var{xdot} is @var{xdot0}.\n\
+\n\
+The first argument, @var{fcn}, is a string that names the function to\n\
+call to compute the vector of residuals for the set of equations.\n\
+It must have the form\n\
+\n\
+@example\n\
+@var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
+@end example\n\
+\n\
+@noindent\n\
+where @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
+scalar.\n\
+\n\
+The second and third arguments to @code{daspk} specify the initial\n\
+condition of the states and their derivatives, and the fourth argument\n\
+specifies a vector of output times at which the solution is desired, \n\
+including the time corresponding to the initial condition.\n\
+\n\
+The set of initial states and derivatives are not strictly required to\n\
+be consistent.  In practice, however, @sc{Dassl} is not very good at\n\
+determining a consistent set for you, so it is best if you ensure that\n\
+the initial values result in the function evaluating to zero.\n\
+\n\
+The fifth argument is optional, and may be used to specify a set of\n\
+times that the DAE solver should not integrate past.  It is useful for\n\
+avoiding difficulties with singularities and points where there is a\n\
+discontinuity in the derivative.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  unwind_protect::begin_frame ("Fdaspk");
+
+  unwind_protect_int (call_depth);
+  call_depth++;
+
+  if (call_depth > 1)
+    DASPK_ABORT1 ("invalid recursive call");
+
+  int nargin = args.length ();
+
+  if (nargin > 3 && nargin < 6)
+    {
+      daspk_fcn = extract_function
+	(args(0), "daspk", "__daspk_fcn__",
+	 "function res = __daspk_fcn__ (x, xdot, t) res = ",
+	 "; endfunction");
+
+      if (! daspk_fcn)
+	DASPK_ABORT ();
+
+      ColumnVector state = ColumnVector (args(1).vector_value ());
+
+      if (error_state)
+	DASPK_ABORT1 ("expecting state vector as second argument");
+
+      ColumnVector deriv (args(2).vector_value ());
+
+      if (error_state)
+	DASPK_ABORT1 ("expecting derivative vector as third argument");
+
+      ColumnVector out_times (args(3).vector_value ());
+
+      if (error_state)
+	DASPK_ABORT1 ("expecting output time vector as fourth argument");
+
+      ColumnVector crit_times;
+      int crit_times_set = 0;
+      if (nargin > 4)
+	{
+	  crit_times = ColumnVector (args(4).vector_value ());
+
+	  if (error_state)
+	    DASPK_ABORT1 ("expecting critical time vector as fifth argument");
+
+	  crit_times_set = 1;
+	}
+
+      if (state.capacity () != deriv.capacity ())
+	DASPK_ABORT1 ("x and xdot must have the same size");
+
+      double tzero = out_times (0);
+
+      DAEFunc func (daspk_user_function);
+      DASPK dae (state, deriv, tzero, func);
+      dae.copy (daspk_opts);
+
+      Matrix output;
+      Matrix deriv_output;
+
+      if (crit_times_set)
+	output = dae.integrate (out_times, deriv_output, crit_times);
+      else
+	output = dae.integrate (out_times, deriv_output);
+
+      if (! error_state)
+	{
+	  retval.resize (2);
+
+	  retval(0) = output;
+	  retval(1) = deriv_output;
+	}
+    }
+  else
+    print_usage ("daspk");
+
+  unwind_protect::run_frame ("Fdaspk");
+
+  return retval;
+}
+
+typedef void (DASPK_options::*d_set_opt_mf) (double);
+typedef double (DASPK_options::*d_get_opt_mf) (void);
+
+#define MAX_TOKENS 3
+
+struct DASPK_OPTIONS
+{
+  const char *keyword;
+  const char *kw_tok[MAX_TOKENS + 1];
+  int min_len[MAX_TOKENS + 1];
+  int min_toks_to_match;
+  d_set_opt_mf d_set_fcn;
+  d_get_opt_mf d_get_fcn;
+};
+
+static DASPK_OPTIONS daspk_option_table [] =
+{
+  { "absolute tolerance",
+    { "absolute", "tolerance", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    &DASPK_options::set_absolute_tolerance,
+    &DASPK_options::absolute_tolerance, },
+
+  { "initial step size",
+    { "initial", "step", "size", 0, },
+    { 1, 0, 0, 0, }, 1,
+    &DASPK_options::set_initial_step_size,
+    &DASPK_options::initial_step_size, },
+
+  { "maximum step size",
+    { "maximum", "step", "size", 0, },
+    { 2, 0, 0, 0, }, 1,
+    &DASPK_options::set_maximum_step_size,
+    &DASPK_options::maximum_step_size, },
+
+  { "relative tolerance",
+    { "relative", "tolerance", 0, 0, },
+    { 1, 0, 0, 0, }, 1,
+    &DASPK_options::set_relative_tolerance,
+    &DASPK_options::relative_tolerance, },
+
+  { 0,
+    { 0, 0, 0, 0, },
+    { 0, 0, 0, 0, }, 0,
+    0, 0, },
+};
+
+static void
+print_daspk_option_list (std::ostream& os)
+{
+  print_usage ("daspk_options", 1);
+
+  os << "\n"
+     << "Options for daspk include:\n\n"
+     << "  keyword                                  value\n"
+     << "  -------                                  -----\n\n";
+
+  DASPK_OPTIONS *list = daspk_option_table;
+
+  const char *keyword;
+  while ((keyword = list->keyword) != 0)
+    {
+      os << "  "
+	 << std::setiosflags (std::ios::left) << std::setw (40)
+	 << keyword
+	 << std::resetiosflags (std::ios::left)
+	 << " ";
+
+      double val = (daspk_opts.*list->d_get_fcn) ();
+      if (val < 0.0)
+	os << "computed automatically";
+      else
+	os << val;
+
+      os << "\n";
+      list++;
+    }
+
+  os << "\n";
+}
+
+static void
+set_daspk_option (const std::string& keyword, double val)
+{
+  DASPK_OPTIONS *list = daspk_option_table;
+
+  while (list->keyword != 0)
+    {
+      if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
+				list->min_toks_to_match, MAX_TOKENS))
+	{
+	  (daspk_opts.*list->d_set_fcn) (val);
+
+	  return;
+	}
+      list++;
+    }
+
+  warning ("daspk_options: no match for `%s'", keyword.c_str ());
+}
+
+static octave_value_list
+show_daspk_option (const std::string& keyword)
+{
+  octave_value retval;
+
+  DASPK_OPTIONS *list = daspk_option_table;
+
+  while (list->keyword != 0)
+    {
+      if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
+				list->min_toks_to_match, MAX_TOKENS))
+	{
+	  double val = (daspk_opts.*list->d_get_fcn) ();
+	  if (val < 0.0)
+	    retval = "computed automatically";
+	  else
+	    retval = val;
+
+	  return retval;
+	}
+      list++;
+    }
+
+  warning ("daspk_options: no match for `%s'", keyword.c_str ());
+
+  return retval;
+}
+
+DEFUN_DLD (daspk_options, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} daspk_options (@var{opt}, @var{val})\n\
+When called with two arguments, this function allows you set options\n\
+parameters for the function @code{lsode}.  Given one argument,\n\
+@code{daspk_options} returns the value of the corresponding option.  If\n\
+no arguments are supplied, the names of all the available options and\n\
+their current values are displayed.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+
+  if (nargin == 0)
+    {
+      print_daspk_option_list (octave_stdout);
+      return retval;
+    }
+  else if (nargin == 1 || nargin == 2)
+    {
+      std::string keyword = args(0).string_value ();
+
+      if (! error_state)
+	{
+	  if (nargin == 1)
+	    return show_daspk_option (keyword);
+	  else
+	    {
+	      double val = args(1).double_value ();
+
+	      if (! error_state)
+		{
+		  set_daspk_option (keyword, val);
+		  return retval;
+		}
+	    }
+	}
+    }
+
+  print_usage ("daspk_options");
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -39,8 +39,8 @@
   endif
 endif
 
-DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc dassl.cc \
-	det.cc eig.cc expm.cc fft.cc fft2.cc filter.cc find.cc \
+DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc daspk.cc \
+	dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc filter.cc find.cc \
 	fsolve.cc gammainc.cc getgrent.cc getpwent.cc getrusage.cc \
 	givens.cc hess.cc ifft.cc ifft2.cc inv.cc kron.cc log.cc lpsolve.cc \
 	lsode.cc lu.cc minmax.cc pinv.cc qr.cc quad.cc qz.cc rand.cc \