* Bugs #5728 10229 invalid - Differential_equations: Quadpack update 87/11487/22
Paul BIGNIER [Tue, 14 May 2013 09:16:20 +0000 (11:16 +0200)]
'intg.f' calls Netlib's Quadpack.
The upstream is now matched.

Sole input: added a warning message in order to prevent edge missing.
This makes the bug invalid because it shows that it is too much to ask from the solver.

Updated tests and doc.

Change-Id: Ia92fe0b437717590b633af8691242e4aa77f62ac

29 files changed:
scilab/CHANGES_5.5.X
scilab/modules/differential_equations/Makefile.am
scilab/modules/differential_equations/Makefile.in
scilab/modules/differential_equations/help/en_US/intg.xml
scilab/modules/differential_equations/help/fr_FR/intg.xml
scilab/modules/differential_equations/help/ja_JP/intg.xml
scilab/modules/differential_equations/help/pt_BR/intg.xml
scilab/modules/differential_equations/help/ru_RU/intg.xml
scilab/modules/differential_equations/sci_gateway/fortran/intg.f
scilab/modules/differential_equations/src/c/DllmainDifferential_equations.c
scilab/modules/differential_equations/src/fortran/Elementary_functions_Import.def
scilab/modules/differential_equations/src/fortran/Output_stream_Import.def
scilab/modules/differential_equations/src/fortran/core_Import.def
scilab/modules/differential_equations/src/fortran/differential_equations_Import.def
scilab/modules/differential_equations/src/fortran/differential_equations_f.vfproj
scilab/modules/differential_equations/src/fortran/differential_equations_f2c.vcxproj
scilab/modules/differential_equations/src/fortran/differential_equations_f2c.vcxproj.filters
scilab/modules/differential_equations/src/fortran/dqag0.f [deleted file]
scilab/modules/differential_equations/src/fortran/dqags.f
scilab/modules/differential_equations/src/fortran/dqagse.f [new file with mode: 0644]
scilab/modules/differential_equations/src/fortran/dqelg.f [new file with mode: 0644]
scilab/modules/differential_equations/src/fortran/dqk21.f [new file with mode: 0644]
scilab/modules/differential_equations/src/fortran/dqpsrt.f [new file with mode: 0644]
scilab/modules/differential_equations/src/fortran/epsalg.f [deleted file]
scilab/modules/differential_equations/src/fortran/order.f [deleted file]
scilab/modules/differential_equations/src/fortran/quarul.f [deleted file]
scilab/modules/differential_equations/tests/unit_tests/integrate.dia.ref
scilab/modules/differential_equations/tests/unit_tests/intg.dia.ref
scilab/modules/differential_equations/tests/unit_tests/intg.tst

index 1d6ba4a..3292414 100644 (file)
@@ -78,6 +78,12 @@ Obsolete
   Use nicholschart() instead.
 
 
+Differential_equations
+======================
+
+* Netlib's Quadpack, used for definite integration, has been updated to match the upstream.
+
+
 Xcos
 ====
 
index f53c094..82435fe 100644 (file)
@@ -29,12 +29,12 @@ src/fortran/lsodar.f \
 src/fortran/ainvg.f \
 src/fortran/lsode.f \
 src/fortran/svcom1.f \
-src/fortran/quarul.f \
+src/fortran/dqk21.f \
 src/fortran/solsy.f \
 src/fortran/lsodi.f \
 src/fortran/ddassl.f \
-src/fortran/order.f \
-src/fortran/epsalg.f \
+src/fortran/dqpsrt.f \
+src/fortran/dqelg.f \
 src/fortran/cfode.f \
 src/fortran/rscma1.f \
 src/fortran/colnew.f \
@@ -43,10 +43,10 @@ src/fortran/xsetf.f \
 src/fortran/dgbsl.f \
 src/fortran/rkf45.f \
 src/fortran/rchek.f \
-src/fortran/dqag0.f \
+src/fortran/dqags.f \
 src/fortran/xerrwv.f \
 src/fortran/twodq.f \
-src/fortran/dqags.f \
+src/fortran/dqagse.f \
 src/fortran/greatr.f \
 src/fortran/hpdel.f \
 src/fortran/hpins.f \
index ad11141..f29b3fa 100644 (file)
@@ -148,10 +148,10 @@ am__objects_1 = libscidifferential_equations_algo_la-dassl.lo \
        libscidifferential_equations_algo_la-arnol.lo \
        libscidifferential_equations_algo_la-rk4.lo
 am__objects_2 = rscar1.lo bcomp.lo lcomp.lo loren.lo prja.lo vnorm.lo \
-       lsoda.lo lsodar.lo ainvg.lo lsode.lo svcom1.lo quarul.lo \
-       solsy.lo lsodi.lo ddassl.lo order.lo epsalg.lo cfode.lo \
+       lsoda.lo lsodar.lo ainvg.lo lsode.lo svcom1.lo dqk21.lo \
+       solsy.lo lsodi.lo ddassl.lo dqpsrt.lo dqelg.lo cfode.lo \
        rscma1.lo colnew.lo dcutet.lo xsetf.lo dgbsl.lo rkf45.lo \
-       rchek.lo dqag0.lo xerrwv.lo twodq.lo dqags.lo greatr.lo \
+       rchek.lo dqags.lo xerrwv.lo twodq.lo dqagse.lo greatr.lo \
        hpdel.lo hpins.lo svcar1.lo rscom1.lo rksimp.lo roots.lo \
        stoda.lo bnorm.lo rchek2.lo stode.lo vmnorm.lo prepj.lo \
        lsdisc.lo fnorm.lo ddaskr.lo daux.lo ddasrt.lo stodi.lo \
@@ -511,12 +511,12 @@ src/fortran/lsodar.f \
 src/fortran/ainvg.f \
 src/fortran/lsode.f \
 src/fortran/svcom1.f \
-src/fortran/quarul.f \
+src/fortran/dqk21.f \
 src/fortran/solsy.f \
 src/fortran/lsodi.f \
 src/fortran/ddassl.f \
-src/fortran/order.f \
-src/fortran/epsalg.f \
+src/fortran/dqpsrt.f \
+src/fortran/dqelg.f \
 src/fortran/cfode.f \
 src/fortran/rscma1.f \
 src/fortran/colnew.f \
@@ -525,10 +525,10 @@ src/fortran/xsetf.f \
 src/fortran/dgbsl.f \
 src/fortran/rkf45.f \
 src/fortran/rchek.f \
-src/fortran/dqag0.f \
+src/fortran/dqags.f \
 src/fortran/xerrwv.f \
 src/fortran/twodq.f \
-src/fortran/dqags.f \
+src/fortran/dqagse.f \
 src/fortran/greatr.f \
 src/fortran/hpdel.f \
 src/fortran/hpins.f \
@@ -1079,8 +1079,8 @@ lsode.lo: src/fortran/lsode.f
 svcom1.lo: src/fortran/svcom1.f
        $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o svcom1.lo `test -f 'src/fortran/svcom1.f' || echo '$(srcdir)/'`src/fortran/svcom1.f
 
-quarul.lo: src/fortran/quarul.f
-       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o quarul.lo `test -f 'src/fortran/quarul.f' || echo '$(srcdir)/'`src/fortran/quarul.f
+dqk21.lo: src/fortran/dqk21.f
+       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o dqk21.lo `test -f 'src/fortran/dqk21.f' || echo '$(srcdir)/'`src/fortran/dqk21.f
 
 solsy.lo: src/fortran/solsy.f
        $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o solsy.lo `test -f 'src/fortran/solsy.f' || echo '$(srcdir)/'`src/fortran/solsy.f
@@ -1091,11 +1091,11 @@ lsodi.lo: src/fortran/lsodi.f
 ddassl.lo: src/fortran/ddassl.f
        $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o ddassl.lo `test -f 'src/fortran/ddassl.f' || echo '$(srcdir)/'`src/fortran/ddassl.f
 
-order.lo: src/fortran/order.f
-       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o order.lo `test -f 'src/fortran/order.f' || echo '$(srcdir)/'`src/fortran/order.f
+dqpsrt.lo: src/fortran/dqpsrt.f
+       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o dqpsrt.lo `test -f 'src/fortran/dqpsrt.f' || echo '$(srcdir)/'`src/fortran/dqpsrt.f
 
-epsalg.lo: src/fortran/epsalg.f
-       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o epsalg.lo `test -f 'src/fortran/epsalg.f' || echo '$(srcdir)/'`src/fortran/epsalg.f
+dqelg.lo: src/fortran/dqelg.f
+       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o dqelg.lo `test -f 'src/fortran/dqelg.f' || echo '$(srcdir)/'`src/fortran/dqelg.f
 
 cfode.lo: src/fortran/cfode.f
        $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o cfode.lo `test -f 'src/fortran/cfode.f' || echo '$(srcdir)/'`src/fortran/cfode.f
@@ -1121,8 +1121,8 @@ rkf45.lo: src/fortran/rkf45.f
 rchek.lo: src/fortran/rchek.f
        $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o rchek.lo `test -f 'src/fortran/rchek.f' || echo '$(srcdir)/'`src/fortran/rchek.f
 
-dqag0.lo: src/fortran/dqag0.f
-       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o dqag0.lo `test -f 'src/fortran/dqag0.f' || echo '$(srcdir)/'`src/fortran/dqag0.f
+dqags.lo: src/fortran/dqags.f
+       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o dqags.lo `test -f 'src/fortran/dqags.f' || echo '$(srcdir)/'`src/fortran/dqags.f
 
 xerrwv.lo: src/fortran/xerrwv.f
        $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o xerrwv.lo `test -f 'src/fortran/xerrwv.f' || echo '$(srcdir)/'`src/fortran/xerrwv.f
@@ -1130,8 +1130,8 @@ xerrwv.lo: src/fortran/xerrwv.f
 twodq.lo: src/fortran/twodq.f
        $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o twodq.lo `test -f 'src/fortran/twodq.f' || echo '$(srcdir)/'`src/fortran/twodq.f
 
-dqags.lo: src/fortran/dqags.f
-       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o dqags.lo `test -f 'src/fortran/dqags.f' || echo '$(srcdir)/'`src/fortran/dqags.f
+dqagse.lo: src/fortran/dqagse.f
+       $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o dqagse.lo `test -f 'src/fortran/dqagse.f' || echo '$(srcdir)/'`src/fortran/dqagse.f
 
 greatr.lo: src/fortran/greatr.f
        $(LIBTOOL)  --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o greatr.lo `test -f 'src/fortran/greatr.f' || echo '$(srcdir)/'`src/fortran/greatr.f
index c89b570..c301d86 100644 (file)
@@ -18,7 +18,7 @@
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
-        <synopsis>[v,err]=intg(a,b,f [,ea [,er])</synopsis>
+        <synopsis>[v,err [,ierr]]=intg(a,b,f [,ea [,er])</synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Arguments</title>
                     <para>estimated absolute error on the result.</para>
                 </listitem>
             </varlistentry>
+            <varlistentry>
+                <term>ierr</term>
+                <listitem>
+                    <para>error flag number (= 0 if no error occured).</para>
+                </listitem>
+            </varlistentry>
         </variablelist>
     </refsection>
     <refsection>
         </para>
     </refsection>
     <refsection>
+        <title>Known Limitation</title>
+        <para>
+            Like all the integrators, <literal>intg</literal> is subject to spike missing.
+        </para>
+        <para>
+            A flat function with a spike will be seen as a fully flat function
+            if the spike is stiff enough.
+        </para>
+        <para>
+            This cannot be bypassed, it is easy to understand why when we know how the integrator operates.
+            Indeed, <literal>intg</literal> uses the 21-point Gauss-Kronrod rule,
+            so if there is a spike in-between two consecutive integration points,
+            then it will go undetected, the function will be considered smooth.
+        </para>
+        <para>
+            However, a warning message will be issued if the function is considered very smooth.
+            The user will then be suggested to reduce the integration interval,
+            should he think that spikes were missed.
+        </para>
+        <para>
+            The following graphs illustrate that phenomenon.
+        </para>
+        <scilab:image>
+            x = 0:.1:22;
+            y = zeros(1,221); y(1) = 3; y(96) = 1;
+            plot(x, y);
+            xtitle("Spike missed");
+        </scilab:image>
+        <para>
+            Being in-between the 9th and 10th integration points,
+            that spike is not detected and
+            <literal>intg</literal> considers the function flat.
+            In the next image, the spike is large enough to be covered by the integration points.
+        </para>
+        <scilab:image>
+            x = 0:21;
+            y = zeros(1,22); y(1) = 3; y(10) = 1;
+            plot(x, y);
+            xtitle("Spike detected");
+        </scilab:image>
+        <para>
+            If the user wants to display the computed solution even if the solver has encountered an error,
+            he should add the third output argument <literal>ierr</literal>, that will transform the
+            errors into warnings. This is mostly used in the case of rounding errors.
+        </para>
+    </refsection>
+    <refsection>
         <title>Examples</title>
         <programlisting role="example"><![CDATA[ 
 // Function written in the Scilab language
@@ -189,6 +242,6 @@ abs(exact-I)
         <para>The associated routines can be found in SCI/modules/differential_equations/src/fortran directory
             :
         </para>
-        <para>dqag0.f and dqags.f from quadpack</para>
+        <para>dqags.f and dqagse.f from quadpack</para>
     </refsection>
 </refentry>
index 93ed580..a493d4d 100644 (file)
@@ -6,7 +6,7 @@
     </refnamediv>
     <refsynopsisdiv>
         <title>Séquence d'appel</title>
-        <synopsis>[v,err]=intg(a,b,f [,ea [,er])</synopsis>
+        <synopsis>[v,err [,ierr]]=intg(a,b,f [,ea [,er])</synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Paramètres</title>
                     <para>estimation de l'erreur absolue sur le résultat.</para>
                 </listitem>
             </varlistentry>
+            <varlistentry>
+                <term>ierr</term>
+                <listitem>
+                    <para>numéro d'erreur (= 0 si absence d'erreur).</para>
+                </listitem>
+            </varlistentry>
         </variablelist>
     </refsection>
     <refsection>
         </para>
     </refsection>
     <refsection>
+        <title>Limitation connue</title>
+        <para>
+            Comme tous les intégrateurs, <literal>intg</literal> est susceptible de manquer des pics.
+        </para>
+        <para>
+            Une fonction plate contenant un pic sera considérée totalement plate si le pic est suffisamment raide.
+        </para>
+        <para>
+            Cela ne peut pas être évité, on comprend facilement pourquoi dès que l'on comprend comment l'intégrateur opère.
+            En effet, <literal>intg</literal> utilise la méthode des 21 points de Gauss-Kronrod,
+            donc s'il y a un pic entre deux points d'intégration consécutifs,
+            il ne sera pas détecté, la fonction sera considérée lisse.
+        </para>
+        <para>
+            Cependant, un warning s'affichera si la fonction est considérée très lisse,
+            qui suggèrera à l'utilisateur de réduire l'intervalle d'intégration,
+            s'il pense que des pics ont été manqués.
+        </para>
+        <para>
+            Les graphes suivants illustrent ce phénomène.
+        </para>
+        <scilab:image>
+            x = 0:.1:22;
+            y = zeros(1,221); y(1) = 3; y(96) = 1;
+            plot(x, y);
+            xtitle("Pic manqué");
+        </scilab:image>
+        <para>
+            Se situant entre les 9è and 10è points d'intégration, le pic n'est pas détecté et
+            <literal>intg</literal> considère la fonction comme plate.
+            Dans l'image suivante, le pic est suffisamment large pour être détecté.
+        </para>
+        <scilab:image>
+            x = 0:21;
+            y = zeros(1,22); y(1) = 3; y(10) = 1;
+            plot(x, y);
+            xtitle("Pic détecté");
+        </scilab:image>
+        <para>
+            Si l'utilisateur veut afficher la solution même si le solveur a rencontré une erreur,
+            il doit ajouter le troisième argument de sortie <literal>ierr</literal>, ce qui transformera
+            les erreurs en warnings. Ceci est surtout utilisé dans les cas d'erreurs d'arrondi.
+        </para>
+    </refsection>
+    <refsection>
         <title>Exemples</title>
         <programlisting role="example"><![CDATA[ 
 //External écrit en Scilab
@@ -177,7 +228,7 @@ abs(exact-I)
     </refsection>
     <refsection>
         <title>Fonctions Utilisées</title>
-        <para>Les programmes correspondants (dqag0.f et dqags.f de quadpack) se
+        <para>Les programmes correspondants (dqags.f et dqagse.f de quadpack) se
             trouvent dans le répertoire SCI/modules/differential_equations/src/fortran/.
         </para>
     </refsection>
index 22650be..4b55568 100644 (file)
@@ -195,6 +195,6 @@ abs(exact-I)
             SCI/modules/differential_equations/src/fortran ディレクトリにあります
             :
         </para>
-        <para>quadpackのdqag0.f および dqags.f</para>
+        <para>quadpackのdqags.f および dqagse.f</para>
     </refsection>
 </refentry>
index 3f3f701..d0982ed 100644 (file)
@@ -190,6 +190,6 @@ abs(exact-I)
         <para>As rotinas associadas podem ser encontradas no diretório
             SCI/modules/differential_equations/src/fortran :
         </para>
-        <para>dqag0.f e dqags.f de quadpack</para>
+        <para>dqags.f e dqagse.f de quadpack</para>
     </refsection>
 </refentry>
index 1722bec..22f7349 100644 (file)
@@ -206,6 +206,6 @@ abs(exact-I)
         <para>Связанные подпрограммы можно найти в директории
             SCI/modules/differential_equations/src/fortran:
         </para>
-        <para>dqag0.f и dqags.f из quadpack</para>
+        <para>dqags.f и dqagse.f из quadpack</para>
     </refsection>
 </refentry>
index a7bea46..a3c1cd6 100644 (file)
@@ -1,29 +1,29 @@
 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 c Copyright (C) INRIA
 c ...
-c 
+c
 c This file must be used under the terms of the CeCILL.
 c This source file is licensed as described in the file COPYING, which
 c you should have received as part of this distribution.  The terms
-c are also available at    
+c are also available at
 c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 c
       subroutine intg
 c     --------------------------------------------
-c     Scilab intg 
+c     Scilab intg
 c      implicit undefined (a-z)
       character*(4) fname
       character*6   namef
       include 'stack.h'
-      integer iero 
+      integer iero
       common/ierajf/iero
       common/cintg/namef
       external bintg,fintg
       double precision epsa,epsr,a,b,val,abserr
       logical getexternal, getscalar,type ,cremat
-      integer topk,lr,katop,kydot,top2,lra,lrb,lc
-      integer iipal,lpal,lw,liw,lpali,ifail
-      integer iadr,sadr
+      integer topk,lr,katop,kydot,top2,lra,lrb,lc,lier
+      integer iipal,lpal,lw,liw,lpali,ier
+      integer iadr,sadr,vfinite
       external setfintg
 c
       iadr(l)=l+l-1
@@ -33,6 +33,10 @@ c
          call error(39)
          return
       endif
+      if(lhs.gt.3) then
+         call error(59)
+         return
+      endif
       type=.false.
       top2=top
       topk=top
@@ -43,7 +47,7 @@ c
       else
          epsr=1.0d-8
       endif
-      if (rhs.ge.4) then 
+      if (rhs.ge.4) then
          if (.not.getscalar(fname,topk,top,lr)) return
          epsa=stk(lr)
          top=top-1
@@ -57,16 +61,26 @@ c     cas standard
       top=top-1
       if (.not.getscalar(fname,topk,top,lrb)) return
       b=stk(lrb)
+      if (isanan(b).eq.1.or.vfinite(1,b).eq.0) then
+         err = 2
+         call error(264)
+         return
+      endif
       top=top-1
       katop=top
       if (.not.getscalar(fname,topk,top,lra)) return
       a=stk(lra)
-c     tableaux de travail 
+      if (isanan(a).eq.1.or.vfinite(1,a).eq.0) then
+         err = 1
+         call error(264)
+         return
+      endif
+c     tableaux de travail
       top=top2+1
       lw=3000
       if (.not.cremat(fname,top,0,1,lw,lpal,lc)) return
       top=top+1
-c     tableau de travail entier necessaire 
+c     tableau de travail entier necessaire
       liw=3000/8+2
       if (.not.cremat(fname,top,0,1,iadr(liw)+1,lpali,lc)) return
       top=top+1
@@ -79,23 +93,85 @@ c
       istk(iipal+2)=kydot
       istk(iipal+3)=katop
       lstk(top+1)=sadr(iipal+4)
-      if(type) then 
-         call dqag0(fintg,a,b,epsa,epsr,val,abserr,
-     +        stk(lpal),lw,stk(lpali),liw,ifail)
+      if(type) then
+         call dqags(fintg,a,b,epsa,epsr,val,abserr,
+     +        neval,ier,lw/4,lw,last,stk(lpali),stk(lpal))
       else
-         call dqag0(bintg,a,b,epsa,epsr,val,abserr,
-     +        stk(lpal),lw,stk(lpali),liw,ifail)
+         call dqags(bintg,a,b,epsa,epsr,val,abserr,
+     +        neval,ier,lw/4,lw,last,stk(lpali),stk(lpal))
       endif
-      if(err.gt.0.or.err1.gt.0)return
-      if(ifail.gt.0) then
-         call error(24)
-         return
+      if(err.gt.0.or.err1.gt.0) return
+      if(ier.gt.0) then
+         if (lhs.le.2) then
+            if (ier.eq.1) then
+               call erro('Error: Maximum number of subdivisons '//
+     &            'achieved. Splitting the interval might help.')
+               return
+            elseif (ier.eq.2) then
+               call erro('Error: Round-off error detected, the '//
+     &         'requested tolerance (or default) cannot be achieved.'//
+     &         ' Try using bigger tolerances.')
+               return
+            elseif (ier.eq.3) then
+               call erro('Error: Bad integrand behavior occurs at '//
+     &            'some points of the integration interval.')
+               return
+            elseif (ier.eq.4) then
+               call erro('Error: Convergence problem, round-off '//
+     &            'error detected. Try using bigger tolerances.')
+               return
+            elseif (ier.eq.5) then
+               call erro('Error: The integral is probably '//
+     &            'divergent, or slowly convergent.')
+               return
+            else
+c           ier = 6
+               call erro('Error: Invalid input, absolute '//
+     &            'tolerance <= 0 and relative tolerance < 2.e-14.')
+               return
+            endif
+         else
+            if (ier.eq.1) then
+               call msgstxt('Warning: Maximum number of '//
+     &            'subdivisions achieved. '//
+     &            'Splitting the interval might help.')
+            elseif (ier.eq.2) then
+               call msgstxt('Warning: Round-off error detected, '//
+     &         'the requested tolerance (or default) cannot be '//
+     &         'achieved. Try using bigger tolerances.')
+            elseif (ier.eq.3) then
+               call msgstxt('Warning: Bad integrand behavior '//
+     &            'occurs at some points of the integration interval.')
+            elseif (ier.eq.4) then
+               call msgstxt('Warning: Convergence problem, round-off'//
+     &            'error detected. Try using bigger tolerances.')
+            elseif (ier.eq.5) then
+               call msgstxt('Warning: The integral is probably '//
+     &            'divergent, or slowly convergent.')
+             else
+c           ier = 6
+               call msgstxt('Warning: Invalid inputs, absolute '//
+     &            'tolerance <= 0 and relative tolerance < 2.e-14.')
+            endif
+        endif
       endif
       top=top2-rhs+1
       stk(lra)=val
-      if(lhs.eq.2) then
+      if(lhs.ge.2) then
          top=top+1
          stk(lrb)=abserr
+         if (lhs.eq.3) then
+            top=top+1
+            ilier=iadr(lstk(top))
+            istk(ilier)=1
+            istk(ilier+1)=1
+            istk(ilier+2)=1
+            istk(ilier+3)=0
+            lier=sadr(ilier+4)
+            stk(lier)=ier
+            lstk(top+1)=lier+1
+            return
+         endif
          return
       endif
       return
index 1e374e3..01df65c 100644 (file)
@@ -159,11 +159,6 @@ DIFFERENTIAL_EQUATIONS_IMPEXP struct
     int iero;
 } C2F(ierdcu);
 
-DIFFERENTIAL_EQUATIONS_IMPEXP struct
-{
-    int jupbnd;
-} C2F(dqa001);
-
 
 DIFFERENTIAL_EQUATIONS_IMPEXP struct
 {
index 5d1f540..b5ed092 100644 (file)
@@ -91,9 +91,9 @@
                <File RelativePath=".\ddasrt.f"/>
                <File RelativePath=".\ddassl.f"/>
                <File RelativePath=".\dgbsl.f"/>
-               <File RelativePath=".\dqag0.f"/>
                <File RelativePath=".\dqags.f"/>
-               <File RelativePath=".\epsalg.f"/>
+               <File RelativePath=".\dqagse.f"/>
+               <File RelativePath=".\dqelg.f"/>
                <File RelativePath=".\ewset.f"/>
                <File RelativePath="..\..\sci_gateway\fortran\Ex-bvode.f"/>
                <File RelativePath="..\..\sci_gateway\fortran\Ex-dasrt.f"/>
                <File RelativePath=".\lsodi.f"/>
                <File RelativePath=".\lsrgk.f"/>
                <File RelativePath=".\odeint.f"/>
-               <File RelativePath=".\order.f"/>
+               <File RelativePath=".\dqpsrt.f"/>
                <File RelativePath=".\prepj.f"/>
                <File RelativePath=".\prepji.f"/>
                <File RelativePath=".\prja.f"/>
-               <File RelativePath=".\quarul.f"/>
+               <File RelativePath=".\dqk21.f"/>
                <File RelativePath=".\rchek.f"/>
                <File RelativePath=".\rchek2.f"/>
                <File RelativePath=".\rkf45.f"/>
index f991fff..5f8f952 100644 (file)
@@ -309,9 +309,9 @@ cd ..
     <ClCompile Include="ddasrt.c" />
     <ClCompile Include="ddassl.c" />
     <ClCompile Include="dgbsl.c" />
-    <ClCompile Include="dqag0.c" />
     <ClCompile Include="dqags.c" />
-    <ClCompile Include="epsalg.c" />
+    <ClCompile Include="dqagse.c" />
+    <ClCompile Include="dqelg.c" />
     <ClCompile Include="ewset.c" />
     <ClCompile Include="..\..\sci_gateway\fortran\Ex-bvode.c" />
     <ClCompile Include="..\..\sci_gateway\fortran\Ex-dasrt.c" />
@@ -339,11 +339,11 @@ cd ..
     <ClCompile Include="lsodi.c" />
     <ClCompile Include="lsrgk.c" />
     <ClCompile Include="odeint.c" />
-    <ClCompile Include="order.c" />
+    <ClCompile Include="dqpsrt.c" />
     <ClCompile Include="prepj.c" />
     <ClCompile Include="prepji.c" />
     <ClCompile Include="prja.c" />
-    <ClCompile Include="quarul.c" />
+    <ClCompile Include="dqk21.c" />
     <ClCompile Include="rchek.c" />
     <ClCompile Include="rchek2.c" />
     <ClCompile Include="rkf45.c" />
@@ -406,9 +406,9 @@ cd ..
     <f2c_rule Include="ddasrt.f" />
     <f2c_rule Include="ddassl.f" />
     <f2c_rule Include="dgbsl.f" />
-    <f2c_rule Include="dqag0.f" />
     <f2c_rule Include="dqags.f" />
-    <f2c_rule Include="epsalg.f" />
+    <f2c_rule Include="dqagse.f" />
+    <f2c_rule Include="dqelg.f" />
     <f2c_rule Include="ewset.f" />
     <f2c_rule Include="..\..\sci_gateway\fortran\Ex-bvode.f" />
     <f2c_rule Include="..\..\sci_gateway\fortran\Ex-dasrt.f" />
@@ -436,11 +436,11 @@ cd ..
     <f2c_rule Include="lsodi.f" />
     <f2c_rule Include="lsrgk.f" />
     <f2c_rule Include="odeint.f" />
-    <f2c_rule Include="order.f" />
+    <f2c_rule Include="dqpsrt.f" />
     <f2c_rule Include="prepj.f" />
     <f2c_rule Include="prepji.f" />
     <f2c_rule Include="prja.f" />
-    <f2c_rule Include="quarul.f" />
+    <f2c_rule Include="dqk21.f" />
     <f2c_rule Include="rchek.f" />
     <f2c_rule Include="rchek2.f" />
     <f2c_rule Include="rkf45.f" />
@@ -489,4 +489,4 @@ cd ..
   <ImportGroup Label="ExtensionTargets">
     <Import Project="..\..\..\..\Visual-Studio-settings\f2c.targets" />
   </ImportGroup>
-</Project>
\ No newline at end of file
+</Project>
index 33970c5..540f288 100644 (file)
     <ClCompile Include="dgbsl.c">
       <Filter>Source Files</Filter>
     </ClCompile>
-    <ClCompile Include="dqag0.c">
+    <ClCompile Include="dqags.c">
       <Filter>Source Files</Filter>
     </ClCompile>
-    <ClCompile Include="dqags.c">
+    <ClCompile Include="dqagse.c">
       <Filter>Source Files</Filter>
     </ClCompile>
-    <ClCompile Include="epsalg.c">
+    <ClCompile Include="dqelg.c">
       <Filter>Source Files</Filter>
     </ClCompile>
     <ClCompile Include="ewset.c">
     <ClCompile Include="odeint.c">
       <Filter>Source Files</Filter>
     </ClCompile>
-    <ClCompile Include="order.c">
+    <ClCompile Include="dqpsrt.c">
       <Filter>Source Files</Filter>
     </ClCompile>
     <ClCompile Include="prepj.c">
     <ClCompile Include="prja.c">
       <Filter>Source Files</Filter>
     </ClCompile>
-    <ClCompile Include="quarul.c">
+    <ClCompile Include="dqk21.c">
       <Filter>Source Files</Filter>
     </ClCompile>
     <ClCompile Include="rchek.c">
     <f2c_rule Include="dgbsl.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
-    <f2c_rule Include="dqag0.f">
+    <f2c_rule Include="dqags.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
-    <f2c_rule Include="dqags.f">
+    <f2c_rule Include="dqagse.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
-    <f2c_rule Include="epsalg.f">
+    <f2c_rule Include="dqelg.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
     <f2c_rule Include="ewset.f">
     <f2c_rule Include="odeint.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
-    <f2c_rule Include="order.f">
+    <f2c_rule Include="dqpsrt.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
     <f2c_rule Include="prepj.f">
     <f2c_rule Include="prja.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
-    <f2c_rule Include="quarul.f">
+    <f2c_rule Include="dqk21.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
     <f2c_rule Include="rchek.f">
       <Filter>Libraries Dependencies</Filter>
     </None>
   </ItemGroup>
-</Project>
\ No newline at end of file
+</Project>
diff --git a/scilab/modules/differential_equations/src/fortran/dqag0.f b/scilab/modules/differential_equations/src/fortran/dqag0.f
deleted file mode 100644 (file)
index 1086ede..0000000
+++ /dev/null
@@ -1,51 +0,0 @@
-      subroutine dqag0(f, a, b, epsabs, epsrel, result, abserr,work,
-     * lwork, iwork, liwork, ifail)
-c
-c     based on quadpack routine
-c
-c     dqag0 is a general purpose integrator which calculates
-c     an approximation to the integral of a function over a finite
-c     interval (a,b)
-c
-c     dqag0 itself is essentially a dummy routine whose function is to
-c     partition the work arrays work and iwork for use by dqags.
-c     work is partitioned into 4 arrays each of size lwork/4.
-c     iwork is a single array in dqags.
-c
-c     .. scalar arguments ..
-      double precision a, abserr, b, epsabs, epsrel, result
-      integer ifail, liwork, lwork
-c     .. array arguments ..
-      double precision work(lwork)
-      integer iwork(liwork)
-c     .. function arguments ..
-      double precision f
-c     ..
-c     .. local scalars ..
-c
-      integer ibl, iel, ier, irl, limit
-c     .. function references ..
-c     .. subroutine references ..
-c     dqags
-c     ..
-      external f
-c     check that minimum workspace requirements are met
-      if (lwork.lt.4) go to 20
-      if (liwork.lt.lwork/8+2) go to 20
-c     limit = upper bound on number of subintervals
-      limit = lwork/4
-c     set up base addresses for work arrays
-      ibl = limit + 1
-      iel = limit + ibl
-      irl = limit + iel
-c     perform integration
-      call dqags(f, a, b, abs(epsabs), abs(epsrel), work(1),work(ibl)
-     *, work(iel), work(irl), limit, iwork, liwork,result, abserr, ier)
-      if (ier.ne.0) go to 40
-      ifail = 0
-      go to 60
-c     error 6 = insufficient workspace
-   20 ier = 6
-   40 ifail = 1
-   60 return
-      end
index 94248b3..ff4c086 100644 (file)
-      subroutine dqags(f, a, b, epsabs, epsrel, alist, blist,elist,
-     * rlist, limit, iord, liord, result, abserr, ier)
+      subroutine dqags(f,a,b,epsabs,epsrel,result,abserr,neval,ier,
+     *   limit,lenw,last,iwork,work)
+c***begin prologue  dqags
+c***date written   800101   (yymmdd)
+c***revision date  830518   (yymmdd)
+c***category no.  h2a1a1
+c***keywords  automatic integrator, general-purpose,
+c             (end-point) singularities, extrapolation,
+c             globally adaptive
+c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
+c           de doncker,elise,appl. math. & prog. div. - k.u.leuven
+c***purpose  the routine calculates an approximation result to a given
+c            definite integral  i = integral of f over (a,b),
+c            hopefully satisfying following claim for accuracy
+c            abs(i-result).le.max(epsabs,epsrel*abs(i)).
+c***description
+c
+c        computation of a definite integral
+c        standard fortran subroutine
+c        double precision version
 c
-c     based on quadpack routine dqags (formerly qags)
-c     **********************************************************
-c
-c        purpose
-c           the routine calculates an approximation
-c           /result/ to a given definite integral   i =
-c           integral of /f/ over (a,b), hopefully
-c           satisfying following claim for accuracy .
-c           abs(i-result) .le. max(epsabs,epsrel*abs(i)).
-c
-c          calling sequence
-c           call dqags (f,a,b,epsabs,epsrel,alist,blist,elist,
-c                        rlist,limit,iord,liord,result,abserr,ier)
 c
 c        parameters
-c            f      - function subprogram defining the integrand
-c                     function f(x). the actual name for f
-c                     needs to be declared e x t e r n a l
-c                     in the driver program
-c
-c            a      - lower limit of integration
-c
-c            b      - upper limit of integration
-c
-c            epsabs - absolute accuracy requested
-c
-c            epsrel - relative accuracy requested
-c
-c            alist,blist,elist,rlist
-c                   - work arrays (functions described below)
-c
-c            limit  - upper bound for number of subintervals
-c
-c            iord   - work array
-c
-c            liord  - length of iord (at least limit/2 + 2)
-c
-c            result - approximation to the integral
-c
-c            abserr - estimate of the modulus of the absolute error,
+c         on entry
+c            f      - double precision
+c                     function subprogram defining the integrand
+c                     function f(x). the actual name for f needs to be
+c                     declared e x t e r n a l in the driver program.
+c
+c            a      - double precision
+c                     lower limit of integration
+c
+c            b      - double precision
+c                     upper limit of integration
+c
+c            epsabs - double precision
+c                     absolute accuracy requested
+c            epsrel - double precision
+c                     relative accuracy requested
+c                     if  epsabs.le.0
+c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+c                     the routine will end with ier = 6.
+c
+c         on return
+c            result - double precision
+c                     approximation to the integral
+c
+c            abserr - double precision
+c                     estimate of the modulus of the absolute error,
 c                     which should equal or exceed abs(i-result)
 c
-c            ier    - ier   = 0 normal and reliable
-c                             termination of the routine.
-c                             it is assumed that the
-c                             requested  accuracy has been
-c                             achieved.
-c                   - ier   .ne. 0 abnormal termination of
-c                             the routine. the estimates
-c                             for integral and error are
-c                             less reliable. it is assumed
-c                             that the  requested accuracy
-c                             has not been achieved.
-c                         = 1 maximum number of subdivisions allowed
-c                             has been achieved. the user can
-c                             allow more sub divisions by
-c                             increasing the dimensions of the
-c                             work arrays work and iwork.
-c                             however, this may
-c                             yield no  improvement, and it
-c                             is rather advised to have a
-c                             close look at the integrand,
-c                             in order to determine the
-c                             integration  difficulties. if
-c                             the position of a local
-c                             difficulty can be determined
-c                             (i.e.  singularity,
-c                             discontinuity within the
-c                             interval) one will probably
-c                             gain from  splitting up the
-c                             interval at this point and
-c                             calling the integrator on the
-c                             sub-ranges. if possible, an
-c                             appropriate special-purpose
-c                             integrator should be used
-c                             which is designed for
-c                             handling the type  of
-c                             difficulty involved.
-c                         = 2 the occurrence of roundoff
-c                             error is detected which
-c                             prevents the requested
-c                             tolerance  from being
-c                             achieved. the error may be
-c                             under-estimated.
+c            neval  - integer
+c                     number of integrand evaluations
+c
+c            ier    - integer
+c                     ier = 0 normal and reliable termination of the
+c                             routine. it is assumed that the requested
+c                             accuracy has been achieved.
+c                     ier.gt.0 abnormal termination of the routine
+c                             the estimates for integral and error are
+c                             less reliable. it is assumed that the
+c                             requested accuracy has not been achieved.
+c            error messages
+c                     ier = 1 maximum number of subdivisions allowed
+c                             has been achieved. one can allow more sub-
+c                             divisions by increasing the value of limit
+c                             (and taking the according dimension
+c                             adjustments into account. however, if
+c                             this yields no improvement it is advised
+c                             to analyze the integrand in order to
+c                             determine the integration difficulties. if
+c                             the position of a local difficulty can be
+c                             determined (e.g. singularity,
+c                             discontinuity within the interval) one
+c                             will probably gain from splitting up the
+c                             interval at this point and calling the
+c                             integrator on the subranges. if possible,
+c                             an appropriate special-purpose integrator
+c                             should be used, which is designed for
+c                             handling the type of difficulty involved.
+c                         = 2 the occurrence of roundoff error is detec-
+c                             ted, which prevents the requested
+c                             tolerance from being achieved.
+c                             the error may be under-estimated.
 c                         = 3 extremely bad integrand behaviour
-c                             occurs at some interior points of the
-c                             integration interval.
-c                         = 4 it is presumed that the requested
-c                             tolerance cannot be achieved,
-c                             and that the returned result
-c                             is the best which can be
-c                             obtained.
+c                             occurs at some points of the integration
+c                             interval.
+c                         = 4 the algorithm does not converge.
+c                             roundoff error is detected in the
+c                             extrapolation table. it is presumed that
+c                             the requested tolerance cannot be
+c                             achieved, and that the returned result is
+c                             the best which can be obtained.
 c                         = 5 the integral is probably divergent, or
-c                             slowly convergent. it must be noted
-c                             that divergency can occur
-c                             with any other value of ier.
-c                         = -1 an error occurs during the evaluation of f
-c     **********************************************************
-c     .. scalar arguments ..
-      double precision a, abserr, b, epsabs, epsrel, result
-      integer ier, limit, liord
-c     .. array arguments ..
-      double precision alist(limit), blist(limit), elist(limit),
-     * rlist(limit)
-      integer iord(liord)
-c     .. function arguments ..
-      double precision f
-c     ..
-c     .. scalars in common ..
-      integer jupbnd
-      integer       iero
-      common/ierajf/iero
-c     ..
-c     .. local scalars ..
-      double precision a1, a2, abseps, area12, area1, area2, area, b1,
-     * b2,correc, defab1, defab2, defabs, dres, epmach, erlarg,erlast,
-     * errbnd, errmax, erro12, error1, error2, errsum,ertest, oflow,
-     * resabs, reseps, small, uflow
-      integer id, ierro, iroff1, iroff2, iroff3, k, ksgn, ktmin,last1,
-     * last, maxerr, nres, nrmax, numrl2
-      logical extrap, noext
-c     .. local arrays ..
-      double precision res3la(3), rlist2(52)
-c     .. function references ..
-      double precision dlamch
-c     .. subroutine references ..
-c     order, epsalg, quarul
-c     ..
-      external f
-      common /dqa001/ jupbnd
-c
-c            the dimension of /rlist2/ is determined by
-c            data /limexp/ in subroutine epsalg (/rlist2/
-c            should be of dimension (limexp+2) at least).
-c
-      epmach=dlamch('p')
-      uflow=dlamch('u')
-      oflow=dlamch('o')
-      iero=0
-c
-c            list of major variables
-c            -----------------------
-c
-c           alist     - list of left end-points of all subintervals
-c                       considered up to now
-c
-c           blist     - list of right end-points of all subintervals
-c                       considered up to now
-c
-c           rlist(i)  - approximation to the integral over
-c                       (alist(i),blist(i))
-c
-c           rlist2    - array of dimension at least limexp+2
-c                       containing the part of the epsilon table
-c                       which is still needed for further
-c                       computations
-c
-c           elist(i)  - error estimate applying to rlist(i)
-c
-c           maxerr    - pointer to the interval with largest error
-c                       estimate
-c
-c           errmax    - elist(maxerr)
-c
-c           erlast    - error on the interval currently subdivided
-c                       (before that subdivision has taken place)
-c
-c           area      - sum of the integrals over the subintervals
-c
-c           errsum    - sum of the errors over the subintervals
-c
-c           errbnd    - requested accuracy max(epsabs,epsrel*
-c                       abs(result))
-c
-c           *****1    - variable for the left interval
-c
-c           *****2    - variable for the right interval
-c
-c           last      - index for subdivision
-c
-c           nres      - number of calls to the extrapolation routine
-c
-c           numrl2    - number of elements currently  in
-c                       rlist2. if an appropriate
-c                       approximation to the compounded
-c                       integral has been obtained it is
-c                       put in  rlist2(numrl2) after numrl2
-c                       has been increased by one.
+c                             slowly convergent. it must be noted that
+c                             divergence can occur with any other value
+c                             of ier.
+c                         = 6 the input is invalid, because
+c                             (epsabs.le.0 and
+c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28)
+c                             or limit.lt.1 or lenw.lt.limit*4.
+c                             result, abserr, neval, last are set to
+c                             zero.except when limit or lenw is invalid,
+c                             iwork(1), work(limit*2+1) and
+c                             work(limit*3+1) are set to zero, work(1)
+c                             is set to a and work(limit+1) to b.
+c
+c         dimensioning parameters
+c            limit - integer
+c                    dimensioning parameter for iwork
+c                    limit determines the maximum number of subintervals
+c                    in the partition of the given integration interval
+c                    (a,b), limit.ge.1.
+c                    if limit.lt.1, the routine will end with ier = 6.
+c
+c            lenw  - integer
+c                    dimensioning parameter for work
+c                    lenw must be at least limit*4.
+c                    if lenw.lt.limit*4, the routine will end
+c                    with ier = 6.
+c
+c            last  - integer
+c                    on return, last equals the number of subintervals
+c                    produced in the subdivision process, detemines the
+c                    number of significant elements actually in the work
+c                    arrays.
+c
+c         work arrays
+c            iwork - integer
+c                    vector of dimension at least limit, the first k
+c                    elements of which contain pointers
+c                    to the error estimates over the subintervals
+c                    such that work(limit*3+iwork(1)),... ,
+c                    work(limit*3+iwork(k)) form a decreasing
+c                    sequence, with k = last if last.le.(limit/2+2),
+c                    and k = limit+1-last otherwise
+c
+c            work  - double precision
+c                    vector of dimension at least lenw
+c                    on return
+c                    work(1), ..., work(last) contain the left
+c                     end-points of the subintervals in the
+c                     partition of (a,b),
+c                    work(limit+1), ..., work(limit+last) contain
+c                     the right end-points,
+c                    work(limit*2+1), ..., work(limit*2+last) contain
+c                     the integral approximations over the subintervals,
+c                    work(limit*3+1), ..., work(limit*3+last)
+c                     contain the error estimates.
+c
+c***references  (none)
+c***routines called  dqagse,xerror
+c***end prologue  dqags
+c
+c
+      double precision a,abserr,b,epsabs,epsrel,f,result,work
+      integer ier,iwork,last,lenw,limit,lvl,l1,l2,l3,neval
+c
+      dimension iwork(limit),work(lenw)
 c
-c           small     - length of the smallest interval considered
-c                       up to now, multiplied by 1.5
-c
-c           erlarg    - sum of the errors over the intervals larger
-c                       than the smallest interval
-c                       considered up to now
-c           extrap    - logical variable denoting that the
-c                       routine is attempting to perform
-c                       extrapolation.  i.e. before
-c                       subdividing the smallest interval
-c                       we try to decrease the value of
-c                       erlarg
-c           noext     - logical variable denoting that extrapolation
-c                       is no longer allowed(/true/ value)
-c
-c           first approximation to the integral
-c           -----------------------------------
-c
-      last1 = 1
-      ier = 0
-      ierro = 0
-      call quarul(f, a, b, result, abserr, defabs, resabs)
-      if(iero.gt.0) then
-         ier=6
-         return
-      endif
-c
-c           test on accuracy
-c
-      dres = abs(result)
-      errbnd = max(epsabs,epsrel*dres)
-      if (abserr.le.1.0d+02*epmach*defabs .and. abserr.gt.errbnd)ier = 2
-      if (limit.lt.2 .and. abserr.gt.errbnd) ier = 1
-      if (ier.ne.0 .or. abserr.le.errbnd) go to 320
-c
-c           initialization
-c           --------------
-c
-      alist(1) = a
-      blist(1) = b
-      rlist(1) = result
-      rlist2(1) = result
-      errmax = abserr
-      maxerr = 1
-      area = result
-      errsum = abserr
-      abserr = oflow
-      nrmax = 1
-      nres = 0
-      numrl2 = 2
-      ktmin = 0
-      extrap = .false.
-      noext = .false.
-      iroff1 = 0
-      iroff2 = 0
-      iroff3 = 0
-      ksgn = -1
-      if (dres.ge.(0.10d+01-0.50d+02*epmach)*defabs) ksgn = 1
-c
-c           main do-loop
-c           ------------
-c
-      if (limit.lt.2) go to 220
-      do 200 last=2,limit
-c
-c           bisect the subinterval with the nrmax-th largest
-c           error estimate
-c
-         last1 = last
-         a1 = alist(maxerr)
-         b1 = 0.50d+00*(alist(maxerr)+blist(maxerr))
-         a2 = b1
-         b2 = blist(maxerr)
-         erlast = errmax
-         call quarul(f, a1, b1, area1, error1, resabs, defab1)
-         if(iero.gt.0) then
-            ier=6
-            return
-         endif
-         call quarul(f, a2, b2, area2, error2, resabs, defab2)
-         if(iero.gt.0) then
-            ier=6
-            return
-         endif
-c
-c           improve previous approximation of integral
-c           and error and test for accuracy
-c
-         area12 = area1 + area2
-         erro12 = error1 + error2
-         errsum = errsum + erro12 - errmax
-         area = area + area12 - rlist(maxerr)
-         if (defab1.eq.error1 .or. defab2.eq.error2) go to 40
-         if (abs(rlist(maxerr)-area12).gt.0.10d-04*abs(area12) .or.
-     *   erro12.lt.0.990d+00*errmax) go to 20
-         if (extrap) iroff2 = iroff2 + 1
-         if (.not.extrap) iroff1 = iroff1 + 1
-   20    if (last.gt.10 .and. erro12.gt.errmax) iroff3 = iroff3 + 1
-   40    rlist(maxerr) = area1
-         rlist(last) = area2
-         errbnd = max(epsabs,epsrel*abs(area))
-         if (errsum.le.errbnd) go to 280
-c
-c           test for roundoff error and eventually
-c           set error flag
-c
-         if (iroff1+iroff2.ge.10 .or. iroff3.ge.20) ier = 2
-         if (iroff2.ge.5) ierro = 3
-c
-c           set error flag in the case that the number of interval
-c            bisections exceeds /limit/
-c
-         if (last.eq.limit) ier = 1
-c
-c           set error flag in the case of bad integrand behaviour
-c           at interior points of integration range
-c
-         if (max(abs(a1),abs(b2)).le.(0.10d+01+0.10d+03*epmach)*
-     *   (abs(a2)+0.10d+04*uflow)) ier = 4
-         if (ier.ne.0) go to 220
-c
-c           append the newly-created intervals to the list
-c
-         if (error2.gt.error1) go to 60
-         alist(last) = a2
-         blist(maxerr) = b1
-         blist(last) = b2
-         elist(maxerr) = error1
-         elist(last) = error2
-         go to 80
-   60    alist(maxerr) = a2
-         alist(last) = a1
-         blist(last) = b1
-         rlist(maxerr) = area2
-         rlist(last) = area1
-         elist(maxerr) = error2
-         elist(last) = error1
-c
-c           call subroutine order to maintain the
-c           descending ordering in the list of error
-c           estimates and select the subinterval with
-c           nrmax-th largest error estimate (to be bisected
-c           next)
-c
-   80    call order(limit, last, maxerr, errmax, elist, iord,liord,
-     *    nrmax)
-         if (last.eq.2) go to 180
-         if (noext) go to 200
-         erlarg = erlarg - erlast
-         if (abs(b1-a1).gt.small) erlarg = erlarg + erro12
-         if (extrap) go to 100
-c
-c           test whether the interval to be bisected next is the
-c           smallest interval
-c
-         if (abs(blist(maxerr)-alist(maxerr)).gt.small) go to 200
-         extrap = .true.
-         nrmax = 2
-  100    if (ierro.eq.3 .or. erlarg.le.ertest) go to 140
-c
-c           the smallest interval has the largest error.
-c           before bisecting decrease the sum of the errors
-c           over the larger intervals(erlarg) and perform
-c           extrapolation
-c
-         id = nrmax
-         do 120 k=id,jupbnd
-            maxerr = iord(nrmax)
-            errmax = elist(maxerr)
-            if (abs(blist(maxerr)-alist(maxerr)).gt.small) go to 200
-            nrmax = nrmax + 1
-  120    continue
-c
-c           perform extrapolation
-c
-  140    numrl2 = numrl2 + 1
-         rlist2(numrl2) = area
-         call epsalg(numrl2, rlist2, reseps, abseps, res3la, nres)
-         ktmin = ktmin + 1
-         if (ktmin.gt.5 .and. abserr.lt.0.10d-02*errsum) ier = 5
-         if (abseps.ge.abserr) go to 160
-         ktmin = 0
-         abserr = abseps
-         result = reseps
-         correc = erlarg
-         ertest = max(epsabs,epsrel*abs(reseps))
-         if (abserr.le.ertest) go to 220
-c
-c           prepare  bisection of the smallest interval
+      external f
 c
-  160    if (numrl2.eq.1) noext = .true.
-         if (ier.eq.5) go to 220
-         maxerr = iord(1)
-         errmax = elist(maxerr)
-         nrmax = 1
-         extrap = .false.
-         small = small*0.50d+00
-         erlarg = errsum
-         go to 200
-  180    small = abs(b-a)*0.3750d+00
-         erlarg = errsum
-         ertest = errbnd
-         rlist2(2) = area
-  200 continue
+c         check validity of limit and lenw.
 c
-c           set  final result and error estimate
-c           ------------------------------------
+c***first executable statement  dqags
+      ier = 6
+      neval = 0
+      last = 0
+      result = 0.0d+00
+      abserr = 0.0d+00
+      if(limit.lt.1.or.lenw.lt.limit*4) go to 10
 c
-  220 if (abserr.eq.oflow) go to 280
-      if (ier+ierro.eq.0) go to 260
-      if (ierro.eq.3) abserr = abserr + correc
-      if (ier.eq.0) ier = 3
-      if (result.ne.0.0d+00.and .area. ne.0.0d+00) go to 240
-      if (abserr.gt.errsum) go to 280
-      if (area.eq.0.0d+00) go to 320
-      go to 260
-  240 if (abserr/abs(result).gt.errsum/abs(area)) go to 280
+c         prepare call for dqagse.
 c
-c           test on divergency
+      l1 = limit+1
+      l2 = limit+l1
+      l3 = limit+l2
 c
-  260 if (ksgn.eq.-1 .and. max(abs(result),abs(area)).le.defabs*
-     *0.10d-01) go to 320
-      if (0.10d-01.gt.(result/area) .or. (result/area).gt.0.10d+03.or.
-     * errsum.gt.abs(area)) ier = 6
-      go to 320
+      call dqagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval,
+     *  ier,work(1),work(l1),work(l2),work(l3),iwork,last)
 c
-c           compute global integral sum
+c         call error handler if necessary.
 c
-  280 result = 0.0d+00
-      do 300 k=1,last
-         result = result + rlist(k)
-  300 continue
-      abserr = errsum
-  320 if (ier.gt.2) ier = ier - 1
-      iord(1) = 4*last1
-      return
+c     Scilab input: comment next lines to not call xerror here,
+c     in order to let intg.f do the error handling, using ier.
+c      lvl = 0
+c10    if(ier.eq.6) lvl = 1
+c      if(ier.ne.0) call xerror(26habnormal return from dqags,26,ier,lvl)
+10      return
       end
diff --git a/scilab/modules/differential_equations/src/fortran/dqagse.f b/scilab/modules/differential_equations/src/fortran/dqagse.f
new file mode 100644 (file)
index 0000000..ca5362d
--- /dev/null
@@ -0,0 +1,467 @@
+      subroutine dqagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval,
+     *   ier,alist,blist,rlist,elist,iord,last)
+c***begin prologue  dqagse
+c***date written   800101   (yymmdd)
+c***revision date  830518   (yymmdd)
+c***category no.  h2a1a1
+c***keywords  automatic integrator, general-purpose,
+c             (end point) singularities, extrapolation,
+c             globally adaptive
+c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
+c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
+c***purpose  the routine calculates an approximation result to a given
+c            definite integral i = integral of f over (a,b),
+c            hopefully satisfying following claim for accuracy
+c            abs(i-result).le.max(epsabs,epsrel*abs(i)).
+c***description
+c
+c        computation of a definite integral
+c        standard fortran subroutine
+c        double precision version
+c
+c        parameters
+c         on entry
+c            f      - double precision
+c                     function subprogram defining the integrand
+c                     function f(x). the actual name for f needs to be
+c                     declared e x t e r n a l in the driver program.
+c
+c            a      - double precision
+c                     lower limit of integration
+c
+c            b      - double precision
+c                     upper limit of integration
+c
+c            epsabs - double precision
+c                     absolute accuracy requested
+c            epsrel - double precision
+c                     relative accuracy requested
+c                     if  epsabs.le.0
+c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+c                     the routine will end with ier = 6.
+c
+c            limit  - integer
+c                     gives an upperbound on the number of subintervals
+c                     in the partition of (a,b)
+c
+c         on return
+c            result - double precision
+c                     approximation to the integral
+c
+c            abserr - double precision
+c                     estimate of the modulus of the absolute error,
+c                     which should equal or exceed abs(i-result)
+c
+c            neval  - integer
+c                     number of integrand evaluations
+c
+c            ier    - integer
+c                     ier = 0 normal and reliable termination of the
+c                             routine. it is assumed that the requested
+c                             accuracy has been achieved.
+c                     ier.gt.0 abnormal termination of the routine
+c                             the estimates for integral and error are
+c                             less reliable. it is assumed that the
+c                             requested accuracy has not been achieved.
+c            error messages
+c                         = 1 maximum number of subdivisions allowed
+c                             has been achieved. one can allow more sub-
+c                             divisions by increasing the value of limit
+c                             (and taking the according dimension
+c                             adjustments into account). however, if
+c                             this yields no improvement it is advised
+c                             to analyze the integrand in order to
+c                             determine the integration difficulties. if
+c                             the position of a local difficulty can be
+c                             determined (e.g. singularity,
+c                             discontinuity within the interval) one
+c                             will probably gain from splitting up the
+c                             interval at this point and calling the
+c                             integrator on the subranges. if possible,
+c                             an appropriate special-purpose integrator
+c                             should be used, which is designed for
+c                             handling the type of difficulty involved.
+c                         = 2 the occurrence of roundoff error is detec-
+c                             ted, which prevents the requested
+c                             tolerance from being achieved.
+c                             the error may be under-estimated.
+c                         = 3 extremely bad integrand behaviour
+c                             occurs at some points of the integration
+c                             interval.
+c                         = 4 the algorithm does not converge.
+c                             roundoff error is detected in the
+c                             extrapolation table.
+c                             it is presumed that the requested
+c                             tolerance cannot be achieved, and that the
+c                             returned result is the best which can be
+c                             obtained.
+c                         = 5 the integral is probably divergent, or
+c                             slowly convergent. it must be noted that
+c                             divergence can occur with any other value
+c                             of ier.
+c                         = 6 the input is invalid, because
+c                             epsabs.le.0 and
+c                             epsrel.lt.max(50*rel.mach.acc.,0.5d-28).
+c                             result, abserr, neval, last, rlist(1),
+c                             iord(1) and elist(1) are set to zero.
+c                             alist(1) and blist(1) are set to a and b
+c                             respectively.
+c
+c            alist  - double precision
+c                     vector of dimension at least limit, the first
+c                      last  elements of which are the left end points
+c                     of the subintervals in the partition of the
+c                     given integration range (a,b)
+c
+c            blist  - double precision
+c                     vector of dimension at least limit, the first
+c                      last  elements of which are the right end points
+c                     of the subintervals in the partition of the given
+c                     integration range (a,b)
+c
+c            rlist  - double precision
+c                     vector of dimension at least limit, the first
+c                      last  elements of which are the integral
+c                     approximations on the subintervals
+c
+c            elist  - double precision
+c                     vector of dimension at least limit, the first
+c                      last  elements of which are the moduli of the
+c                     absolute error estimates on the subintervals
+c
+c            iord   - integer
+c                     vector of dimension at least limit, the first k
+c                     elements of which are pointers to the
+c                     error estimates over the subintervals,
+c                     such that elist(iord(1)), ..., elist(iord(k))
+c                     form a decreasing sequence, with k = last
+c                     if last.le.(limit/2+2), and k = limit+1-last
+c                     otherwise
+c
+c            last   - integer
+c                     number of subintervals actually produced in the
+c                     subdivision process
+c
+c***references  (none)
+c***routines called  d1mach,dqelg,dqk21,dqpsrt
+c***Scilab Enterprises input
+c        2013 - Paul Bignier
+c        internal parameter nbisect is the number of bisections carried out.
+c        if this parameter exceeds two, then a warning is issued to
+c        suggest the user to reduce the integration interval if he thinks
+c        edges have been missed.
+c***end prologue  dqagse
+c
+      double precision a,abseps,abserr,alist,area,area1,area12,area2,a1,
+     *  a2,b,blist,b1,b2,correc,dabs,defabs,defab1,defab2,d1mach,dmax1,
+     *  dres,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd,errmax,
+     *  error1,error2,erro12,errsum,ertest,f,oflow,resabs,reseps,result,
+     *  res3la,rlist,rlist2,small,uflow
+      integer id,ier,ierro,iord,iroff1,iroff2,iroff3,jupbnd,k,ksgn,
+     *  ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2
+c     /* Scilab Enterprises input
+      integer nbisect
+c     */
+      logical extrap,noext
+c
+      dimension alist(limit),blist(limit),elist(limit),iord(limit),
+     * res3la(3),rlist(limit),rlist2(52)
+c
+      external f
+c
+c            the dimension of rlist2 is determined by the value of
+c            limexp in subroutine dqelg (rlist2 should be of dimension
+c            (limexp+2) at least).
+c
+c            list of major variables
+c            -----------------------
+c
+c           alist     - list of left end points of all subintervals
+c                       considered up to now
+c           blist     - list of right end points of all subintervals
+c                       considered up to now
+c           rlist(i)  - approximation to the integral over
+c                       (alist(i),blist(i))
+c           rlist2    - array of dimension at least limexp+2 containing
+c                       the part of the epsilon table which is still
+c                       needed for further computations
+c           elist(i)  - error estimate applying to rlist(i)
+c           maxerr    - pointer to the interval with largest error
+c                       estimate
+c           errmax    - elist(maxerr)
+c           erlast    - error on the interval currently subdivided
+c                       (before that subdivision has taken place)
+c           area      - sum of the integrals over the subintervals
+c           errsum    - sum of the errors over the subintervals
+c           errbnd    - requested accuracy max(epsabs,epsrel*
+c                       abs(result))
+c           *****1    - variable for the left interval
+c           *****2    - variable for the right interval
+c           last      - index for subdivision
+c           nres      - number of calls to the extrapolation routine
+c           numrl2    - number of elements currently in rlist2. if an
+c                       appropriate approximation to the compounded
+c                       integral has been obtained it is put in
+c                       rlist2(numrl2) after numrl2 has been increased
+c                       by one.
+c           small     - length of the smallest interval considered up
+c                       to now, multiplied by 1.5
+c           erlarg    - sum of the errors over the intervals larger
+c                       than the smallest interval considered up to now
+c           extrap    - logical variable denoting that the routine is
+c                       attempting to perform extrapolation i.e. before
+c                       subdividing the smallest interval we try to
+c                       decrease the value of erlarg.
+c           noext     - logical variable denoting that extrapolation
+c                       is no longer allowed (true value)
+c
+c            machine dependent constants
+c            ---------------------------
+c
+c           epmach is the largest relative spacing.
+c           uflow is the smallest positive magnitude.
+c           oflow is the largest positive magnitude.
+c
+c***first executable statement  dqagse
+      epmach = d1mach(4)
+c
+c            test on validity of parameters
+c            ------------------------------
+      ier = 0
+      neval = 0
+      last = 0
+      result = 0.0d+00
+      abserr = 0.0d+00
+      alist(1) = a
+      blist(1) = b
+      rlist(1) = 0.0d+00
+      elist(1) = 0.0d+00
+      if(epsabs.le.0.0d+00.and.epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28))
+     *   ier = 6
+      if(ier.eq.6) go to 999
+c
+c           first approximation to the integral
+c           -----------------------------------
+c
+      uflow = d1mach(1)
+      oflow = d1mach(2)
+      ierro = 0
+c     /* Scilab Enterprises input
+      nbisect = 0
+c     */
+      call dqk21(f,a,b,result,abserr,defabs,resabs)
+c
+c           test on accuracy.
+c
+      dres = dabs(result)
+      errbnd = dmax1(epsabs,epsrel*dres)
+      last = 1
+      rlist(1) = result
+      elist(1) = abserr
+      iord(1) = 1
+      if(abserr.le.1.0d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2
+      if(limit.eq.1) ier = 1
+      if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or.
+     *  abserr.eq.0.0d+00) go to 140
+c
+c           initialization
+c           --------------
+c
+      rlist2(1) = result
+      errmax = abserr
+      maxerr = 1
+      area = result
+      errsum = abserr
+      abserr = oflow
+      nrmax = 1
+      nres = 0
+      numrl2 = 2
+      ktmin = 0
+      extrap = .false.
+      noext = .false.
+      iroff1 = 0
+      iroff2 = 0
+      iroff3 = 0
+      ksgn = -1
+      if(dres.ge.(0.1d+01-0.5d+02*epmach)*defabs) ksgn = 1
+c
+c           main do-loop
+c           ------------
+c
+      do 90 last = 2,limit
+c
+c           bisect the subinterval with the nrmax-th largest error
+c           estimate.
+c
+c       /* Scilab Enterprises input
+        nbisect = nbisect + 1
+c       */
+        a1 = alist(maxerr)
+        b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
+        a2 = b1
+        b2 = blist(maxerr)
+        erlast = errmax
+        call dqk21(f,a1,b1,area1,error1,resabs,defab1)
+        call dqk21(f,a2,b2,area2,error2,resabs,defab2)
+c
+c           improve previous approximations to integral
+c           and error and test for accuracy.
+c
+        area12 = area1+area2
+        erro12 = error1+error2
+        errsum = errsum+erro12-errmax
+        area = area+area12-rlist(maxerr)
+        if(defab1.eq.error1.or.defab2.eq.error2) go to 15
+        if(dabs(rlist(maxerr)-area12).gt.0.1d-04*dabs(area12)
+     *  .or.erro12.lt.0.99d+00*errmax) go to 10
+        if(extrap) iroff2 = iroff2+1
+        if(.not.extrap) iroff1 = iroff1+1
+   10   if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1
+   15   rlist(maxerr) = area1
+        rlist(last) = area2
+        errbnd = dmax1(epsabs,epsrel*dabs(area))
+c
+c           test for roundoff error and eventually set error flag.
+c
+        if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2
+        if(iroff2.ge.5) ierro = 3
+c
+c           set error flag in the case that the number of subintervals
+c           equals limit.
+c
+        if(last.eq.limit) ier = 1
+c
+c           set error flag in the case of bad integrand behaviour
+c           at a point of the integration range.
+c
+        if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)*
+     *  (dabs(a2)+0.1d+04*uflow)) ier = 4
+c
+c           append the newly-created intervals to the list.
+c
+        if(error2.gt.error1) go to 20
+        alist(last) = a2
+        blist(maxerr) = b1
+        blist(last) = b2
+        elist(maxerr) = error1
+        elist(last) = error2
+        go to 30
+   20   alist(maxerr) = a2
+        alist(last) = a1
+        blist(last) = b1
+        rlist(maxerr) = area2
+        rlist(last) = area1
+        elist(maxerr) = error2
+        elist(last) = error1
+c
+c           call subroutine dqpsrt to maintain the descending ordering
+c           in the list of error estimates and select the subinterval
+c           with nrmax-th largest error estimate (to be bisected next).
+c
+   30   call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
+c ***jump out of do-loop
+        if(errsum.le.errbnd) go to 115
+c ***jump out of do-loop
+        if(ier.ne.0) go to 100
+        if(last.eq.2) go to 80
+        if(noext) go to 90
+        erlarg = erlarg-erlast
+        if(dabs(b1-a1).gt.small) erlarg = erlarg+erro12
+        if(extrap) go to 40
+c
+c           test whether the interval to be bisected next is the
+c           smallest interval.
+c
+        if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90
+        extrap = .true.
+        nrmax = 2
+   40   if(ierro.eq.3.or.erlarg.le.ertest) go to 60
+c
+c           the smallest interval has the largest error.
+c           before bisecting decrease the sum of the errors over the
+c           larger intervals (erlarg) and perform extrapolation.
+c
+        id = nrmax
+        jupbnd = last
+        if(last.gt.(2+limit/2)) jupbnd = limit+3-last
+        do 50 k = id,jupbnd
+          maxerr = iord(nrmax)
+          errmax = elist(maxerr)
+c ***jump out of do-loop
+          if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90
+          nrmax = nrmax+1
+   50   continue
+c
+c           perform extrapolation.
+c
+   60   numrl2 = numrl2+1
+        rlist2(numrl2) = area
+        call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
+        ktmin = ktmin+1
+        if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5
+        if(abseps.ge.abserr) go to 70
+        ktmin = 0
+        abserr = abseps
+        result = reseps
+        correc = erlarg
+        ertest = dmax1(epsabs,epsrel*dabs(reseps))
+c ***jump out of do-loop
+        if(abserr.le.ertest) go to 100
+c
+c           prepare bisection of the smallest interval.
+c
+   70   if(numrl2.eq.1) noext = .true.
+        if(ier.eq.5) go to 100
+        maxerr = iord(1)
+        errmax = elist(maxerr)
+        nrmax = 1
+        extrap = .false.
+        small = small*0.5d+00
+        erlarg = errsum
+        go to 90
+   80   small = dabs(b-a)*0.375d+00
+        erlarg = errsum
+        ertest = errbnd
+        rlist2(2) = area
+   90 continue
+c
+c           set final result and error estimate.
+c           ------------------------------------
+c
+  100 if(abserr.eq.oflow) go to 115
+      if(ier+ierro.eq.0) go to 110
+      if(ierro.eq.3) abserr = abserr+correc
+      if(ier.eq.0) ier = 3
+      if(result.ne.0.0d+00.and.area.ne.0.0d+00) go to 105
+      if(abserr.gt.errsum) go to 115
+      if(area.eq.0.0d+00) go to 130
+      go to 110
+  105 if(abserr/dabs(result).gt.errsum/dabs(area)) go to 115
+c
+c           test on divergence.
+c
+  110 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le.
+     * defabs*0.1d-01) go to 130
+      if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03
+     * .or.errsum.gt.dabs(area)) ier = 6
+      go to 130
+c
+c           compute global integral sum.
+c
+  115 result = 0.0d+00
+      do 120 k = 1,last
+         result = result+rlist(k)
+  120 continue
+      abserr = errsum
+  130 if(ier.gt.2) ier = ier-1
+  140 neval = 42*last-21
+c     /* Scilab Enterprises input
+      if (nbisect.le.1) then
+         call msgstxt('Warning: argument function is detected to be '//
+     &      'very smooth.')
+         call msgstxt('Reduce integration interval if you think '//
+     &      'edges have been missed.')
+      endif
+c     */
+  999 return
+      end
diff --git a/scilab/modules/differential_equations/src/fortran/dqelg.f b/scilab/modules/differential_equations/src/fortran/dqelg.f
new file mode 100644 (file)
index 0000000..6266682
--- /dev/null
@@ -0,0 +1,184 @@
+      subroutine dqelg(n,epstab,result,abserr,res3la,nres)
+c***begin prologue  dqelg
+c***refer to  dqagie,dqagoe,dqagpe,dqagse
+c***routines called  d1mach
+c***revision date  830518   (yymmdd)
+c***keywords  epsilon algorithm, convergence acceleration,
+c             extrapolation
+c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
+c           de doncker,elise,appl. math & progr. div. - k.u.leuven
+c***purpose  the routine determines the limit of a given sequence of
+c            approximations, by means of the epsilon algorithm of
+c            p.wynn. an estimate of the absolute error is also given.
+c            the condensed epsilon table is computed. only those
+c            elements needed for the computation of the next diagonal
+c            are preserved.
+c***description
+c
+c           epsilon algorithm
+c           standard fortran subroutine
+c           double precision version
+c
+c           parameters
+c              n      - integer
+c                       epstab(n) contains the new element in the
+c                       first column of the epsilon table.
+c
+c              epstab - double precision
+c                       vector of dimension 52 containing the elements
+c                       of the two lower diagonals of the triangular
+c                       epsilon table. the elements are numbered
+c                       starting at the right-hand corner of the
+c                       triangle.
+c
+c              result - double precision
+c                       resulting approximation to the integral
+c
+c              abserr - double precision
+c                       estimate of the absolute error computed from
+c                       result and the 3 previous results
+c
+c              res3la - double precision
+c                       vector of dimension 3 containing the last 3
+c                       results
+c
+c              nres   - integer
+c                       number of calls to the routine
+c                       (should be zero at first call)
+c
+c***end prologue  dqelg
+c
+      double precision abserr,dabs,delta1,delta2,delta3,dmax1,d1mach,
+     *  epmach,epsinf,epstab,error,err1,err2,err3,e0,e1,e1abs,e2,e3,
+     *  oflow,res,result,res3la,ss,tol1,tol2,tol3
+      integer i,ib,ib2,ie,indx,k1,k2,k3,limexp,n,newelm,nres,num
+      dimension epstab(52),res3la(3)
+c
+c           list of major variables
+c           -----------------------
+c
+c           e0     - the 4 elements on which the computation of a new
+c           e1       element in the epsilon table is based
+c           e2
+c           e3                 e0
+c                        e3    e1    new
+c                              e2
+c           newelm - number of elements to be computed in the new
+c                    diagonal
+c           error  - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
+c           result - the element in the new diagonal with least value
+c                    of error
+c
+c           machine dependent constants
+c           ---------------------------
+c
+c           epmach is the largest relative spacing.
+c           oflow is the largest positive magnitude.
+c           limexp is the maximum number of elements the epsilon
+c           table can contain. if this number is reached, the upper
+c           diagonal of the epsilon table is deleted.
+c
+c***first executable statement  dqelg
+      epmach = d1mach(4)
+      oflow = d1mach(2)
+      nres = nres+1
+      abserr = oflow
+      result = epstab(n)
+      if(n.lt.3) go to 100
+      limexp = 50
+      epstab(n+2) = epstab(n)
+      newelm = (n-1)/2
+      epstab(n) = oflow
+      num = n
+      k1 = n
+      do 40 i = 1,newelm
+        k2 = k1-1
+        k3 = k1-2
+        res = epstab(k1+2)
+        e0 = epstab(k3)
+        e1 = epstab(k2)
+        e2 = res
+        e1abs = dabs(e1)
+        delta2 = e2-e1
+        err2 = dabs(delta2)
+        tol2 = dmax1(dabs(e2),e1abs)*epmach
+        delta3 = e1-e0
+        err3 = dabs(delta3)
+        tol3 = dmax1(e1abs,dabs(e0))*epmach
+        if(err2.gt.tol2.or.err3.gt.tol3) go to 10
+c
+c           if e0, e1 and e2 are equal to within machine
+c           accuracy, convergence is assumed.
+c           result = e2
+c           abserr = abs(e1-e0)+abs(e2-e1)
+c
+        result = res
+        abserr = err2+err3
+c ***jump out of do-loop
+        go to 100
+   10   e3 = epstab(k1)
+        epstab(k1) = e1
+        delta1 = e1-e3
+        err1 = dabs(delta1)
+        tol1 = dmax1(e1abs,dabs(e3))*epmach
+c
+c           if two elements are very close to each other, omit
+c           a part of the table by adjusting the value of n
+c
+        if(err1.le.tol1.or.err2.le.tol2.or.err3.le.tol3) go to 20
+        ss = 0.1d+01/delta1+0.1d+01/delta2-0.1d+01/delta3
+        epsinf = dabs(ss*e1)
+c
+c           test to detect irregular behaviour in the table, and
+c           eventually omit a part of the table adjusting the value
+c           of n.
+c
+        if(epsinf.gt.0.1d-03) go to 30
+   20   n = i+i-1
+c ***jump out of do-loop
+        go to 50
+c
+c           compute a new element and eventually adjust
+c           the value of result.
+c
+   30   res = e1+0.1d+01/ss
+        epstab(k1) = res
+        k1 = k1-2
+        error = err2+dabs(res-e2)+err3
+        if(error.gt.abserr) go to 40
+        abserr = error
+        result = res
+   40 continue
+c
+c           shift the table.
+c
+   50 if(n.eq.limexp) n = 2*(limexp/2)-1
+      ib = 1
+      if((num/2)*2.eq.num) ib = 2
+      ie = newelm+1
+      do 60 i=1,ie
+        ib2 = ib+2
+        epstab(ib) = epstab(ib2)
+        ib = ib2
+   60 continue
+      if(num.eq.n) go to 80
+      indx = num-n+1
+      do 70 i = 1,n
+        epstab(i)= epstab(indx)
+        indx = indx+1
+   70 continue
+   80 if(nres.ge.4) go to 90
+      res3la(nres) = result
+      abserr = oflow
+      go to 100
+c
+c           compute error estimate
+c
+   90 abserr = dabs(result-res3la(3))+dabs(result-res3la(2))
+     *  +dabs(result-res3la(1))
+      res3la(1) = res3la(2)
+      res3la(2) = res3la(3)
+      res3la(3) = result
+  100 abserr = dmax1(abserr,0.5d+01*epmach*dabs(result))
+      return
+      end
diff --git a/scilab/modules/differential_equations/src/fortran/dqk21.f b/scilab/modules/differential_equations/src/fortran/dqk21.f
new file mode 100644 (file)
index 0000000..9ea1458
--- /dev/null
@@ -0,0 +1,182 @@
+      subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
+c***begin prologue  dqk21
+c***date written   800101   (yymmdd)
+c***revision date  830518   (yymmdd)
+c***category no.  h2a1a2
+c***keywords  21-point gauss-kronrod rules
+c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
+c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
+c***purpose  to compute i = integral of f over (a,b), with error
+c                           estimate
+c                       j = integral of abs(f) over (a,b)
+c***description
+c
+c           integration rules
+c           standard fortran subroutine
+c           double precision version
+c
+c           parameters
+c            on entry
+c              f      - double precision
+c                       function subprogram defining the integrand
+c                       function f(x). the actual name for f needs to be
+c                       declared e x t e r n a l in the driver program.
+c
+c              a      - double precision
+c                       lower limit of integration
+c
+c              b      - double precision
+c                       upper limit of integration
+c
+c            on return
+c              result - double precision
+c                       approximation to the integral i
+c                       result is computed by applying the 21-point
+c                       kronrod rule (resk) obtained by optimal addition
+c                       of abscissae to the 10-point gauss rule (resg).
+c
+c              abserr - double precision
+c                       estimate of the modulus of the absolute error,
+c                       which should not exceed abs(i-result)
+c
+c              resabs - double precision
+c                       approximation to the integral j
+c
+c              resasc - double precision
+c                       approximation to the integral of abs(f-i/(b-a))
+c                       over (a,b)
+c
+c***references  (none)
+c***routines called  d1mach
+c***end prologue  dqk21
+c
+      double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
+     *  d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
+     *  resg,resk,reskh,result,uflow,wg,wgk,xgk
+      integer j,jtw,jtwm1
+      external f
+c
+      dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
+c
+c           the abscissae and weights are given for the interval (-1,1).
+c           because of symmetry only the positive abscissae and their
+c           corresponding weights are given.
+c
+c           xgk    - abscissae of the 21-point kronrod rule
+c                    xgk(2), xgk(4), ...  abscissae of the 10-point
+c                    gauss rule
+c                    xgk(1), xgk(3), ...  abscissae which are optimally
+c                    added to the 10-point gauss rule
+c
+c           wgk    - weights of the 21-point kronrod rule
+c
+c           wg     - weights of the 10-point gauss rule
+c
+c
+c gauss quadrature weights and kronron quadrature abscissae and weights
+c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
+c bell labs, nov. 1981.
+c
+      data wg  (  1) / 0.0666713443 0868813759 3568809893 332 d0 /
+      data wg  (  2) / 0.1494513491 5058059314 5776339657 697 d0 /
+      data wg  (  3) / 0.2190863625 1598204399 5534934228 163 d0 /
+      data wg  (  4) / 0.2692667193 0999635509 1226921569 469 d0 /
+      data wg  (  5) / 0.2955242247 1475287017 3892994651 338 d0 /
+c
+      data xgk (  1) / 0.9956571630 2580808073 5527280689 003 d0 /
+      data xgk (  2) / 0.9739065285 1717172007 7964012084 452 d0 /
+      data xgk (  3) / 0.9301574913 5570822600 1207180059 508 d0 /
+      data xgk (  4) / 0.8650633666 8898451073 2096688423 493 d0 /
+      data xgk (  5) / 0.7808177265 8641689706 3717578345 042 d0 /
+      data xgk (  6) / 0.6794095682 9902440623 4327365114 874 d0 /
+      data xgk (  7) / 0.5627571346 6860468333 9000099272 694 d0 /
+      data xgk (  8) / 0.4333953941 2924719079 9265943165 784 d0 /
+      data xgk (  9) / 0.2943928627 0146019813 1126603103 866 d0 /
+      data xgk ( 10) / 0.1488743389 8163121088 4826001129 720 d0 /
+      data xgk ( 11) / 0.0000000000 0000000000 0000000000 000 d0 /
+c
+      data wgk (  1) / 0.0116946388 6737187427 8064396062 192 d0 /
+      data wgk (  2) / 0.0325581623 0796472747 8818972459 390 d0 /
+      data wgk (  3) / 0.0547558965 7435199603 1381300244 580 d0 /
+      data wgk (  4) / 0.0750396748 1091995276 7043140916 190 d0 /
+      data wgk (  5) / 0.0931254545 8369760553 5065465083 366 d0 /
+      data wgk (  6) / 0.1093871588 0229764189 9210590325 805 d0 /
+      data wgk (  7) / 0.1234919762 6206585107 7958109831 074 d0 /
+      data wgk (  8) / 0.1347092173 1147332592 8054001771 707 d0 /
+      data wgk (  9) / 0.1427759385 7706008079 7094273138 717 d0 /
+      data wgk ( 10) / 0.1477391049 0133849137 4841515972 068 d0 /
+      data wgk ( 11) / 0.1494455540 0291690566 4936468389 821 d0 /
+c
+c
+c           list of major variables
+c           -----------------------
+c
+c           centr  - mid point of the interval
+c           hlgth  - half-length of the interval
+c           absc   - abscissa
+c           fval*  - function value
+c           resg   - result of the 10-point gauss formula
+c           resk   - result of the 21-point kronrod formula
+c           reskh  - approximation to the mean value of f over (a,b),
+c                    i.e. to i/(b-a)
+c
+c
+c           machine dependent constants
+c           ---------------------------
+c
+c           epmach is the largest relative spacing.
+c           uflow is the smallest positive magnitude.
+c
+c***first executable statement  dqk21
+      epmach = d1mach(4)
+      uflow = d1mach(1)
+c
+      centr = 0.5d+00*(a+b)
+      hlgth = 0.5d+00*(b-a)
+      dhlgth = dabs(hlgth)
+c
+c           compute the 21-point kronrod approximation to
+c           the integral, and estimate the absolute error.
+c
+      resg = 0.0d+00
+      fc = f(centr)
+      resk = wgk(11)*fc
+      resabs = dabs(resk)
+      do 10 j=1,5
+        jtw = 2*j
+        absc = hlgth*xgk(jtw)
+        fval1 = f(centr-absc)
+        fval2 = f(centr+absc)
+        fv1(jtw) = fval1
+        fv2(jtw) = fval2
+        fsum = fval1+fval2
+        resg = resg+wg(j)*fsum
+        resk = resk+wgk(jtw)*fsum
+        resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
+   10 continue
+      do 15 j = 1,5
+        jtwm1 = 2*j-1
+        absc = hlgth*xgk(jtwm1)
+        fval1 = f(centr-absc)
+        fval2 = f(centr+absc)
+        fv1(jtwm1) = fval1
+        fv2(jtwm1) = fval2
+        fsum = fval1+fval2
+        resk = resk+wgk(jtwm1)*fsum
+        resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
+   15 continue
+      reskh = resk*0.5d+00
+      resasc = wgk(11)*dabs(fc-reskh)
+      do 20 j=1,10
+        resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
+   20 continue
+      result = resk*hlgth
+      resabs = resabs*dhlgth
+      resasc = resasc*dhlgth
+      abserr = dabs((resk-resg)*hlgth)
+      if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
+     *  abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
+      if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
+     *  ((epmach*0.5d+02)*resabs,abserr)
+      return
+      end
diff --git a/scilab/modules/differential_equations/src/fortran/dqpsrt.f b/scilab/modules/differential_equations/src/fortran/dqpsrt.f
new file mode 100644 (file)
index 0000000..2f8bf8e
--- /dev/null
@@ -0,0 +1,129 @@
+      subroutine dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax)
+c***begin prologue  dqpsrt
+c***refer to  dqage,dqagie,dqagpe,dqawse
+c***routines called  (none)
+c***revision date  810101   (yymmdd)
+c***keywords  sequential sorting
+c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
+c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
+c***purpose  this routine maintains the descending ordering in the
+c            list of the local error estimated resulting from the
+c            interval subdivision process. at each call two error
+c            estimates are inserted using the sequential search
+c            method, top-down for the largest error estimate and
+c            bottom-up for the smallest error estimate.
+c***description
+c
+c           ordering routine
+c           standard fortran subroutine
+c           double precision version
+c
+c           parameters (meaning at output)
+c              limit  - integer
+c                       maximum number of error estimates the list
+c                       can contain
+c
+c              last   - integer
+c                       number of error estimates currently in the list
+c
+c              maxerr - integer
+c                       maxerr points to the nrmax-th largest error
+c                       estimate currently in the list
+c
+c              ermax  - double precision
+c                       nrmax-th largest error estimate
+c                       ermax = elist(maxerr)
+c
+c              elist  - double precision
+c                       vector of dimension last containing
+c                       the error estimates
+c
+c              iord   - integer
+c                       vector of dimension last, the first k elements
+c                       of which contain pointers to the error
+c                       estimates, such that
+c                       elist(iord(1)),...,  elist(iord(k))
+c                       form a decreasing sequence, with
+c                       k = last if last.le.(limit/2+2), and
+c                       k = limit+1-last otherwise
+c
+c              nrmax  - integer
+c                       maxerr = iord(nrmax)
+c
+c***end prologue  dqpsrt
+c
+      double precision elist,ermax,errmax,errmin
+      integer i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last,limit,maxerr,
+     *  nrmax
+      dimension elist(last),iord(last)
+c
+c           check whether the list contains more than
+c           two error estimates.
+c
+c***first executable statement  dqpsrt
+      if(last.gt.2) go to 10
+      iord(1) = 1
+      iord(2) = 2
+      go to 90
+c
+c           this part of the routine is only executed if, due to a
+c           difficult integrand, subdivision increased the error
+c           estimate. in the normal case the insert procedure should
+c           start after the nrmax-th largest error estimate.
+c
+   10 errmax = elist(maxerr)
+      if(nrmax.eq.1) go to 30
+      ido = nrmax-1
+      do 20 i = 1,ido
+        isucc = iord(nrmax-1)
+c ***jump out of do-loop
+        if(errmax.le.elist(isucc)) go to 30
+        iord(nrmax) = isucc
+        nrmax = nrmax-1
+   20    continue
+c
+c           compute the number of elements in the list to be maintained
+c           in descending order. this number depends on the number of
+c           subdivisions still allowed.
+c
+   30 jupbn = last
+      if(last.gt.(limit/2+2)) jupbn = limit+3-last
+      errmin = elist(last)
+c
+c           insert errmax by traversing the list top-down,
+c           starting comparison from the element elist(iord(nrmax+1)).
+c
+      jbnd = jupbn-1
+      ibeg = nrmax+1
+      if(ibeg.gt.jbnd) go to 50
+      do 40 i=ibeg,jbnd
+        isucc = iord(i)
+c ***jump out of do-loop
+        if(errmax.ge.elist(isucc)) go to 60
+        iord(i-1) = isucc
+   40 continue
+   50 iord(jbnd) = maxerr
+      iord(jupbn) = last
+      go to 90
+c
+c           insert errmin by traversing the list bottom-up.
+c
+   60 iord(i-1) = maxerr
+      k = jbnd
+      do 70 j=i,jbnd
+        isucc = iord(k)
+c ***jump out of do-loop
+        if(errmin.lt.elist(isucc)) go to 80
+        iord(k+1) = isucc
+        k = k-1
+   70 continue
+      iord(i) = last
+      go to 90
+   80 iord(k+1) = last
+c
+c           set maxerr and ermax.
+c
+   90 maxerr = iord(nrmax)
+      ermax = elist(maxerr)
+      return
+      end
diff --git a/scilab/modules/differential_equations/src/fortran/epsalg.f b/scilab/modules/differential_equations/src/fortran/epsalg.f
deleted file mode 100644 (file)
index 47c4171..0000000
+++ /dev/null
@@ -1,178 +0,0 @@
-C/MEMBR ADD NAME=EPSALG,SSI=0
-      subroutine epsalg(n, epstab, result, abserr, res3la, nres)
-c
-c     based on quadpack routine epsalg
-c     ******************************************************
-c
-c           purpose
-c              the routine transforms a given sequence of
-c              approximations, by means of the epsilon
-c              algorithm of p. wynn.
-c
-c              an estimate of the absolute error is also given.
-c              the condensed epsilon table is computed. only those
-c              elements needed for the computation of the
-c              next diagonal are preserved.
-c
-c           calling sequence
-c              call epsalg (n,epstab,result,abserr,res3la,nres)
-c
-c           parameters
-c              n      - epstab(n) contains the new element in the
-c                       first column of the epsilon table.
-c
-c              epstab - one dimensional array containing the
-c                       elements of the two lower diagonals of
-c                       the triangular epsilon table.
-c                       the elements are numbered starting at the
-c                       right-hand corner of the triangle.
-c                       the dimension should be at least n+2.
-c
-c              result - resulting approximation to the integral
-c
-c              abserr - estimate of the absolute error computed from
-c                       result and the 3 previous /results/
-c
-c              res3la - array containing the last 3 /results/
-c
-c              nres   - number of calls to the routine
-c                       (should be zero at first call)
-c
-c     ******************************************************
-c     .. scalar arguments ..
-      double precision abserr, result
-      integer n, nres
-c     .. array arguments ..
-      double precision epstab(52), res3la(3)
-c     ..
-c     .. local scalars ..
-      double precision delta1, delta2, delta3, e0, e1, e1abs, e2, e3,
-     * epmach,epsinf, err1, err2, err3, error, oflow, res, ss, tol1,
-     * tol2,tol3
-      integer i, ib2, ib, ie, ind, k1, k2, k3, limexp, newelm, num
-c     .. function references ..
-      double precision dlamch
-c     ..
-c
-c            machine dependent constants
-c             -------------------------
-c            /limexp/ is the maximum number of elements the epsilon
-c            table can contain. if this number is reached, the upper
-c            diagonal of the epsilon table is deleted.
-c
-      data limexp /50/
-      epmach=dlamch('p')
-      oflow=dlamch('o')
-c
-c           list of major variables
-c           -----------------------
-c           e0     - the 4 elements on which the
-c           e1       computation of a new element in
-c           e2       the epsilon table is based
-c           e3                 e0
-c                        e3    e1    new
-c                              e2
-c           newelm - number of elements to be computed in the new
-c                    diagonal
-c           error  - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
-c           result - the element in the new diagonal with least
-c                    error
-c
-      nres = nres + 1
-      abserr = oflow
-      result = epstab(n)
-      if (n.lt.3) go to 200
-      epstab(n+2) = epstab(n)
-      newelm = (n-1)/2
-      epstab(n) = oflow
-      num = n
-      k1 = n
-      do 80 i=1,newelm
-         k2 = k1 - 1
-         k3 = k1 - 2
-         res = epstab(k1+2)
-         e0 = epstab(k3)
-         e1 = epstab(k2)
-         e2 = res
-         e1abs = abs(e1)
-         delta2 = e2 - e1
-         err2 = abs(delta2)
-         tol2 = max(abs(e2),e1abs)*epmach
-         delta3 = e1 - e0
-         err3 = abs(delta3)
-         tol3 = max(e1abs,abs(e0))*epmach
-         if (err2.gt.tol2 .or. err3.gt.tol3) go to 20
-c
-c           if e0, e1 and e2 are equal to within machine
-c           accuracy, convergence is assumed
-c           result = e2
-c           abserr = abs(e1-e0)+abs(e2-e1)
-c
-         result = res
-         abserr = err2 + err3
-         go to 200
-   20    e3 = epstab(k1)
-         epstab(k1) = e1
-         delta1 = e1 - e3
-         err1 = abs(delta1)
-         tol1 = max(e1abs,abs(e3))*epmach
-c
-c           if two elements are very close to each other, omit
-c           a part of the table by adjusting the value of n
-c
-         if (err1.lt.tol1 .or. err2.lt.tol2 .or. err3.lt.tol3) goto 40
-         ss = 0.10d+01/delta1 + 0.10d+01/delta2 - 0.10d+01/delta3
-         epsinf = abs(ss*e1)
-c
-c           test to detect irregular behaviour in the table, and
-c           eventually omit a part of the table adjusting the value
-c           of n
-c
-         if (epsinf.gt.0.10d-03) go to 60
-   40    n = i + i - 1
-         go to 100
-c
-c           compute a new element and eventually adjust
-c           the value of result
-c
-   60    res = e1 + 0.10d+01/ss
-         epstab(k1) = res
-         k1 = k1 - 2
-         error = err2 + abs(res-e2) + err3
-         if (error.gt.abserr) go to 80
-         abserr = error
-         result = res
-   80 continue
-c
-c           shift the table
-c
-  100 if (n.eq.limexp) n = 2*(limexp/2) - 1
-      ib = 1
-      if ((num/2)*2.eq.num) ib = 2
-      ie = newelm + 1
-      do 120 i=1,ie
-         ib2 = ib + 2
-         epstab(ib) = epstab(ib2)
-         ib = ib2
-  120 continue
-      if (num.eq.n) go to 160
-      ind = num - n + 1
-      do 140 i=1,n
-         epstab(i) = epstab(ind)
-         ind = ind + 1
-  140 continue
-  160 if (nres.ge.4) go to 180
-      res3la(nres) = result
-      abserr = oflow
-      go to 200
-c
-c           compute error estimate
-c
-  180 abserr = abs(result-res3la(3)) + abs(result-res3la(2)) +
-     *abs(result-res3la(1))
-      res3la(1) = res3la(2)
-      res3la(2) = res3la(3)
-      res3la(3) = result
-  200 abserr = max(abserr,5.0d+00*epmach*abs(result))
-      return
-      end
diff --git a/scilab/modules/differential_equations/src/fortran/order.f b/scilab/modules/differential_equations/src/fortran/order.f
deleted file mode 100644 (file)
index 1ef9776..0000000
+++ /dev/null
@@ -1,134 +0,0 @@
-C/MEMBR ADD NAME=ORDER,SSI=0
-      subroutine order(limit, last, maxerr, ermax, elist, iord,liord,
-     * nrmax)
-c
-c     based on quadpack routine order
-c     ******************************************************
-c
-c           purpose
-c              this routine maintains the descending ordering
-c              in the list of the local error estimates
-c              resulting from the interval subdivision
-c              process. at each call two error estimates
-c              are inserted using the sequential search
-c              method . top-down for the largest error
-c              estimate,  bottom-up for the smallest error
-c              estimate.
-c
-c           calling sequence
-c              call order
-c              (limit,last,maxerr,ermax,elist,iord,liord,nrmax)
-c
-c             parameters (meaning at output)
-c              limit  - maximum number of error estimates the list
-c                       can contain
-c
-c              last   - number of error estimates currently
-c                       in the list. elist(last) contains
-c                       the smallest error estimate.
-c
-c              maxerr - maxerr points to the nrmax-th largest error
-c                       estimate currently in the list.
-c
-c              ermax  - nrmax-th largest error estimate
-c                       ermax = elist(maxerr)
-c
-c              elist  - array of dimension last containing
-c                       the error estimates
-c
-c              iord   - array containing pointers to elist so
-c                       that iord(1) points to the largest
-c                       error estimate,...,iord(last) to the
-c                       smallest error estimate
-c
-c              liord  - dimension of iord
-c
-c              nrmax  - maxerr = iord(nrmax)
-c
-c     ******************************************************
-c
-c     .. scalar arguments ..
-      double precision ermax
-      integer last, limit, liord, maxerr, nrmax
-c     .. array arguments ..
-      double precision elist(last)
-      integer iord(liord)
-c     ..
-c     .. scalars in common ..
-      integer jupbnd
-c     ..
-c     .. local scalars ..
-      double precision errmax, errmin
-      integer i, ibeg, ido, isucc, j, jbnd, k
-c     ..
-      common /dqa001/ jupbnd
-c
-c            check whether the list contains more than
-c            two error estimates
-c
-      if (last.gt.2) go to 20
-      iord(1) = 1
-      iord(2) = 2
-      go to 180
-c
-c           this part of the routine is only executed
-c           if, due to a difficult integrand, subdivision
-c           increased the error estimate. in the normal case
-c           the insert procedure should start after the
-c           nrmax-th largest error estimate.
-c
-   20 errmax = elist(maxerr)
-      if (nrmax.eq.1) go to 60
-      ido = nrmax - 1
-      do 40 i=1,ido
-         isucc = iord(nrmax-1)
-         if (errmax.le.elist(isucc)) go to 60
-         iord(nrmax) = isucc
-         nrmax = nrmax - 1
-   40 continue
-c
-c           compute the number of elements in the list to
-c           be maintained in descending order. this number
-c           depends on the number of subdivisions still
-c           allowed
-c
-   60 jupbnd = last
-      if (last.gt.(limit/2+2)) jupbnd = limit + 3 - last
-      errmin = elist(last)
-c
-c           insert errmax by traversing the list top-down
-c           starting comparison from the element
-c           elist(iord(nrmax+1))
-c
-      jbnd = jupbnd - 1
-      ibeg = nrmax + 1
-      if (ibeg.gt.jbnd) go to 100
-      do 80 i=ibeg,jbnd
-         isucc = iord(i)
-         if (errmax.ge.elist(isucc)) go to 120
-         iord(i-1) = isucc
-   80 continue
-  100 iord(jbnd) = maxerr
-      iord(jupbnd) = last
-      go to 180
-c
-c           insert errmin by traversing the list bottom-up
-c
-  120 iord(i-1) = maxerr
-      k = jbnd
-      do 140 j=i,jbnd
-         isucc = iord(k)
-         if (errmin.lt.elist(isucc)) go to 160
-         iord(k+1) = isucc
-         k = k - 1
-  140 continue
-      iord(i) = last
-      go to 180
-  160 iord(k+1) = last
-c
-c           set maxerr and ermax
-c
-  180 maxerr = iord(nrmax)
-      ermax = elist(maxerr)
-      return
-      end
diff --git a/scilab/modules/differential_equations/src/fortran/quarul.f b/scilab/modules/differential_equations/src/fortran/quarul.f
deleted file mode 100644 (file)
index b35ec0f..0000000
+++ /dev/null
@@ -1,159 +0,0 @@
-C/MEMBR ADD NAME=QUARUL,SSI=0
-      subroutine quarul(f, a, b, result, abserr, resabs, resasc)
-c
-c     based on quadpack routine quarul
-c     ******************************************************
-c
-c           purpose
-c              to compute i = integral of f over (a,b), with error
-c                             estimate
-c                         j = integral of abs(f) over (a,b)
-c
-c           calling sequence
-c              call quarul (f,a,b,result,abserr,resabs,resasc)
-c
-c           parameters
-c              f      - function subprogram defining the integrand
-c                       function f(x). the actual name for f needs
-c                       to be declared e x t e r n a l in the
-c                       calling program
-c
-c              a      - lower limit of integration
-c
-c              b      - upper limit of integration
-c
-c              result - approximation to the integral i.
-c                       result is calculated by applying
-c                       the 21-point gauss-kronrod rule
-c                       (resk), obtained by optimal
-c                       addition of abscissae to the
-c                       10-point gauss  rule (resg).
-c
-c              abserr - estimate of the modulus of the
-c                       absolute error, which should not
-c                       exceed abs(i-result)
-c              resabs - approximation to the integral j
-c
-c              resasc - approximation to the integral of
-c                       abs(f-i/(b-a)) over (a,b)
-c
-c     ******************************************************
-c     .. scalar arguments ..
-      double precision a, abserr, b, resabs, resasc, result
-c     .. function arguments ..
-      double precision f
-      external f
-c recuperation d'erreur
-      integer       iero
-      common/ierajf/iero
-c     ..
-c     .. local scalars ..
-      double precision absc, centre, dhlgth, epmach, fc, fsum, fval1,
-     * fval2,hlgth, resg, resk, reskh, uflow
-      integer j
-c     .. local arrays ..
-      double precision fv1(10), fv2(10), wg(10), wgk(11), xgk(11)
-c     .. function references ..
-      double precision dlamch
-c     ..
-c
-c            the abscissae and weights are given for the
-c            interval (-1,1) . because of symmetry only the
-c            positive abscissae and their corresponding
-c            weights are given.
-c            xgk    - abscissae of the 21-point gauss-kronrod rule
-c                     xgk(2), xgk(4), .... abscissae of the 10-point
-c                     gauss rule
-c                     xgk(1), xgk(3), .... abscissae which
-c                     are optimally added to the 10-point
-c                     gauss rule
-c            wgk    - weights of the 21-point gauss-kronrod rule
-c            wg     - weights of the 10-point gauss rule,
-c                     corresponding to the abscissae xgk(2),
-c                     xgk(4), ... wg(1), wg(3), ... are set
-c                     to zero.
-c
-      data xgk(1), xgk(2), xgk(3), xgk(4), xgk(5), xgk(6), xgk(7),xgk(8)
-     *, xgk(9), xgk(10), xgk(11) /0.99565716302580808073552728070d+00,
-     *0.97390652851717172007796401210d+00,
-     *0.93015749135570822600120718010d+00,
-     *0.86506336668898451073209668840d+00,
-     *0.78081772658641689706371757830d+00,
-     *0.67940956829902440623432736510d+00,
-     *0.56275713466860468333900009930d+00,
-     *0.43339539412924719079926594320d+00,
-     *0.29439286270146019813112660310d+00,
-     *0.14887433898163121088482600110d+00,0.0d+0/
-      data wgk(1), wgk(2), wgk(3), wgk(4), wgk(5), wgk(6), wgk(7),wgk(8)
-     *, wgk(9), wgk(10), wgk(11) /0.11694638867371874278064396060d-01,
-     *0.32558162307964727478818972460d-01,
-     *0.54755896574351996031381300240d-01,
-     *0.75039674810919952767043140920d-01,
-     *0.93125454583697605535065465080d-01,
-     *0.10938715880229764189921059030d+00,
-     *0.12349197626206585107795810980d+00,
-     *0.13470921731147332592805400180d+00,
-     *0.14277593857706008079709427310d+00,
-     *0.14773910490133849137484151600d+00,
-     *0.14944555400291690566493646840d+00/
-      data wg(1), wg(2), wg(3), wg(4), wg(5), wg(6), wg(7), wg(8),wg(9),
-     * wg(10) /0.0d+0,0.66671344308688137593568809890d-01,0.0d+0,
-     *0.14945134915058059314577633970d+00,0.0d+0,
-     *0.21908636251598204399553493420d+00,0.0d+0,
-     *0.26926671930999635509122692160d+00,0.0d+0,
-     *0.29552422471475287017389299470d+00/
-      epmach=dlamch('p')
-      uflow=dlamch('u')
-c
-c           list of major variables
-c           ----------------------
-c           centre - mid point of the interval
-c           hlgth  - half length of the interval
-c           absc   - abscissa
-c           fval*  - function value
-c           resg   - 10-point gauss formula
-c           resk   - 21-point gauss-kronrod formula
-c           reskh  - approximation to mean value of f over
-c                    (a,b), i.e. to i/(b-a)
-c
-      centre = 0.50d+00*(a+b)
-      hlgth = 0.50d+00*(b-a)
-      dhlgth = abs(hlgth)
-c
-c           compute the 21-point gauss-kronrod approximation to
-c           the integral, and estimate the absolute error
-c
-      resg = 0.0d+00
-      fc = f(centre)
-      if(iero.ne.0) return
-      resk = wgk(11)*fc
-      resabs = abs(resk)
-      do 20 j=1,10
-         absc = hlgth*xgk(j)
-         fval1 = f(centre-absc)
-         if(iero.ne.0) return
-         fval2 = f(centre+absc)
-         if(iero.ne.0) return
-         fv1(j) = fval1
-         fv2(j) = fval2
-         fsum = fval1 + fval2
-         resg = resg + wg(j)*fsum
-         resk = resk + wgk(j)*fsum
-         resabs = resabs + wgk(j)*(abs(fval1)+abs(fval2))
-   20 continue
-      reskh = resk*0.50d+00
-      resasc = wgk(11)*abs(fc-reskh)
-      do 40 j=1,10
-         resasc = resasc + wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)
-     *   )
-   40 continue
-      result = resk*hlgth
-      resabs = resabs*dhlgth
-      resasc = resasc*dhlgth
-      abserr = abs((resk-resg)*hlgth)
-      if (resasc.ne.0.0d+00 .and. abserr.ne.0.0d+00) abserr =resasc*
-     *min(0.10d+01,(0.20d+03*abserr/resasc)**1.50d+0)
-      if (resabs.gt.uflow/(0.50d+02*epmach)) abserr =max(epmach*resabs*
-     *0.50d+02,abserr)
-      return
-      end
index 37165a2..d29b6fa 100644 (file)
 x0 = 0;
 x1=0:0.1:2*%pi;
 computed = integrate('sin(x)','x',x0,x1);
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
 expected = 1 - cos(x1);
 if (norm(computed-expected) > 10* %eps) then bugmes();quit;end
index 0a6772e..99eb60c 100644 (file)
@@ -1,5 +1,91 @@
-// to check that intg works
-function y=f(x),y=x*sin(30*x)/sqrt(1-((x/(2*%pi))^2)),endfunction
-exact=-2.5432596188;
-I=intg(0,2*%pi,f);
-if abs(exact-I) > 1 then bugmes();quit;end
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2007-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- ENGLISH IMPOSED -->
+// Run with test_run('differential_equations','intg',['no_check_error_output'])
+// Function written in the Scilab language
+function y = f(x), y = x*sin(30*x)/sqrt(1-((x/(2*%pi))^2)), endfunction
+exact = -2.5432596188;
+I = intg(0, 2*%pi, f);
+if abs(exact-I) > 1e-9 then bugmes();quit;end
+// Function with an argument written in the Scilab language
+function y = f1(x, w), y = x*sin(w*x)/sqrt(1-((x/(2*%pi))^2)), endfunction
+I = intg(0, 2*%pi, list(f1, 30));
+if abs(exact-I) > 1e-9 then bugmes();quit;end
+// Function written in Fortran (a Fortran compiler is required)
+// define a Fortran function
+cd TMPDIR;
+F=['      double precision function ffun(x)'
+   '      double precision x, pi'
+   '      pi = 3.14159265358979312d+0'
+   '      ffun = x*sin(30.0d+0*x)/sqrt(1.0d+0-(x/(2.0d+0*pi))**2)'
+   '      return'
+   '      end'];
+mputl(F, fullfile(TMPDIR, 'ffun.f'));
+// compile the function
+l = ilib_for_link('ffun', fullfile(TMPDIR, 'ffun.f'), [], 'f');
+   Generate a loader file
+   Generate a Makefile
+   ilib_gen_Make: Copy compilation files (Makefile*, libtool...) to TMPDIR
+   ilib_gen_Make: Copy TMPDIR/ffun.f to TMPDIR
+   ilib_gen_Make: Modification of the Makefile in TMPDIR.
+   Running the Makefile
+   Generate a cleaner file
+// add the function to the working environment
+link(l, 'ffun', 'f');
+Shared archive loaded.
+Link done.
+// integrate the function
+I = intg(0, 2*%pi, 'ffun');
+abs(exact-I);
+if abs(exact-I) > 1e-9 then bugmes();quit;end
+// Function written in C (a C compiler is required)
+// define a C function
+C=['#include <math.h>'
+   'double cfun(double *x)'
+   '{'
+   '  double y, pi = 3.14159265358979312;'
+   '  y = *x/(2.0e0*pi);'
+   '  return *x*sin(30.0e0**x)/sqrt(1.0e0-y*y);'
+   '}'];
+mputl(C, fullfile(TMPDIR, 'cfun.c'));
+// compile the function
+l = ilib_for_link('cfun', fullfile(TMPDIR, 'cfun.c'), [], 'c');
+   Generate a loader file
+   Generate a Makefile
+   ilib_gen_Make: Copy compilation files (Makefile*, libtool...) to TMPDIR
+   ilib_gen_Make: Copy TMPDIR/cfun.c to TMPDIR
+   ilib_gen_Make: Modification of the Makefile in TMPDIR.
+   Running the Makefile
+   Generate a cleaner file
+// add the function to the working environment
+link(l, 'cfun', 'c');
+Shared archive loaded.
+Link done.
+// integrate the function
+I = intg(0, 2*%pi, 'cfun');
+if abs(exact-I) > 1e-9 then bugmes();quit;end
+// Test third output argument
+[i, err, ierr] = intg(0, 1, f);
+if abs(ierr) <> 0 then bugmes();quit;end
+function y = f(x), y = cos(x); endfunction
+Warning : redefining function: f                       . Use funcprot(0) to avoid this message
+
+[i, err, ierr] = intg(0, %pi, f);
+Warning: argument function is detected to be very smooth.
+Reduce integration interval if you think edges have been missed.
+Warning: Round-off error detected, the requested tolerance (or default) cannot be achieved. Try using bigger tolerances.
+if abs(ierr) <> 2 then bugmes();quit;end
+// IEEE compatibility
+// Error 264: "Wrong value for argument #i: Must not contain NaN or Inf."
+if execstr("I = intg(%inf, 0, f)", 'errcatch')    <> 264 then bugmes();quit;end
+if execstr("I = intg(-%inf, 0, f)", 'errcatch')   <> 264 then bugmes();quit;end
+if execstr("I = intg(%nan, 0, f)", 'errcatch')    <> 264 then bugmes();quit;end
+if execstr("I = intg(0, %inf, f)", 'errcatch')    <> 264 then bugmes();quit;end
+if execstr("I = intg(0, -%inf, f)", 'errcatch')   <> 264 then bugmes();quit;end
+if execstr("I = intg(0, %nan, f)", 'errcatch')    <> 264 then bugmes();quit;end
+if execstr("I = intg(%nan, %nan, f)", 'errcatch') <> 264 then bugmes();quit;end
index ba9f957..aca1280 100644 (file)
@@ -6,9 +6,78 @@
 // =============================================================================
 
 // <-- CLI SHELL MODE -->
+// <-- ENGLISH IMPOSED -->
 
-// to check that intg works
-function y=f(x),y=x*sin(30*x)/sqrt(1-((x/(2*%pi))^2)),endfunction
-exact=-2.5432596188;
-I=intg(0,2*%pi,f);
-if abs(exact-I) > 1 then pause,end
\ No newline at end of file
+// Run with test_run('differential_equations','intg',['no_check_error_output'])
+
+// Function written in the Scilab language
+function y = f(x), y = x*sin(30*x)/sqrt(1-((x/(2*%pi))^2)), endfunction
+exact = -2.5432596188;
+I = intg(0, 2*%pi, f);
+if abs(exact-I) > 1e-9 then pause, end
+
+// Function with an argument written in the Scilab language
+function y = f1(x, w), y = x*sin(w*x)/sqrt(1-((x/(2*%pi))^2)), endfunction
+I = intg(0, 2*%pi, list(f1, 30));
+if abs(exact-I) > 1e-9 then pause, end
+
+// Function written in Fortran (a Fortran compiler is required)
+// define a Fortran function
+cd TMPDIR;
+F=['      double precision function ffun(x)'
+   '      double precision x, pi'
+   '      pi = 3.14159265358979312d+0'
+   '      ffun = x*sin(30.0d+0*x)/sqrt(1.0d+0-(x/(2.0d+0*pi))**2)'
+   '      return'
+   '      end'];
+mputl(F, fullfile(TMPDIR, 'ffun.f'));
+
+// compile the function
+l = ilib_for_link('ffun', fullfile(TMPDIR, 'ffun.f'), [], 'f');
+
+// add the function to the working environment
+link(l, 'ffun', 'f');
+
+// integrate the function
+I = intg(0, 2*%pi, 'ffun');
+abs(exact-I);
+if abs(exact-I) > 1e-9 then pause,end
+
+// Function written in C (a C compiler is required)
+// define a C function
+C=['#include <math.h>'
+   'double cfun(double *x)'
+   '{'
+   '  double y, pi = 3.14159265358979312;'
+   '  y = *x/(2.0e0*pi);'
+   '  return *x*sin(30.0e0**x)/sqrt(1.0e0-y*y);'
+   '}'];
+mputl(C, fullfile(TMPDIR, 'cfun.c'));
+
+// compile the function
+l = ilib_for_link('cfun', fullfile(TMPDIR, 'cfun.c'), [], 'c');
+
+// add the function to the working environment
+link(l, 'cfun', 'c');
+
+// integrate the function
+I = intg(0, 2*%pi, 'cfun');
+if abs(exact-I) > 1e-9 then pause,end
+
+// Test third output argument
+[i, err, ierr] = intg(0, 1, f);
+if abs(ierr) <> 0 then pause,end
+
+function y = f(x), y = cos(x); endfunction
+[i, err, ierr] = intg(0, %pi, f);
+if abs(ierr) <> 2 then pause,end
+
+// IEEE compatibility
+// Error 264: "Wrong value for argument #i: Must not contain NaN or Inf."
+if execstr("I = intg(%inf, 0, f)", 'errcatch')    <> 264 then pause,end
+if execstr("I = intg(-%inf, 0, f)", 'errcatch')   <> 264 then pause,end
+if execstr("I = intg(%nan, 0, f)", 'errcatch')    <> 264 then pause,end
+if execstr("I = intg(0, %inf, f)", 'errcatch')    <> 264 then pause,end
+if execstr("I = intg(0, -%inf, f)", 'errcatch')   <> 264 then pause,end
+if execstr("I = intg(0, %nan, f)", 'errcatch')    <> 264 then pause,end
+if execstr("I = intg(%nan, %nan, f)", 'errcatch') <> 264 then pause,end