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