Mercurial > octave-libgccjit

changeset 3119:f3e1da120048

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1997-11-30 21:13:58 by jwe]
author jwe
date Sun, 30 Nov 1997 21:17:04 +0000
parents 74cc8e2fe2c0
children a2371c4a1d50
files libcruft/ChangeLog libcruft/specfun/ribesl.f libcruft/specfun/rjbesl.f libcruft/specfun/rkbesl.f libcruft/specfun/rybesl.f liboctave/ChangeLog liboctave/Makefile.in liboctave/dbleBessel.cc liboctave/dbleBessel.h src/ChangeLog src/DLD-FUNCTIONS/bessel.cc src/Makefile.in
diffstat 12 files changed, 526 insertions(+), 10 deletions(-) [+]
line wrap: on
line diff
--- a/libcruft/ChangeLog	Sat Nov 29 19:47:53 1997 +0000
+++ b/libcruft/ChangeLog	Sun Nov 30 21:17:04 1997 +0000
@@ -1,3 +1,11 @@
+Sat Nov 29 13:02:14 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* specfun/ribesl.f (ribesl): Use d1mach to get machine parameters.
+	SAVE static data between calls.  Use PARAMETERS where possible.
+	* specfun/rjbesl.f (rjbesl): Likewise.
+	* specfun/rkbesl.f (rkbesl): Likewise.
+	* specfun/rybesl.f (rybesl): Likewise.
+
 Fri Nov 28 14:05:12 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* specfun: New subdirectory.
--- a/libcruft/specfun/ribesl.f	Sat Nov 29 19:47:53 1997 +0000
+++ b/libcruft/specfun/ribesl.f	Sun Nov 30 21:17:04 1997 +0000
@@ -182,6 +182,8 @@
 C-------------------------------------------------------------------
       DATA FIRST /.TRUE./
 C-------------------------------------------------------------------
+      SAVE FIRST, NSIG, ENTEN, ENSIG, RTNSIG, ENMTEN, EXPARG, XLARGE
+C-------------------------------------------------------------------
 C  Statement functions for conversion
 C-------------------------------------------------------------------
       CONV(N) = DBLE(N)
--- a/libcruft/specfun/rjbesl.f	Sat Nov 29 19:47:53 1997 +0000
+++ b/libcruft/specfun/rjbesl.f	Sun Nov 30 21:17:04 1997 +0000
@@ -173,6 +173,8 @@
 C---------------------------------------------------------------------
       DATA FIRST /.TRUE./
 C---------------------------------------------------------------------
+      SAVE FACT, FIRST, NSIG, ENTEN, ENSIG, RTNSIG, ENMTEN, XLARGE
+C---------------------------------------------------------------------
 C Statement functions for conversion and the gamma function.
 C---------------------------------------------------------------------
       CONV(I) = DBLE(I)
--- a/libcruft/specfun/rkbesl.f	Sat Nov 29 19:47:53 1997 +0000
+++ b/libcruft/specfun/rkbesl.f	Sun Nov 30 21:17:04 1997 +0000
@@ -187,9 +187,13 @@
       DATA ESTF/4.18341D1, 7.1075D0, 6.4306D0, 4.25110D1, 1.35633D0,
      1          8.45096D1, 2.0D1/
 C---------------------------------------------------------------------
+      DATA FIRST /.TRUE./
+C---------------------------------------------------------------------
+      SAVE P, Q, R, S, T, ESTM, ESTF
+      SAVE FIRST, EPS, XINF, XMIN, SQXMIN, XMAX
+C---------------------------------------------------------------------
 C  Machine dependent parameters
 C---------------------------------------------------------------------
-      DATA FIRST /.TRUE./
       IF (FIRST) THEN
         EPS = D1MACH (4)
         XINF = D1MACH (2)
--- a/libcruft/specfun/rybesl.f	Sat Nov 29 19:47:53 1997 +0000
+++ b/libcruft/specfun/rybesl.f	Sun Nov 30 21:17:04 1997 +0000
@@ -176,9 +176,12 @@
      9        -0.76852840844786673690D-01,-0.28387654227602353814D+00,
      A         0.92187029365045265648D+00/
 C----------------------------------------------------------------------
+      DATA FIRST /.TRUE./
+C----------------------------------------------------------------------
+      SAVE CH, FIRST, EPS, XINF, XMIN, DEL, XLARGE, THRESH
+C----------------------------------------------------------------------
 C  Machine-dependent constants
 C----------------------------------------------------------------------
-      DATA FIRST /.TRUE./
       IF (FIRST) THEN
         EPS = D1MACH (4)
         XINF = D1MACH (2)
--- a/liboctave/ChangeLog	Sat Nov 29 19:47:53 1997 +0000
+++ b/liboctave/ChangeLog	Sun Nov 30 21:17:04 1997 +0000
@@ -1,3 +1,8 @@
+Sun Nov 30 14:59:12 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* dbleBessel.h, dbleBessel.cc: New files.
+	* Makefile.in (INCLUDES, SOURCES): Add them to the lists.
+
 Wed Nov 26 20:02:13 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* lo-sysdep.cc (octave_getcwd): Prefer getcwd over getwd.
--- a/liboctave/Makefile.in	Sat Nov 29 19:47:53 1997 +0000
+++ b/liboctave/Makefile.in	Sun Nov 30 21:17:04 1997 +0000
@@ -42,9 +42,9 @@
 
 INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.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 dir-ops.h file-ops.h file-stat.h getopt.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 \
+	dbleBessel.h dir-ops.h file-ops.h file-stat.h getopt.h \
 	glob-match.h idx-vector.h lo-ieee.h lo-mappers.h lo-sysdep.h \
 	lo-utils.h mach-info.h oct-alloc.h oct-cmplx.h oct-env.h \
 	oct-math.h oct-group.h oct-passwd.h oct-syscalls.h pathlen.h \
@@ -76,9 +76,9 @@
 	mx-m-cs.cc mx-m-dm.cc mx-s-cdm.cc mx-s-cm.cc mx-s-dm.cc
 
 SOURCES := Bounds.cc CollocWt.cc DAE.cc DASSL.cc FEGrid.cc \
-	LinConst.cc LPsolve.cc LSODE.cc NLEqn.cc \
-	Quad.cc Range.cc acosh.c asinh.c atanh.c cmd-edit.cc \
-	cmd-hist.cc data-conv.cc dir-ops.cc erf.c erfc.c f2c-main.c \
+	LinConst.cc LPsolve.cc LSODE.cc NLEqn.cc Quad.cc Range.cc \
+	acosh.c asinh.c atanh.c cmd-edit.cc cmd-hist.cc data-conv.cc \
+	dbleBessel.cc dir-ops.cc erf.c erfc.c f2c-main.c \
 	file-ops.cc file-stat.cc filemode.c gamma.c getopt.c getopt1.c \
 	glob-match.cc idx-vector.cc lgamma.c lo-ieee.cc lo-mappers.cc \
 	lo-sysdep.cc lo-utils.cc mach-info.cc mkdir.c oct-alloc.cc \
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/dbleBessel.cc	Sun Nov 30 21:17:04 1997 +0000
@@ -0,0 +1,221 @@
+/*
+
+Copyright (C) 1996 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 "Range.h"
+#include "dColVector.h"
+#include "dMatrix.h"
+#include "dbleBessel.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+#include "mx-inlines.cc"
+
+extern "C"
+{
+  int F77_FCN (rjbesl, RJBESL) (const double&, const double&,
+				const int&, double*, int&);
+
+  int F77_FCN (rybesl, RYBESL) (const double&, const double&,
+				const int&, double*, int&);
+
+  int F77_FCN (ribesl, RIBESL) (const double&, const double&,
+				const int&, const int&, double*, int&);
+
+  int F77_FCN (rkbesl, RKBESL) (const double&, const double&,
+				const int&, const int&, double*, int&);
+}
+
+int
+F77_FCN (ribesl, RIBESL) (const double& x, const double& alpha,
+			  const int& nb, double *result, int& ncalc)
+{
+  int retval = 0;
+  F77_FCN (ribesl, RIBESL) (x, alpha, nb, 1, result, ncalc);
+  return retval;
+}
+
+int
+F77_FCN (rkbesl, RKBESL) (const double& x, const double& alpha,
+			  const int& nb, double *result, int& ncalc)
+{
+  int retval = 0;
+  F77_FCN (rkbesl, RKBESL) (x, alpha, nb, 1, result, ncalc);
+  return retval;
+}
+
+typedef int (*fptr) (const double&, const double&, const int&, double*, int&);
+
+static Matrix
+do_bessel (fptr f, const char *fn, double alpha, const Matrix& x)
+{
+  Matrix retval;
+
+  if (alpha >= 0.0)
+    {
+      int nb = 1;
+
+      if (alpha >= 1.0)
+	{
+	  double integral;
+	  alpha = modf (alpha, &integral);
+	  nb = static_cast<int> (integral) + 1;
+	}
+
+      Array<double> tmp (nb);
+
+      double *ptmp = tmp.fortran_vec ();
+
+      int x_nr = x.rows ();
+      int x_nc = x.cols ();
+
+      retval.resize (x_nr, x_nc);
+
+      int ncalc;
+
+      for (int j = 0; j < x_nc; j++)
+	{
+	  for (int i = 0; i < x_nr; i++)
+	    {
+	      f (x(i,j), alpha, nb, ptmp, ncalc);
+
+	      // XXX FIXME XXX -- check ncalc.
+
+	      retval(i,j) = tmp(nb-1);
+	    }
+	}
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("%s: alpha must be greater than zero", fn);
+
+  return retval;
+}
+
+static Matrix
+do_bessel (fptr f, const char *fn, const Range& alpha_range,
+	   const ColumnVector& x)
+{
+  Matrix retval;
+
+  double alpha_base = alpha_range.base ();
+
+  if (alpha_base >= 0.0)
+    {
+      int nb_orig = alpha_range.nelem ();
+      int nb = nb_orig;
+      int offset = 0;
+
+      if (alpha_base >= 1.0)
+	{
+	  double integral;
+	  alpha_base = modf (alpha_base, &integral);
+	  offset = static_cast<int> (integral);
+	  nb += offset;
+	}
+
+      Array<double> tmp (nb);
+
+      double *ptmp = tmp.fortran_vec ();
+
+      int x_len = x.length ();
+
+      retval.resize (x_len, nb_orig);
+
+      int ncalc;
+
+      for (int i = 0; i < x_len; i++)
+	{
+	  f (x(i), alpha_base, nb, ptmp, ncalc);
+
+	  // XXX FIXME XXX -- check ncalc.
+
+	  for (int j = 0; j < nb_orig; j++)
+	    retval(i,j) = tmp(offset+j);
+	}
+    }
+  else
+    (*current_liboctave_error_handler)
+      ("%s: alpha must be greater than zero", fn);
+
+  return retval;
+}
+
+Matrix
+besselj (double alpha, const Matrix& x)
+{
+  return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x);
+}
+
+Matrix
+bessely (double alpha, const Matrix& x)
+{
+  return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x);
+}
+
+Matrix
+besseli (double alpha, const Matrix& x)
+{
+  return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x);
+}
+
+Matrix
+besselk (double alpha, const Matrix& x)
+{
+  return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x);
+}
+
+Matrix
+besselj (const Range& alpha, const ColumnVector& x)
+{
+  return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x);
+}
+
+Matrix
+bessely (const Range& alpha, const ColumnVector& x)
+{
+  return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x);
+}
+
+Matrix
+besseli (const Range& alpha, const ColumnVector& x)
+{
+  return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x);
+}
+
+Matrix
+besselk (const Range& alpha, const ColumnVector& x)
+{
+  return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x);
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/dbleBessel.h	Sun Nov 30 21:17:04 1997 +0000
@@ -0,0 +1,50 @@
+/*
+
+Copyright (C) 1996 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_BESSEL_h)
+#define octave_BESSEL_h 1
+
+#if defined (__GNUG__)
+#pragma interface
+#endif
+
+class ColumnVector;
+class Matrix;
+class Range;
+
+extern Matrix besselj (double alpha, const Matrix& x);
+extern Matrix bessely (double alpha, const Matrix& x);
+extern Matrix besseli (double alpha, const Matrix& x);
+extern Matrix besselk (double alpha, const Matrix& x);
+
+extern Matrix besselj (const Range& alpha, const ColumnVector& x);
+extern Matrix bessely (const Range& alpha, const ColumnVector& x);
+extern Matrix besseli (const Range& alpha, const ColumnVector& x);
+extern Matrix besselk (const Range& alpha, const ColumnVector& x);
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- a/src/ChangeLog	Sat Nov 29 19:47:53 1997 +0000
+++ b/src/ChangeLog	Sun Nov 30 21:17:04 1997 +0000
@@ -1,3 +1,8 @@
+Sun Nov 30 14:58:56 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* DLD-FUNCTIONS/bessel.cc: New file.
+	* Makefile.in (DLD_XSRC): Add it to the list.
+
 Thu Nov 27 23:28:59 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* lex.l (handle_string): Constructor for string class takes
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/bessel.cc	Sun Nov 30 21:17:04 1997 +0000
@@ -0,0 +1,216 @@
+/*
+
+Copyright (C) 1997 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 "dbleBessel.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "help.h"
+#include "oct-obj.h"
+#include "utils.h"
+
+#define DO_BESSEL(type, alpha, x) \
+  do \
+    { \
+      switch (type) \
+	{ \
+	  case 'j': \
+	    retval = besselj (alpha, x); \
+	    break; \
+ \
+	  case 'y': \
+	    retval = bessely (alpha, x); \
+	    break; \
+ \
+	  case 'i': \
+	    retval = besseli (alpha, x); \
+	    break; \
+ \
+	  case 'k': \
+	    retval = besselk (alpha, x); \
+	    break; \
+ \
+	  default: \
+	    break; \
+        } \
+    } \
+  while (0)
+
+static void
+gripe_bessel_arg_1 (const char *fn)
+{
+  error ("%s: alpha must be scalar or range with increment equal to 1", fn);
+}
+
+octave_value_list
+do_bessel (char type, const char *fn, const octave_value_list& args)
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  if (nargin == 2)
+    {
+      octave_value alpha_arg = args(0);
+
+      if (alpha_arg.is_scalar_type ())
+	{
+	  double alpha = alpha_arg.double_value ();
+
+	  if (! error_state)
+	    {
+	      Matrix x = args(1).matrix_value ();
+
+	      if (! error_state)
+		DO_BESSEL (type, alpha, x);
+	      else
+		error ("%s: expecting matrix as second argument", fn);
+	    }
+	  else
+	    gripe_bessel_arg_1 (fn);
+	}
+      else
+	{
+	  Range alpha;
+
+	  if (! alpha_arg.is_range ())
+	    {
+	      ColumnVector tmp = alpha_arg.vector_value ();
+
+	      if (! error_state)
+		{
+		  int len = tmp.length ();
+
+		  double base = tmp(0);
+
+		  for (int i = 1; i < len; i++)
+		    {
+		      if (tmp(i) != base + i)
+			{
+			  gripe_bessel_arg_1 (fn);
+			  break;
+			}
+		    }
+
+		  if (! error_state)
+		    alpha = Range (tmp(0), tmp(len-1));
+		}
+	    }
+	  else
+	    alpha = alpha_arg.range_value ();
+
+	  if (! error_state)
+	    {
+	      ColumnVector x = args(1).vector_value ();
+
+	      if (! error_state)
+		DO_BESSEL (type, alpha, x);
+	      else
+		error ("%s: expecting vector as second argument", fn);
+	    }
+	}
+    }
+  else
+    print_usage (fn);
+
+  return retval;
+}
+
+DEFUN_DLD (besselj, args, ,
+  "besselj (alpha, x)\n\
+\n\
+Compute Bessel functions of the first kind.\n\
+\n\
+X must be a real matrix, vector or scalar.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If ALPHA is a\n\
+range, X must be a vector or scalar, and the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.\n\
+\n\
+ALPHA must be greater than or equal to zero.  If ALPHA is a range, it\n\
+must have an increment equal to one.")
+{
+  return do_bessel ('j', "besselj", args);
+}
+
+DEFUN_DLD (bessely, args, ,
+  "bessely (alpha, x)\n\
+\n\
+Compute Bessel functions of the second kind.\n\
+\n\
+X must be a real matrix, vector or scalar.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If ALPHA is a\n\
+range, X must be a vector or scalar, and the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.\n\
+\n\
+ALPHA must be greater than or equal to zero.  If ALPHA is a range, it\n\
+must have an increment equal to one.")
+{
+  return do_bessel ('y', "bessely", args);
+}
+
+DEFUN_DLD (besseli, args, ,
+  "besseli (alpha, x)\n\
+\n\
+Compute modified Bessel functions of the first kind.\n\
+\n\
+X must be a real matrix, vector or scalar.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If ALPHA is a\n\
+range, X must be a vector or scalar, and the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.\n\
+\n\
+ALPHA must be greater than or equal to zero.  If ALPHA is a range, it\n\
+must have an increment equal to one.")
+{
+  return do_bessel ('i', "besseli", args);
+}
+
+DEFUN_DLD (besselk, args, ,
+  "besselk (alpha, x)\n\
+\n\
+Compute modified Bessel functions of the second kind.\n\
+\n\
+X must be a real matrix, vector or scalar.\n\
+\n\
+If ALPHA is a scalar, the result is the same size as X.  If ALPHA is a\n\
+range, X must be a vector or scalar, and the result is a matrix with\n\
+length(X) rows and length(ALPHA) columns.\n\
+\n\
+ALPHA must be greater than or equal to zero.  If ALPHA is a range, it\n\
+must have an increment equal to one.")
+{
+  return do_bessel ('k', "besselk", args);
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
+
--- a/src/Makefile.in	Sat Nov 29 19:47:53 1997 +0000
+++ b/src/Makefile.in	Sun Nov 30 21:17:04 1997 +0000
@@ -39,8 +39,8 @@
   endif
 endif
 
-DLD_XSRC := balance.cc chol.cc colloc.cc dassl.cc det.cc eig.cc \
-	expm.cc fft.cc fft2.cc filter.cc find.cc fsolve.cc \
+DLD_XSRC := balance.cc bessel.cc chol.cc colloc.cc dassl.cc det.cc \
+	eig.cc expm.cc fft.cc fft2.cc filter.cc find.cc fsolve.cc \
 	getgrent.cc getpwent.cc getrusage.cc givens.cc hess.cc \
 	ifft.cc ifft2.cc inv.cc log.cc lpsolve.cc lsode.cc \
 	lu.cc minmax.cc pinv.cc qr.cc quad.cc qzval.cc rand.cc \