fix import/export of ierode common on Windows
[scilab.git] / scilab / modules / differential_equations / sci_gateway / fortran / sci_f_daskr.f
1 c
2 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 c Copyright (C) Scilab Enterprises - 2013 - Paul Bignier
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.1-en.txt
10 c
11       subroutine daskri(fname)
12 c ====================================================================
13 C     daskr
14 c ====================================================================
15 c
16 c     This code was inspired by sci_f_dasrt.f from the same folder.
17 c     daskr() has more optional parameters than dasrt(), so the gateway
18 c has to adapt to the new options
19       INCLUDE 'stack.h'
20 c
21       character*(*) fname
22       character*(nlgh+1) namjac
23       common/cjac/namjac
24       integer iadr,sadr,gettype
25 c
26       double precision atol,rtol,t0
27       integer info(20),topk,topw
28       integer LENWP, LENIWP, maxl, kmp, nrmax, mxnit, mxnj, mxnh, lsoff
29       logical hotstart,type,getexternal,getrvect
30       logical checkrhs,checklhs,getrmat,cremat,getscalar
31       double precision tout,tstop,maxstep,stepin
32       double precision epli, stptol, epinit
33       character*(nlgh+1) namer,namej,names,namep,namepj
34       common /dassln/ namer,namej,names,namep,namepj
35       external bresd,bjacd,bsurfd,bpsold,bpjacd
36       external setfresd,setfjacd,setfsurfd,setfpsold,setfpjacd
37 c
38       data atol/1.d-7/,rtol/1.d-9/
39 c
40       iadr(l)=l+l-1
41       sadr(l)=(l/2)+1
42
43 c     SCILAB function : daskr
44 c     --------------------------
45 c     [y0,nvs[,hotdata]]=daskr(y0,t0,t1[,atol[,rtol]],res[,jac],nh,h[,info
46 c     [,psol][,pjac]][,hotdata])
47       ierror=0
48       maxord=5
49       lbuf = 1
50       topk=top
51       topw=top+1
52       lw = lstk(topw)
53       l0 = lstk(top+1-rhs)
54       if (.not.checkrhs(fname,6,13)) return
55       if (.not.checklhs(fname,2,3)) return
56 c     Checking variable y0 (number 1)
57 c     -------------------------------
58       ky=top-rhs+1
59       if(.not.getrmat(fname,topk,ky,n1,m1,l1))return
60       neq=n1
61       lydot=l1+n1
62       info(11)=0
63       if (m1 .eq.1) then
64          if (.not.cremat(fname,topw,0,n1,1,lydot,lc)) return
65          topw=topw+1
66          info(11)=1
67          call dset(n1,0.0d0,stk(lydot),1)
68       elseif(m1.ne.2) then
69          err = 1
70          call error(89)
71          return
72       else
73          il1 = iadr(lstk(top-rhs+1))
74          istk(il1+2)=1
75       endif
76 c     Checking variable t0 (number 2)
77 c     ----------------------------
78       kt0=top-rhs+2
79       if(.not.getscalar(fname,topk,kt0,lr2))return
80       t0=stk(lr2)
81 c     Checking variable t1 (number 3)
82 c     -------------------------------
83       if(.not.getrmat(fname,topk,top-rhs+3,m3,n3,l3))return
84       nt=m3*n3
85 c
86 c     Checking variable atol (number 4)
87 c     --------------------------------
88       iskip=0
89       itype = gettype(top-rhs+4)
90       if ( itype .ne. 1) then
91          if (.not.cremat(fname,topw,0,1,1,latol,lc)) return
92          topw=topw+1
93          if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
94          topw=topw+1
95          stk(latol)=atol
96          stk(lrtol)=rtol
97          info(2)=0
98          iskip=iskip+2
99          goto 1105
100       endif
101       if(.not.getrvect(fname,topk,top-rhs+4,m4,n4,latol))return
102       m4 = m4*n4
103 c     Checking variable rtol (number 5)
104 c     --------------------------------
105       itype = gettype(top-rhs+5)
106       if (itype .ne. 1) then
107          if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
108          topw=topw+1
109          stk(lrtol)=lrtol
110          info(2)=0
111          iskip=iskip+1
112          goto 1105
113       endif
114       if(.not.getrvect(fname,topk,top-rhs+5,m5,n5,lrtol))return
115       m5 = m5*n5
116       if(m5.ne.m4.or.(m5.ne.1.and.m5.ne.neq)) then
117          call error(60)
118          return
119       endif
120       if(m5.eq.1) then
121          info(2)=0
122       else
123          info(2)=1
124       endif
125
126 c     Checking variable res (number 6)
127 c     --------------------------------
128  1105 kres=top-rhs+6-iskip
129       if (.not.getexternal(fname,topk,kres,namer,type,
130      $     setfresd)) return
131
132 c     Checking variable number 7
133 c     -----------------------------
134       kjac=top-rhs+7-iskip
135       if(kjac.gt.top) then
136          iskip=iskip+1
137          info(5)=0
138       else
139          is7 = gettype(kjac)
140          if(is7.eq.15) then
141 c     .     info or jac ? get first element type to decide
142             il7=iadr(lstk(kjac))
143             if (istk(il7).lt.0)  il7=istk(il7+1)
144             nelt=istk(il7+1)
145             l71=sadr(il7+3+nelt)
146             if (abs(istk(iadr(l71))).eq.11.or.
147      $           abs(istk(iadr(l71))).eq.13) then
148 c     .        jac
149                is7=istk(iadr(l71))
150             endif
151          endif
152          if((is7.ne.10).and.(is7.ne.11).and.(is7.ne.13)) then
153             iskip=iskip+1
154             info(5)=0
155          else
156             info(5)=1
157             if (.not.getexternal(fname,topk,kjac,namej,type,
158      $           setfjacd)) return
159          endif
160       endif
161 c     DASKR nh, h
162 c     Checking variable number 8
163 c     -----------------------------
164       if(.not.getscalar(fname,topk,top-rhs+8-iskip,lr8))return
165       nh=int(stk(lr8))
166 c     Checking variable number 9
167       ksurf=top-rhs+9-iskip
168       if (.not.getexternal(fname,topk,ksurf,names,type,
169      $        setfsurfd)) return
170 c
171 c     Checking variable info (number 10)
172 c     ------------------------------------
173       kinfo = top-rhs+10-iskip
174       if (kinfo.gt.top) then
175          info(4)=0
176          info(3)=0
177          info(6)=0
178          info(7)=0
179          info(8)=0
180          info(10)=0
181          info(11)=0
182          info(12)=0
183          info(13)=0
184          info(14)=0
185          info(15)=0
186          info(16)=0
187          info(17)=0
188          info(18)=0
189          iskip=iskip+1
190          goto 8
191       endif
192       il10 = iadr(lstk(top-rhs+10-iskip))
193       if (istk(il10) .ne. 15) then
194 c     Default info values
195          info(4)=0
196          info(3)=0
197          info(6)=0
198          info(7)=0
199          info(8)=0
200          info(10)=0
201          info(11)=0
202          info(12)=0
203          info(13)=0
204          info(14)=0
205          info(15)=0
206          info(16)=0
207          info(17)=0
208          info(18)=0
209          iskip=iskip+1
210          goto 8
211       endif
212       n10=istk(il10+1)
213       l10=sadr(il10+n10+3)
214 c
215 c     --   subvariable tstop(info) --
216       il10e1=iadr(l10+istk(il10+1+1)-1)
217       l10e1 = sadr(il10e1+4)
218       m10e1 = istk(il10e1+1)*istk(il10e1+2)
219       if(m10e1.eq.0) then
220          info(4)=0
221       else
222          info(4)=1
223          tstop=stk(l10e1)
224       endif
225
226 c
227 c     --   subvariable imode(info) --
228       il10e2=iadr(l10+istk(il10+1+2)-1)
229       l10e2 = sadr(il10e2+4)
230       info(3)=stk(l10e2)
231
232 c
233 c     --   subvariable band(info) --
234       il10e3=iadr(l10+istk(il10+1+3)-1)
235       m10e3 =istk(il10e3+2)*istk(il10e3+2)
236       l10e3 = sadr(il10e3+4)
237       if(m10e3.eq.0) then
238          info(6)=0
239       elseif(m10e3.eq.2) then
240          info(6)=1
241          ml=stk(l10e3)
242          mu=stk(l10e3+1)
243       else
244          err=10-iskip
245          call error(89)
246          return
247       endif
248
249 c
250 c     --   subvariable maxstep(info) --
251       il10e4=iadr(l10+istk(il10+1+4)-1)
252       m10e4 =istk(il10e4+2)*istk(il10e4+2)
253       l10e4 = sadr(il10e4+4)
254       if(m10e4.eq.0) then
255          info(7)=0
256       else
257          info(7)=1
258          maxstep=stk(l10e4)
259       endif
260
261 c
262 c     --   subvariable stepin(info) --
263       il10e5=iadr(l10+istk(il10+1+5)-1)
264       m10e5 =istk(il10e5+2)*istk(il10e5+2)
265       l10e5 = sadr(il10e5+4)
266       if(m10e5.eq.0) then
267          info(8)=0
268       else
269          info(8)=1
270          stepin=stk(l10e5)
271       endif
272
273 c
274 c     --   subvariable nonneg(info) --
275       il10e6=iadr(l10+istk(il10+1+6)-1)
276       l10e6 = sadr(il10e6+4)
277       info(10)=stk(l10e6)
278
279 c
280 c     --   subvariable consistent(info) --
281       il10e7=iadr(l10+istk(il10+1+7)-1)
282       m10e7 =istk(il10e7+2)*istk(il10e7+2)
283       l10e7 = sadr(il10e7+4)
284       if(m10e7.eq.0.or.(m10e7.eq.1.and.stk(l10e7).eq.0)) then
285 c        info is then [] or [0]
286          info(11)=0
287       else
288 c        info then looks like list(..., [+-1 +-1 ... +-1 +-1],...)
289          info(11)=1
290          if (info(10).eq.0.or.info(10).eq.2) then
291             LID = 40
292          else
293             LID = 40+neq
294          endif
295 c        Make no further changes here, but copy the [+-1] tab
296 c in 'iwork' once its size 'liw' is defined.
297       endif
298
299 c
300 c     --   subvariable iteration(info) --
301       il10e8=iadr(l10+istk(il10+1+8)-1)
302       l10e8 = sadr(il10e8+4)
303       iter=stk(l10e8)
304       if(iter.eq.1) then
305          info(12)=1
306       else
307          info(12)=0
308       endif
309
310 c
311 c     --   subvariable defaultKrylov(info) --
312       il10e9=iadr(l10+istk(il10+1+9)-1)
313       m10e9 =istk(il10e9+2)*istk(il10e9+2)
314       l10e9 = sadr(il10e9+4)
315       if(m10e9.eq.0) then
316          info(13)=0
317 c        maxl and kmp need default values
318          maxl  = min(5,neq)
319          kmp   = maxl
320       else
321 c     info then looks like list(..., [maxl kmp nrmax epli],...)
322          info(13)=1
323          maxl  = stk(l10e9)
324          kmp   = stk(l10e9+1)
325          nrmax = stk(l10e9+2)
326          epli  = stk(l10e9+3)
327       endif
328
329 c
330 c     --   subvariable justConsistentComp(info) --
331       il10e10=iadr(l10+istk(il10+1+10)-1)
332       l10e10 = sadr(il10e10+4)
333       proceed=stk(l10e10)
334       if(proceed.eq.0) then
335          info(14)=0
336       else
337 c        Check that info(11) = 1, meaning that the provided initial values
338 c        are not consistent
339          if (info(11).eq.1)  info(14) = 1
340       endif
341
342 c
343 c     --   subvariable psolJac(info) --
344       il10e11=iadr(l10+istk(il10+1+11)-1)
345       l10e11 = sadr(il10e11+4)
346       psolJac=stk(l10e11)
347       if (psolJac.eq.0) then
348          info(15)=0
349       else
350          info(15)=1
351       endif
352
353 c
354 c     --   subvariable excludeAlgebraic(info) --
355       il10e12=iadr(l10+istk(il10+1+12)-1)
356       m10e12 =istk(il10e12+2)*istk(il10e12+2)
357       l10e12 = sadr(il10e12+4)
358       if (m10e12.eq.0) then
359          info(16)=0
360       else
361 c     info then looks like list(..., [+-1 +-1 ... +-1 +-1],...)
362          info(16)=1
363          if (info(10).eq.0.or.info(10).eq.2) then
364             LID = 40
365          else
366             LID = 40+neq
367          endif
368 c        Make no further changes here, but copy the [+-1] tab
369 c in 'iwork' once its size 'liw' is defined.
370       endif
371
372 c
373 c     --   subvariable defaultHeuristic(info) --
374       il10e13=iadr(l10+istk(il10+1+13)-1)
375       m10e13 =istk(il10e13+2)*istk(il10e13+2)
376       l10e13 = sadr(il10e13+4)
377       if (m10e13.eq.0) then
378          info(17)=0
379       else
380 c  info then looks like list(..., [mxnit mxnj lsoff stptol epinit],...)
381          info(17)=1
382          mxnit  = stk(l10e9)
383          mxnj   = stk(l10e9+1)
384          mxnh   = stk(l10e9+2)
385          lsoff  = stk(l10e9+3)
386          stptol = stk(l10e9+4)
387          epinit = stk(l10e9+5)
388       endif
389
390 c
391 c     --   subvariable verbosity(info) --
392       il10e14=iadr(l10+istk(il10+1+14)-1)
393       l10e14 = sadr(il10e14+4)
394       verbosity=stk(l10e14)
395       if (verbosity.eq.1) then
396          info(18)=1
397       elseif(verbosity.eq.2) then
398          info(18)=2
399       else
400          info(18)=0
401       endif
402 c
403 c     Checking variable psol (number 11)
404 c     --------------------------------------
405 8      kpsol=top-rhs+11-iskip
406       if (kpsol.gt.top) then
407          iskip = iskip + 1
408          go to 9
409       endif
410       il11  = iadr(lstk(top-rhs+11-iskip))
411       pType = istk(il11)
412       if((pType.ne.10).and.(pType.ne.11).and.(pType.ne.13)) then
413          iskip = iskip + 1
414          go to 9
415       endif
416       if (.not.getexternal(fname,topk,kpsol,namep,type,
417      $        setfpsold)) return
418 c
419 c     Checking variable pjac (number 12)
420 c     --------------------------------------
421 9      kpjac=top-rhs+12-iskip
422       if (kpjac.gt.top) then
423          iskip = iskip + 1
424          go to 10
425       endif
426       il12  = iadr(lstk(top-rhs+12-iskip))
427       pType = istk(il12)
428       if((pType.ne.10).and.(pType.ne.11).and.(pType.ne.13)) then
429          iskip = iskip + 1
430          go to 10
431       endif
432       if (.not.getexternal(fname,topk,kpjac,namepj,type,
433      $        setfpjacd)) return
434
435  10   hotstart=.false.
436       if(rhs.eq.13-iskip) then
437          hotstart=.true.
438 c
439 c     Checking variable hotdata (number 13)
440 c     --------------------------------------
441          il13 = iadr(lstk(top-rhs+13-iskip))
442          if (istk(il13) .ne. 1) then
443             err = 13-iskip
444             call error(53)
445             return
446          endif
447          n13 = istk(il13+1)*istk(il13+2)
448          lhot = sadr(il13+4)
449       elseif(rhs.ne.12-iskip) then
450          call error(39)
451          return
452       endif
453
454 c     --------------------Work Tables
455       if (.not.cremat(fname,topw,0,1,2,lw15,lc)) return
456       topw=topw+1
457       if (.not.cremat(fname,topw,0,1,30,lw17,lc)) return
458       topw=topw+1
459       il17=iadr(lw17)
460 c     Set ipar and rpar to 0 by default
461       do 98 i=0,1
462          istk(iadr(lw15)+i) = 0
463 98    continue
464       do 99 i=0,29
465          istk(il17+i) = 0
466 99    continue
467 c     daskr needs more
468       if (.not.cremat(fname,topw,0,1,nh,lgr,lc)) return
469       topw=topw+1
470       lgroot=iadr(lgr)
471 c
472 c     base = lrw = 60 then augment size
473 c     according to the case (full dense, banded, ...)
474       base = 60
475       if (info(12).eq.0) then
476          lrw = base + max(maxord+4,7)*neq + 3*nh
477          if (info(6).eq.0) then
478 c        For the full (dense) JACOBIAN case
479             lrw = lrw + neq**2
480          elseif(info(5).eq.1) then
481 c        For the banded user-defined JACOBIAN case
482             lrw = lrw + (2*ml+mu+1)*neq
483          elseif(info(5).eq.0) then
484 c        For the banded finite-difference-generated JACOBIAN case
485             lrw = lrw + (2*ml+mu+1)*neq + 2*(neq/(ml+mu+1)+1)
486          endif
487       elseif(info(12).eq.1) then
488 c        LENWP is the length ot the rwork segment containing
489 c        the matrix elements of the preconditioner P
490          LENWP = neq*neq
491          lrw = base + (maxord+5)*neq + 3*nh
492      $                   + (maxl + 3 + min(1,maxl-kmp))*neq
493      $                   + (maxl+3)*maxl + 1 + LENWP
494       endif
495       if(info(16).eq.1) lrw = lrw + neq
496
497 c
498 c     base = liw = 40, then augment size according to the case
499       base = 40
500       if (info(12).eq.0) then
501          liw = base + neq
502       elseif(info(12).eq.1) then
503 c        LENIWP is the length ot the iwork segment containing
504 c        the matrix indexes of the preconditioner P
505 c        (compressed sparse row format)
506          LENIWP = 2*neq*neq
507          liw = base + LENIWP
508       endif
509       if(info(10).eq.1.or.info(10).eq.3) liw = liw + neq
510       if(info(11).eq.1.or.info(16).eq.1) liw = liw + neq
511       if(.not.hotstart) then
512          if (.not.cremat(fname,topw,0,1,lrw,lrwork,lc)) return
513          topw=topw+1
514          if (.not.cremat(fname,topw,0,1,sadr(liw)+1,liwork,lc)) return
515          topw=topw+1
516       else
517          if(lrw+liw.gt.n13) then
518 c           Test if output (/input) array 'hotdata' is large enough
519 c           to hold at least rwork and iwork
520             err=12-iskip
521             call error(89)
522             return
523          endif
524          lrwork=lhot
525          liwork=lhot+lrw
526          call entier(liw,stk(liwork),istk(iadr(liwork)))
527       endif
528 c
529       if(info(4).eq.1) then
530          stk(lrwork)=tstop
531       endif
532       if(info(7).eq.1) then
533          stk(lrwork+1)=maxstep
534       endif
535       if(info(8).eq.1) then
536          stk(lrwork+2)=stepin
537       endif
538       if(info(6).eq.1) then
539          istk(iadr(liwork))=ml
540          istk(iadr(liwork)+1)=mu
541       endif
542       if(info(11).eq.1) then
543          do 100 i=0,neq-1
544             istk(iadr(liwork)+LID+i) = stk(l10e7+i)
545 100      continue
546       endif
547       if(info(12).eq.1) then
548          istk(iadr(liwork)+26) = LENWP
549          istk(iadr(liwork)+27) = LENIWP
550       endif
551       if(info(13).eq.1) then
552          istk(iadr(liwork)+23) = maxl
553          istk(iadr(liwork)+24) = kmp
554          istk(iadr(liwork)+25) = nrmax
555          stk(lrwork+10)        = epli
556       endif
557       if (info(15).eq.1) then
558 c        Set ipar and rpar
559          istk(iadr(lw15)) = neq
560          istk(iadr(lw15)+1) = neq
561          istk(iadr(lw15)+2) = 2
562          istk(iadr(lw15)+3) = 2
563          istk(iadr(lw15)+4) = 1
564          istk(il17) = 0.005
565          istk(il17+1) = 0.05
566       endif
567       if(info(16).eq.1) then
568          do 101 i=0,neq-1
569             istk(iadr(liwork)+LID+i) = stk(l10e12+i)
570 101      continue
571       endif
572       if(info(17).eq.1) then
573          istk(iadr(liwork)+31) = mxnit
574          istk(iadr(liwork)+32) = mxnj
575          istk(iadr(liwork)+33) = mxnh
576          istk(iadr(liwork)+34) = lsoff
577          stk(lrwork+13)        = stptol
578          stk(lrwork+14)        = epinit
579       endif
580       istk(iadr(liwork)+16) = lrw
581       istk(iadr(liwork)+17) = liw
582
583 c     create header of wp
584       if (.not.cremat(fname,topw,0,neq*neq,1,lr,lc)) return
585       iheaderwp = topw
586       topw=topw+1
587 c     create header of wp and iwp
588       if (.not.cremat(fname,topw,0,neq*neq,2,lr,lc)) return
589       iheaderiwp = topw
590       topw=topw+1
591
592 c     Structure d'info pour les externals
593       top=topw
594       lw=lstk(top)
595       ilext=iadr(lw)
596       istk(ilext)=5
597       istk(ilext+1)=ilext+6
598       istk(ilext+2)=ilext+10
599       istk(ilext+3)=ilext+14
600       istk(ilext+4)=ilext+17
601       istk(ilext+5)=ilext+21
602       istk(ilext+6)=kres
603       istk(ilext+7)=neq
604       istk(ilext+8)=kt0
605       istk(ilext+9)=ky
606       istk(ilext+10)=kjac
607       istk(ilext+11)=neq
608       istk(ilext+12)=kt0
609       istk(ilext+13)=ky
610       istk(ilext+14)=ksurf
611       istk(ilext+15)=kt0
612       istk(ilext+16)=ky
613 c
614 c     function psol
615       istk(ilext+17)=kpsol
616 c     header used to create wp
617       istk(ilext+18)=iheaderwp
618 c     header used to create iwp
619       istk(ilext+19)=iheaderiwp
620 c     header of y used to create b
621       istk(ilext+20)=ky
622 c
623 c     function pjac
624       istk(ilext+21)=kpjac
625 c     header of t0 used to create neq, t0, h, cj
626       istk(ilext+22)=kt0
627 c     header of y used to create y, ydot, rewt and savr
628       istk(ilext+23)=ky
629       lw=sadr(ilext)+24
630
631       lw0=lw
632       ilyr=iadr(lw)
633       istk(ilyr)=1
634       istk(ilyr+1)=2*n1+1
635       istk(ilyr+3)=0
636       lyr=sadr(ilyr+4)
637       lyri=lyr-(2*n1+1)
638       k=0
639       info(1)=0
640       if(hotstart) info(1)=1
641       info(9)=0
642       do 1120 i=0,nt-1
643          tout=stk(l3+i)
644 c
645  1115    k=k+1
646          lyri=lyri+(2*n1+1)
647          lw=lyri+(2*n1+1)
648          lstk(top+1)=lw
649          margin=(k-1)*(2*n1+1)+4
650          lw1=lw+margin
651          if(lhs.eq.3) lw1=lw1+4+lrw+liw
652          if(lw1-lstk(bot).gt.0) then
653 c     Not enough memory
654             call msgstxt('Not enough memory to go further')
655             k=k-1
656             goto 1125
657          endif
658          if (tout .eq. t0) then
659             stk(lyri)=tout
660             call unsfdcopy(n1,stk(l1),1,stk(lyri+1),1)
661             call unsfdcopy(n1,stk(lydot),1,stk(lyri+n1+1),1)
662             l1=lyri+1
663             lydot=lyri+n1+1
664             t0=tout
665             goto 1120
666          else
667             stk(lyri)=tout
668             call unsfdcopy(n1,stk(l1),1,stk(lyri+1),1)
669             call unsfdcopy(n1,stk(lydot),1,stk(lyri+n1+1),1)
670             l1=lyri+1
671             lydot=lyri+n1+1
672             if (info(15).eq.1) then
673             call ddaskr(bresd, n1, t0, stk(l1), stk(lydot),
674      &           stk(lyri), info, stk(lrtol), stk(latol), idid,
675      &           stk(lrwork), lrw, istk(iadr(liwork)), liw, stk(lw15),
676      &           istk(il17), bpjacd, bpsold, bsurfd, nh, istk(lgroot))
677             else
678             call ddaskr(bresd, n1, t0, stk(l1), stk(lydot),
679      &           stk(lyri), info, stk(lrtol), stk(latol), idid,
680      &           stk(lrwork), lrw, istk(iadr(liwork)), liw, stk(lw15),
681      &           istk(il17), bjacd, bpsold, bsurfd, nh, istk(lgroot))
682             endif
683 C     SUBROUTINE DDASKR(RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
684 C     *   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL,
685 C     *   RT, NRT, JROOT)
686          endif
687          if(err.gt.0.or.err1.gt.0)  return
688          if(idid.eq.1) then
689 C     A step was successfully taken in the intermediate-output mode.
690 C     The code has not yet reached TOUT.
691             stk(lyri)=t0
692             info(1)=1
693             goto 1115
694
695          elseif(idid.eq.2.or.idid.eq.4) then
696 C     The integration to TSTOP was successfully completed (T=TSTOP)
697             goto 1125
698
699          elseif(idid.eq.3) then
700 C     The integration to TOUT was successfully completed (T=TOUT) by
701 C     stepping past TOUT. Y and Ydot are obtained by interpolation.
702             t0=tout
703             info(1)=1
704             goto 1120
705          elseif(idid.eq.5) then
706 C     One or more root found
707             stk(lyri)=t0
708 C     stk(lrw+41)
709             goto 1125
710          elseif(idid.eq.-1) then
711 C     A large amount of work has been expended (About 500 steps)
712             call msgstxt('Too many steps necessary to reach next '//
713      &           'required time discretization point.')
714             call msgstxt('Change discretization of time vector t '//
715      &           'or decrease accuracy.')
716             stk(lyri)=t0
717             goto 1125
718          elseif(idid.eq.-2) then
719 C     The error tolerances are too stringent.
720             t0=tout
721             info(1)=1
722             goto 1115
723 c     buf='The error tolerances are too stringent'
724 c     call error(9982)
725 c     return
726          elseif(idid.eq.-3) then
727 C     The local error test cannot be satisfied because you specified
728 C     a zero component in ATOL and the corresponding computed solution
729 C     component is zero. Thus, a pure relative error test is impossible
730 C     for this component.
731             buf='atol and computed test value are zero.'
732             call error(9983)
733             return
734          elseif(idid.eq.-5) then
735 C     There were repeated failures in the evaluation
736 C     or processing of the preconditioner (in JAC).
737             buf='Cannot evaluate the preconditioner.'
738             call error(9983)
739             return
740          elseif(idid.eq.-6) then
741 C     Repeated error test failures on the last attempted step.
742             call msgstxt('A singularity in the solution '//
743      &           'may be present.')
744             goto 1125
745          elseif(idid.eq.-7) then
746 C     The corrector could not converge.
747             call msgstxt('May be inaccurate or ill-conditioned '//
748      &           'JACOBIAN.')
749             goto 1125
750          elseif(idid.eq.-8) then
751 C     The matrix of partial derivatives is singular.
752             buf='The matrix of partial derivatives is singular'//
753      &           'Some of your equations may be redundant.'
754             call error(9986)
755             return
756          elseif(idid.eq.-9) then
757 C     The corrector could not converge. there were repeated error
758 c     test failures in this step.
759             call msgstxt('Either ill-posed problem or '//
760      &           'discontinuity or singularity encountered.')
761             goto 1125
762          elseif(idid.eq.-10) then
763             call msgstxt('External ''res'' returned many times'//
764      &           'with ires=-1.')
765             goto 1125
766          elseif(idid.eq.-11) then
767 C     IRES equal to -2 was encountered and control is being returned to the
768 C     calling program.
769             buf='Error in external ''res''.'
770             call error(9989)
771             return
772          elseif(idid.eq.-12) then
773 C     DDASKR failed to compute the initial Yprime.
774             buf='daskr failed to compute the initial Ydot.'
775             call error(9990)
776             return
777          elseif(idid.eq.-13) then
778 C     An unrecoverable error was encountered inside the user's PSOL routine,
779 C     and control is being returned to the calling program.
780             buf='Error in external ''Psol''.'
781             call error(9990)
782             return
783          elseif(idid.eq.-14) then
784 C     The Krylov linear system solver could not achieve convergence.
785             buf='The Krylov linear system solver did not converge.'
786             call error(9990)
787             return
788          elseif(idid.eq.-33) then
789 C     The code has encountered trouble from which
790 C     it cannot recover. A message is printed
791 C     explaining the trouble and control is returned
792 C     to the calling program. For example, this occurs
793 C     when invalid input is detected.
794             call msgstxt('daskr encountered trouble.')
795             goto 1125
796          endif
797          t0=tout
798          info(1)=1
799  1120 continue
800 c
801  1125 top=topk-rhs
802       mv=lw0-l0
803 c
804 c     Variable de sortie: y0
805 c
806       top=top+1
807       if(k.eq.0) istk(ilyr+1)=0
808       istk(ilyr+2)=k
809       lw=lyr+(2*n1+1)*k
810       lstk(top+1)=lw-mv
811 c
812 c     Variable de sortie: roots
813 c
814       top=top+1
815       ilw=iadr(lw)
816       err=lw+4+nh+1-lstk(bot)
817       if (err .gt. 0) then
818          call error(17)
819          return
820       endif
821       istk(ilw)=1
822       istk(ilw+1)=1
823       istk(ilw+2)=1
824       istk(ilw+3)=0
825       l=sadr(ilw+4)
826       stk(l)=t0
827       kkk=1
828       do 1153 i=0,nh-1
829          if(istk(lgroot+i).ne.0) then
830             l=l+1
831             kkk=kkk+1
832             istk(ilw+2)=istk(ilw+2)+1
833             stk(l)=i+1
834          endif
835  1153 continue
836       lw=l+1
837       lstk(top+1)=lw-mv
838       if(lhs.eq.2) goto 1150
839 c
840 c     Variable de sortie: rwork
841 c
842       top=top+1
843       ilw=iadr(lw)
844       err=lw+4+lrw+liw-lstk(bot)
845       if (err .gt. 0) then
846          call error(17)
847          return
848       endif
849       istk(ilw)=1
850       istk(ilw+1)=lrw+liw
851       istk(ilw+2)=1
852       istk(ilw+3)=0
853       lw=sadr(ilw+4)
854       call unsfdcopy(lrw,stk(lrwork),1,stk(lw),1)
855       call int2db(liw,istk(iadr(liwork)),1,stk(lw+lrw),1)
856       lw=lw+lrw+liw
857       lstk(top+1)=lw-mv
858 c
859 c     Remise en place de la pile
860  1150 call unsfdcopy(lw-lw0,stk(lw0),1,stk(l0),1)
861       return
862       end