* Bugs #5728 10229 invalid - Differential_equations: Quadpack update
[scilab.git] / scilab / modules / differential_equations / sci_gateway / fortran / intg.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) INRIA
3 c ...
4 c
5 c This file must be used under the terms of the CeCILL.
6 c This source file is licensed as described in the file COPYING, which
7 c you should have received as part of this distribution.  The terms
8 c are also available at
9 c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 c
11       subroutine intg
12 c     --------------------------------------------
13 c     Scilab intg
14 c      implicit undefined (a-z)
15       character*(4) fname
16       character*6   namef
17       include 'stack.h'
18       integer iero
19       common/ierajf/iero
20       common/cintg/namef
21       external bintg,fintg
22       double precision epsa,epsr,a,b,val,abserr
23       logical getexternal, getscalar,type ,cremat
24       integer topk,lr,katop,kydot,top2,lra,lrb,lc,lier
25       integer iipal,lpal,lw,liw,lpali,ier
26       integer iadr,sadr,vfinite
27       external setfintg
28 c
29       iadr(l)=l+l-1
30       sadr(l)=(l/2)+1
31       fname='intg'
32       if(rhs.lt.3) then
33          call error(39)
34          return
35       endif
36       if(lhs.gt.3) then
37          call error(59)
38          return
39       endif
40       type=.false.
41       top2=top
42       topk=top
43       if(rhs.eq.5) then
44          if (.not.getscalar(fname,topk,top,lr)) return
45          epsr=stk(lr)
46          top=top-1
47       else
48          epsr=1.0d-8
49       endif
50       if (rhs.ge.4) then
51          if (.not.getscalar(fname,topk,top,lr)) return
52          epsa=stk(lr)
53          top=top-1
54       else
55          epsa=1.0d-14
56       endif
57 c     cas standard
58       if (.not.getexternal(fname,topk,top,namef,type,
59      $     setfintg)) return
60       kydot=top
61       top=top-1
62       if (.not.getscalar(fname,topk,top,lrb)) return
63       b=stk(lrb)
64       if (isanan(b).eq.1.or.vfinite(1,b).eq.0) then
65          err = 2
66          call error(264)
67          return
68       endif
69       top=top-1
70       katop=top
71       if (.not.getscalar(fname,topk,top,lra)) return
72       a=stk(lra)
73       if (isanan(a).eq.1.or.vfinite(1,a).eq.0) then
74          err = 1
75          call error(264)
76          return
77       endif
78 c     tableaux de travail
79       top=top2+1
80       lw=3000
81       if (.not.cremat(fname,top,0,1,lw,lpal,lc)) return
82       top=top+1
83 c     tableau de travail entier necessaire
84       liw=3000/8+2
85       if (.not.cremat(fname,top,0,1,iadr(liw)+1,lpali,lc)) return
86       top=top+1
87 c
88 c     external scilab
89 c
90       iipal=iadr(lstk(top))
91       istk(iipal)=1
92       istk(iipal+1)=iipal+2
93       istk(iipal+2)=kydot
94       istk(iipal+3)=katop
95       lstk(top+1)=sadr(iipal+4)
96       if(type) then
97          call dqags(fintg,a,b,epsa,epsr,val,abserr,
98      +        neval,ier,lw/4,lw,last,stk(lpali),stk(lpal))
99       else
100          call dqags(bintg,a,b,epsa,epsr,val,abserr,
101      +        neval,ier,lw/4,lw,last,stk(lpali),stk(lpal))
102       endif
103       if(err.gt.0.or.err1.gt.0) return
104       if(ier.gt.0) then
105          if (lhs.le.2) then
106             if (ier.eq.1) then
107                call erro('Error: Maximum number of subdivisons '//
108      &            'achieved. Splitting the interval might help.')
109                return
110             elseif (ier.eq.2) then
111                call erro('Error: Round-off error detected, the '//
112      &         'requested tolerance (or default) cannot be achieved.'//
113      &         ' Try using bigger tolerances.')
114                return
115             elseif (ier.eq.3) then
116                call erro('Error: Bad integrand behavior occurs at '//
117      &            'some points of the integration interval.')
118                return
119             elseif (ier.eq.4) then
120                call erro('Error: Convergence problem, round-off '//
121      &            'error detected. Try using bigger tolerances.')
122                return
123             elseif (ier.eq.5) then
124                call erro('Error: The integral is probably '//
125      &            'divergent, or slowly convergent.')
126                return
127             else
128 c           ier = 6
129                call erro('Error: Invalid input, absolute '//
130      &            'tolerance <= 0 and relative tolerance < 2.e-14.')
131                return
132             endif
133          else
134             if (ier.eq.1) then
135                call msgstxt('Warning: Maximum number of '//
136      &            'subdivisions achieved. '//
137      &            'Splitting the interval might help.')
138             elseif (ier.eq.2) then
139                call msgstxt('Warning: Round-off error detected, '//
140      &         'the requested tolerance (or default) cannot be '//
141      &         'achieved. Try using bigger tolerances.')
142             elseif (ier.eq.3) then
143                call msgstxt('Warning: Bad integrand behavior '//
144      &            'occurs at some points of the integration interval.')
145             elseif (ier.eq.4) then
146                call msgstxt('Warning: Convergence problem, round-off'//
147      &            'error detected. Try using bigger tolerances.')
148             elseif (ier.eq.5) then
149                call msgstxt('Warning: The integral is probably '//
150      &            'divergent, or slowly convergent.')
151              else
152 c           ier = 6
153                call msgstxt('Warning: Invalid inputs, absolute '//
154      &            'tolerance <= 0 and relative tolerance < 2.e-14.')
155             endif
156         endif
157       endif
158       top=top2-rhs+1
159       stk(lra)=val
160       if(lhs.ge.2) then
161          top=top+1
162          stk(lrb)=abserr
163          if (lhs.eq.3) then
164             top=top+1
165             ilier=iadr(lstk(top))
166             istk(ilier)=1
167             istk(ilier+1)=1
168             istk(ilier+2)=1
169             istk(ilier+3)=0
170             lier=sadr(ilier+4)
171             stk(lier)=ier
172             lstk(top+1)=lier+1
173             return
174          endif
175          return
176       endif
177       return
178       end
179