Bug #10574 fixed - Runge-Kutta-Method failed for vector [x, 1] with x > 9.
[scilab.git] / scilab / modules / differential_equations / src / fortran / odeint.f
1 c     ====================================
2 c .   Runge-Kutta driver with adaptive stepsize control.  Integrate
3 c .   the starting values ystart(1:nvar) from x1 to x2 with accuracy 
4 c .   eps, storing intermediate results in the common block /path/.
5 c .   h1 should be set as a guessed first stepsize, hmin as the 
6 c .   minimum allowed stepsize (can be zero).
7 c .   On output, nok and nbad are the number of good and bad (but
8 c .   retried and fixed) steps taken, and ystart is replaced by 
9 c .   values at the end of the integration interval.
10 c .   derivs is the user-supplied subroutine for calculating the
11 c .   right-hand side derivatives, while rkqs is the name of 
12 c .   the stepper routine to be used.  
13 c .   /path/ contains its own information about how often an
14 c .   intermediate value is to be stored.
15 c .
16 c .   The original version has been modified to replace statically
17 c .   allocated arrays y, yscal and dxdy by parts of rwork argument
18 c     array + blas use. Serge Steer INRIA- feb 2012
19
20 c     ====================================
21       subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,
22      $     rkqc,rwork)
23       external derivs,rkqc
24       integer maxstp,kount,nvar,i,nok,nbad,nstp
25       double precision two,zero,tiny,x,h
26       parameter (maxstp=10000,two=2.0d0,zero=0.0d0,tiny=1.d-30)
27       double precision x1,x2,eps,h1,hmin,hdid,hnext
28       double precision ystart(nvar)
29       double precision rwork(*)
30       character*80 messag
31       integer ly,lyscal,ldydx,lwork
32       integer iero
33       common/ierode/iero
34 c     
35       ly=1
36       lyscal=ly+nvar
37       ldydx=lyscal+nvar
38       lwork=ldydx+nvar
39
40       iero=0
41       if ( abs(x2-x1).le.tiny) return
42       x=x1
43       h=sign(h1,x2-x1)
44       nok=0
45       nbad=0
46       kount=0
47       call dcopy(nvar,ystart,1,rwork(ly),1)
48
49       do 16 nstp=1,maxstp
50          call derivs(nvar,x,rwork(ly),rwork(ldydx))
51          if (iero.gt.0) return 
52          do 12 i=0,nvar-1
53             rwork(lyscal+i)=abs(rwork(ly+i))+abs(h*rwork(ldydx+i))+tiny
54  12      continue
55
56          if((x+h-x2)*(x+h-x1).gt.zero) h=x2-x
57          call rkqc(rwork(ly),rwork(ldydx),nvar,x,h,eps,rwork(lyscal),
58      $        hdid,hnext,derivs,rwork(lwork))
59          if(iero.gt.0) return
60          if(hdid.eq.h)then
61             nok=nok+1
62          else
63             nbad=nbad+1
64          endif
65          if((x-x2)*(x2-x1).ge.zero)then
66             call dcopy(nvar,rwork(ly),1,ystart,1)
67             return
68          endif
69          if(abs(hnext).lt.hmin) then
70             write(messag, 17) hnext,hmin
71             hnext=hmin
72          endif
73          h=hnext
74  16   continue
75       iero=-1
76 c     print *, 'Trop d''iterations a faire pour la precision demandee.'
77       return
78  17   format('stepsize ',e10.3,' smaller than minimum ',e10.3)
79       end
80 c     ====================================
81       
82