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