fix import/export of ierode common on Windows
[scilab.git] / scilab / modules / differential_equations / sci_gateway / fortran / sci_f_dasrt.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) INRIA/ENPC
3 c ...
4
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.1-en.txt
10 c
11       subroutine dasrti(fname)
12 c ====================================================================
13 C     dasrt 
14 c ====================================================================
15 c
16       INCLUDE 'stack.h'
17 c
18       character*(*) fname
19       character*(nlgh+1) namjac
20       common/cjac/namjac
21       integer iadr,sadr,gettype
22 c
23       double precision atol,rtol,t0
24       integer info(15),topk,topw
25       logical hotstart,type,getexternal,getrvect
26       logical checkrhs,checklhs,getrmat,cremat,getscalar
27       double precision tout,tstop,maxstep,stepin
28       character*(nlgh+1) namer,namej,names
29       common /dassln/ namer,namej,names
30       external bresd,bjacd,bsurfd
31       external setfresd,setfjacd,setfsurfd
32 c     
33       data atol/1.d-7/,rtol/1.d-9/
34 c     
35       iadr(l)=l+l-1
36       sadr(l)=(l/2)+1
37
38 c     SCILAB function : dasrt
39 c     --------------------------
40 c     [y0,nvs,[,hotdata]]=dasrt(y0,t0,t1[,atol,rtol],res[,jac],nh,h,info
41 c     [,hotdata])
42       ierror=0
43       maxord=5
44       lbuf = 1
45       topk=top
46       topw=top+1
47       lw = lstk(topw)
48       l0 = lstk(top+1-rhs)
49       if (.not.checkrhs(fname,6,11)) return
50       if (.not.checklhs(fname,2,3)) return
51 c     checking variable y0 (number 1)
52 c     -------------------------------
53       ky=top-rhs+1
54       if(.not.getrmat(fname,topk,ky,n1,m1,l1))return
55       neq=n1
56       lydot=l1+n1
57       info(11)=0
58       if (m1 .eq.1) then
59          if (.not.cremat(fname,topw,0,n1,1,lydot,lc)) return
60          topw=topw+1
61          info(11)=1
62          call dset(n1,0.0d0,stk(lydot),1)
63       elseif(m1.ne.2) then
64          err = 1
65          call error(89)
66          return
67       else 
68          il1 = iadr(lstk(top-rhs+1))
69          istk(il1+2)=1
70       endif
71 c     checking variable t0 (number 2)
72 c     ----------------------------
73       kt0=top-rhs+2
74       if(.not.getscalar(fname,topk,kt0,lr2))return
75       t0=stk(lr2)
76 c     checking variable t1 (number 3)
77 c     -------------------------------
78       if(.not.getrmat(fname,topk,top-rhs+3,m3,n3,l3))return
79       nt=m3*n3
80 c     
81 c     checking variable atol (number 4)
82 c     --------------------------------
83       iskip=0
84       itype = gettype(top-rhs+4)
85       if ( itype .ne. 1) then
86          if (.not.cremat(fname,topw,0,1,1,latol,lc)) return
87          topw=topw+1
88          if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
89          topw=topw+1
90          stk(latol)=atol
91          stk(lrtol)=rtol
92          info(2)=0
93          iskip=iskip+2
94          goto 1105
95       endif
96       if(.not.getrvect(fname,topk,top-rhs+4,m4,n4,latol))return
97       m4 = m4*n4
98 c     checking variable rtol (number 5)
99 c     --------------------------------
100       itype = gettype(top-rhs+5)
101       if (itype .ne. 1) then
102          if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
103          topw=topw+1
104          stk(lrtol)=lrtol
105          info(2)=0
106          iskip=iskip+1
107          goto 1105
108       endif
109       if(.not.getrvect(fname,topk,top-rhs+5,m5,n5,lrtol))return
110       m5 = m5*n5
111       if(m5.ne.m4.or.(m5.ne.1.and.m5.ne.neq)) then
112          call error(60)
113          return
114       endif
115       if(m5.eq.1) then
116          info(2)=0
117       else
118          info(2)=1
119       endif
120       
121 c     checking variable res (number 6)
122 c     --------------------------------
123  1105 kres=top-rhs+6-iskip
124       if (.not.getexternal(fname,topk,kres,namer,type,
125      $     setfresd)) return
126
127 c     checking variable number 7
128 c     -----------------------------
129       kjac=top-rhs+7-iskip
130       if(kjac.gt.top) then
131          iskip=iskip+1
132          info(5)=0
133       else
134          is7 = gettype(kjac)
135          if(is7.eq.15) then
136 c     .     info or jac ? get list list first element type to decide
137             il7=iadr(lstk(kjac))
138             if (istk(il7).lt.0)  il7=istk(il7+1)
139             nelt=istk(il7+1)
140             l71=sadr(il7+3+nelt)
141             if (abs(istk(iadr(l71))).eq.11.or.
142      $           abs(istk(iadr(l71))).eq.13) then
143 c     .        jac
144                is7=istk(iadr(l71))
145             endif
146          endif
147          if((is7.ne.10).and.(is7.ne.11).and.(is7.ne.13)) then
148             iskip=iskip+1
149             info(5)=0
150          else
151             info(5)=1
152             if (.not.getexternal(fname,topk,kjac,namej,type,
153      $           setfjacd)) return
154          endif
155       endif
156 c     DASRT nh,h
157 c     checking variable number 8
158 c     -----------------------------
159       if(.not.getscalar(fname,topk,top-rhs+8-iskip,lr8))return
160       nh=int(stk(lr8))
161 c     checking variable number 9
162       ksurf=top-rhs+9-iskip
163       if (.not.getexternal(fname,topk,ksurf,names,type,
164      $        setfsurfd)) return
165 c     
166 c     checking variable info (number 10)
167 c     ------------------------------------
168       kinfo = top-rhs+10-iskip
169       if (kinfo.gt.top) then
170          info(4)=0
171          info(3)=0
172          info(6)=0
173          info(7)=0
174          info(8)=0
175          info(10)=0
176          info(11)=0
177          iskip=iskip+1
178          goto 10
179       endif
180       il10 = iadr(lstk(top-rhs+10-iskip))
181       if (istk(il10) .ne. 15) then
182 c     default info values
183          info(4)=0
184          info(3)=0
185          info(6)=0
186          info(7)=0
187          info(8)=0
188          info(10)=0
189          info(11)=0
190          iskip=iskip+1
191          goto 10
192       endif
193       n10=istk(il10+1)
194       l10=sadr(il10+n10+3)
195 c     
196 c     --   subvariable tstop(info) --
197       il10e1=iadr(l10+istk(il10+1+1)-1)
198       l10e1 = sadr(il10e1+4)
199       m10e1 = istk(il10e1+1)*istk(il10e1+2)
200       if(m10e1.eq.0) then
201          info(4)=0
202       else
203          info(4)=1
204          tstop=stk(l10e1)
205       endif
206       
207 c     
208 c     --   subvariable imode(info) --
209       il10e2=iadr(l10+istk(il10+1+2)-1)
210       l10e2 = sadr(il10e2+4)
211       info(3)=stk(l10e2)
212       
213 c     
214 c     --   subvariable band(info) --
215       il10e3=iadr(l10+istk(il10+1+3)-1)
216       m10e3 =istk(il10e3+2)*istk(il10e3+2)
217       l10e3 = sadr(il10e3+4)
218       if(m10e3.eq.0) then
219          info(6)=0
220       elseif(m10e3.eq.2) then
221          info(6)=1
222          ml=stk(l10e3)
223          mu=stk(l10e3+1)
224       else
225          err=10-iskip
226          call error(89)
227          return
228       endif
229 c     
230 c     --   subvariable maxstep(info) --
231       il10e4=iadr(l10+istk(il10+1+4)-1)
232       m10e4 =istk(il10e4+2)*istk(il10e4+2)
233       l10e4 = sadr(il10e4+4)
234       if(m10e4.eq.0) then
235          info(7)=0
236       else
237          info(7)=1
238          maxstep=stk(l10e4)
239       endif
240       
241 c     
242 c     --   subvariable stepin(info) --
243       il10e5=iadr(l10+istk(il10+1+5)-1)
244       m10e5 =istk(il10e5+2)*istk(il10e5+2)
245       l10e5 = sadr(il10e5+4)
246       if(m10e5.eq.0) then
247          info(8)=0
248       else
249          info(8)=1
250          stepin=stk(l10e5)
251       endif
252       
253 c     
254 c     --   subvariable nonneg(info) --
255       il10e6=iadr(l10+istk(il10+1+6)-1)
256       l10e6 = sadr(il10e6+4)
257       info(10)=stk(l10e6)
258 c     
259 c     --   subvariable isest(info) --
260       il10e7=iadr(l10+istk(il10+1+7)-1)
261       l10e7 = sadr(il10e7+4)
262       isest=stk(l10e7)
263       if(isest.eq.1) info(11)=1
264       
265       
266  10   hotstart=.false.
267       if(rhs.eq.11-iskip) then
268          hotstart=.true.
269 c     
270 c     checking variable hotdata (number 11)
271 c     --------------------------------------
272          
273          il11 = iadr(lstk(top-rhs+11-iskip))
274          if (istk(il11) .ne. 1) then
275             err = 11-iskip
276             call error(53)
277             return
278          endif
279          n11 = istk(il11+1)*istk(il11+2)
280          lhot = sadr(il11+4)
281       elseif(rhs.ne.10-iskip) then
282          call error(39)
283          return
284       endif
285 c     --------------------Work Tables 
286       if (.not.cremat(fname,topw,0,1,1,lw15,lc)) return
287       topw=topw+1      
288       if (.not.cremat(fname,topw,0,1,1,lw17,lc)) return
289       topw=topw+1      
290       il17=iadr(lw17)
291 c     dasrt needs more
292       if (.not.cremat(fname,topw,0,1,nh,lgr,lc)) return
293       topw=topw+1      
294       lgroot=iadr(lgr)
295 c     
296       if(info(6).eq.0) then
297 C     for the full (dense) JACOBIAN case 
298          lrw = 50+(maxord+4)*neq+neq**2+3*nh
299       elseif(info(5).eq.1) then
300 C     for the banded user-defined JACOBIAN case
301          lrw=50+(maxord+4)*neq+(2*ml+mu+1)*neq+3*nh
302       elseif(info(5).eq.0) then
303 C     for the banded finite-difference-generated JACOBIAN case
304          lrw = 50+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1)+
305      $        3*nh
306       endif
307       liw=20+neq
308       if(.not.hotstart) then
309          if (.not.cremat(fname,topw,0,1,lrw,lrwork,lc)) return
310          topw=topw+1
311          if (.not.cremat(fname,topw,0,1,sadr(liw)+1,liwork,lc)) return
312          topw=topw+1
313       else
314          if(lrw+liw.gt.n11) then
315             err=11-iskip
316             call error(89)
317             return
318          endif
319          lrwork=lhot
320          liwork=lhot+lrw
321          call entier(liw,stk(liwork),istk(iadr(liwork)))
322       endif
323 c     
324       if(info(4).eq.1) then
325          stk(lrwork)=tstop
326       endif
327       if(info(7).eq.1) then
328          stk(lrwork+1)=maxstep
329       endif
330       if(info(8).eq.1) then
331          stk(lrwork+2)=stepin
332       endif
333       if(info(6).eq.1) then
334          istk(iadr(liwork))=ml
335          istk(iadr(liwork+1))=mu
336       endif
337 c     structure d'info pour les externals
338       top=topw
339       lw=lstk(top)
340       ilext=iadr(lw)
341       istk(ilext)=3
342       istk(ilext+1)=ilext+5
343       istk(ilext+2)=ilext+9
344       istk(ilext+3)=ilext+13
345       istk(ilext+4)=ilext+16
346       istk(ilext+5)=kres
347       istk(ilext+6)=neq
348       istk(ilext+7)=kt0
349       istk(ilext+8)=ky
350       istk(ilext+9)=kjac
351       istk(ilext+10)=neq
352       istk(ilext+11)=kt0
353       istk(ilext+12)=ky
354       istk(ilext+13)=ksurf
355       istk(ilext+14)=kt0
356       istk(ilext+15)=ky
357 c     istk(ilext+16)=ky
358       lw=sadr(ilext)+16
359       
360       lw0=lw
361       ilyr=iadr(lw)
362       istk(ilyr)=1
363       istk(ilyr+1)=2*n1+1
364       istk(ilyr+3)=0
365       lyr=sadr(ilyr+4)
366       lyri=lyr-(2*n1+1)
367       k=0
368       info(1)=0
369       if(hotstart) info(1)=1
370       info(9)=0
371       do 1120 i=0,nt-1
372          tout=stk(l3+i)
373 c     
374  1115    k=k+1
375          lyri=lyri+(2*n1+1)
376          lw=lyri+(2*n1+1)
377          lstk(top+1)=lw
378          margin=(k-1)*(2*n1+1)+4
379          lw1=lw+margin
380          if(lhs.eq.3) lw1=lw1+4+lrw+liw
381          if(lw1-lstk(bot).gt.0) then
382 c     not enough memory
383             call msgstxt('Not enough memory to go further')
384             k=k-1
385             goto 1125
386          endif
387          if (tout .eq. t0) then
388             stk(lyri)=tout
389             call unsfdcopy(n1,stk(l1),1,stk(lyri+1),1)
390             call unsfdcopy(n1,stk(lydot),1,stk(lyri+n1+1),1)
391             l1=lyri+1
392             lydot=lyri+n1+1
393             t0=tout
394             goto 1120            
395          else
396             stk(lyri)=tout
397             call unsfdcopy(n1,stk(l1),1,stk(lyri+1),1)
398             call unsfdcopy(n1,stk(lydot),1,stk(lyri+n1+1),1)
399             l1=lyri+1
400             lydot=lyri+n1+1
401             call ddasrt(bresd,n1,t0,stk(l1),stk(lydot),
402      &           stk(lyri),info,stk(lrtol),stk(latol),idid,
403      &           stk(lrwork),lrw,istk(iadr(liwork)),liw,stk(lw15),
404      &           istk(il17),bjacd,bsurfd,nh,istk(lgroot))
405 C     SUBROUTINE DDASRT (RES,NEQ,T,Y,YPRIME,TOUT,
406 C     *  INFO,RTOL,ATOL,IDID,RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC,
407 C     *  G,NG,JROOT)
408          endif
409          if(err.gt.0.or.err1.gt.0)  return
410          if(idid.eq.1) then
411 C     A step was successfully taken in the intermediate-output mode. 
412 C     The code has not yet reached TOUT.
413             stk(lyri)=t0
414             info(1)=1
415             goto 1115
416             
417          elseif(idid.eq.2) then
418 C     The integration to TSTOP was successfully completed (T=TSTOP)
419             goto 1125
420             
421          elseif(idid.eq.3) then
422 C     The integration to TOUT was successfully completed (T=TOUT) by 
423 C     stepping past TOUT. Y,ydot are obtained by interpolation.
424             t0=tout
425             info(1)=1
426             goto 1120
427          elseif(idid.eq.4) then
428 C     one or more root found
429             stk(lyri)=t0
430 C     stk(lrw+41)
431             goto 1125 
432          elseif(idid.eq.-1) then
433 C     A large amount of work has been expended (About 500 steps)
434             call msgstxt('Too many steps necessary to reach next '//
435      &           'required time discretization point')
436             call msgstxt('Change discretization of time vector t '//
437      &           'or decrease accuracy')
438             stk(lyri)=t0
439             goto 1125
440          elseif(idid.eq.-2) then
441 C     The error tolerances are too stringent.
442             t0=tout
443             info(1)=1
444             goto 1115
445 c     buf='The error tolerances are too stringent'
446 c     call error(9982)
447 c     return
448          elseif(idid.eq.-3) then
449 C     The local error test cannot be satisfied because you specified 
450 C     a zero component in ATOL and the corresponding computed solution
451 C     component is zero. Thus, a pure relative error test is impossible 
452 C     for this component.
453             buf='atol and computed test value are zero'
454             call error(9983)
455             return
456          elseif(idid.eq.-6) then
457 C     repeated error test failures on the last attempted step.
458             call msgstxt('A singularity in the solution '//
459      &           'may be present')
460             goto 1125
461          elseif(idid.eq.-7) then
462 C     The corrector could not converge.
463             call msgstxt('May be inaccurate or ill-conditioned '//
464      &           'JACOBIAN')
465             goto 1125
466          elseif(idid.eq.-8) then
467 C     The matrix of partial derivatives is singular.
468             buf='The matrix of partial derivatives is singular'//
469      &           'Some of your equations may be redundant'
470             call error(9986)
471             return
472          elseif(idid.eq.-9) then
473 C     The corrector could not converge. there were repeated error 
474 c     test failures in this step.
475             call msgstxt('Either ill-posed problem or '//
476      &           'discontinuity or singularity encountered')
477             goto 1125
478          elseif(idid.eq.-10) then
479             call msgstxt('external ''res'' return many times'//
480      &           'with ires=-1')
481             goto 1125
482          elseif(idid.eq.-11) then
483 C     IRES equal to -2 was encountered  and control is being returned to the
484 C     calling program.
485             buf='error in external ''res'' '
486             call error(9989)
487             return
488          elseif(idid.eq.-12) then
489 C     DDASSL failed to compute the initial YPRIME.
490             buf='dassrt failed to compute the initial Ydot.'
491             call error(9990)
492             return
493          elseif(idid.eq.-33) then
494 C     The code has encountered trouble from which
495 C     it cannot recover. A message is printed
496 C     explaining the trouble and control is returned
497 C     to the calling program. For example, this occurs
498 C     when invalid input is detected.
499             call msgstxt('dassrt encountered trouble')
500             goto 1125
501          endif
502          t0=tout
503          info(1)=1
504  1120 continue
505 c     
506  1125 top=topk-rhs
507       mv=lw0-l0
508 c     
509 c     Variable de sortie: y0
510 c     
511       top=top+1
512       if(k.eq.0) istk(ilyr+1)=0
513       istk(ilyr+2)=k
514       lw=lyr+(2*n1+1)*k
515       lstk(top+1)=lw-mv
516 c     
517 c     Variable de sortie: roots
518 c     
519       top=top+1
520       ilw=iadr(lw)
521       err=lw+4+nh+1-lstk(bot)
522       if (err .gt. 0) then
523          call error(17)
524          return
525       endif
526       istk(ilw)=1
527       istk(ilw+1)=1
528       istk(ilw+2)=1
529       istk(ilw+3)=0
530       l=sadr(ilw+4)
531       stk(l)=t0
532       kkk=1
533       do 1153 i=0,nh-1
534          if(istk(lgroot+i).ne.0) then
535             l=l+1
536             kkk=kkk+1
537             istk(ilw+2)=istk(ilw+2)+1
538             stk(l)=i+1
539          endif
540  1153 continue
541       lw=l+1
542       lstk(top+1)=lw-mv
543       if(lhs.eq.2) goto 1150
544 c     
545 c     Variable de sortie: rwork
546 c     
547       top=top+1
548       ilw=iadr(lw)
549       err=lw+4+lrw+liw-lstk(bot)
550       if (err .gt. 0) then
551          call error(17)
552          return
553       endif
554       istk(ilw)=1
555       istk(ilw+1)=lrw+liw
556       istk(ilw+2)=1
557       istk(ilw+3)=0
558       lw=sadr(ilw+4)
559       call unsfdcopy(lrw,stk(lrwork),1,stk(lw),1)
560       call int2db(liw,istk(iadr(liwork)),1,stk(lw+lrw),1)
561       lw=lw+lrw+liw
562       lstk(top+1)=lw-mv
563 c     
564 c     Remise en place de la pile
565  1150 call unsfdcopy(lw-lw0,stk(lw0),1,stk(l0),1)      
566       return
567       end       
568