← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/interf/matimp.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/interf/matimp.f
 
1
      subroutine matimp
 
2
c ====================================================================
 
3
c     impl dassl dasrt : simulation  de systeme algebrico-differentiel
 
4
c ====================================================================
 
5
c     Copyright INRIA
 
6
      INCLUDE '../stack.h'
 
7
c     impl     dassl     dasrt
 
8
c     1         2          3
 
9
      if (ddt .eq. 4) then
 
10
         write(buf(1:4),'(i4)') fin
 
11
         call basout(io,wte,' matimp '//buf(1:4))
 
12
      endif
 
13
c     
 
14
      goto(10,100,1000) fin
 
15
      return
 
16
c     impl part 
 
17
 10   call sciimpl("impl")
 
18
      return
 
19
 100  call dassli("dassl")
 
20
      return
 
21
 1000 call dasrti("dasrt")
 
22
      return
 
23
      end
 
24
      
 
25
      subroutine sciimpl(fname)
 
26
c     ==================================================
 
27
      INCLUDE '../stack.h'
 
28
c
 
29
      character*(*) fname
 
30
      integer iadr,sadr
 
31
c
 
32
      double precision atol,rtol,t0,t1
 
33
      integer topk,topw
 
34
      logical jaco,achaud
 
35
      external bresid,badd,bj2
 
36
      external fres,fadda,fj2
 
37
      integer gettype
 
38
      logical getexternal,getrvect,vcopyobj
 
39
      logical checkrhs,checklhs,getrmat,cremat,getscalar
 
40
      logical typej,typea,typer,getsmat,vectsize
 
41
      character*(nlgh+1) namres,namadd,namjac
 
42
      character*1 strf
 
43
      common/cjac/namjac
 
44
      external setfres,setfadda,setfj2
 
45
      common/ierode/iero
 
46
c     
 
47
      data atol/1.d-7/,rtol/1.d-9/
 
48
c     
 
49
      iadr(l)=l+l-1
 
50
      sadr(l)=(l/2)+1
 
51
      if (.not.checkrhs(fname,6,10)) return
 
52
C     XXXXXX : pour l'instant 
 
53
      if (.not.checklhs(fname,1,3)) return
 
54
c     ---------------------------------
 
55
      iero=0
 
56
      topk=top
 
57
      topw=top+1
 
58
      iskip=1
 
59
      mf=20
 
60
c     first argument (check for string )
 
61
c     ----------------------------
 
62
      if(gettype(topk-rhs+1).eq.10) then
 
63
         iskip=0
 
64
         if(.not.getsmat(fname,topk,topk-rhs+1,
 
65
     $        m1,n1,1,1,lr1,nlr1))return
 
66
         call cvstr(1,istk(lr1),strf,1)
 
67
         if ( strf.eq.'a') mf = 10
 
68
         if ( strf.eq.'s') mf = 20
 
69
         if(strf.ne.'a'.and.strf.ne.'s') then
 
70
            call error(42)
 
71
            return
 
72
         endif
 
73
      endif
 
74
c     Initial condition : y0 arg 2 - iskip
 
75
c     -------------------
 
76
      kynew=topk-rhs+2-iskip
 
77
      if(.not.getrvect(fname,topk,kynew,ny,my,ly))return
 
78
      neq=ny*my
 
79
c     
 
80
c     Initial derivative condition : y0dot arg 3 - iskip
 
81
c     -------------------
 
82
      kydtop=topk-rhs+3-iskip
 
83
      if(.not.getrvect(fname,topk,kydtop,nyd,myd,lyd))return
 
84
      if(.not.vectsize(fname,topk,kydtop,ny*my)) return
 
85
c     t0 arg 4 - iskip
 
86
c     ----------------------------
 
87
      kttop=topk-rhs+4-iskip 
 
88
      if(.not.getscalar(fname,topk,kttop,lr4))return
 
89
      t0=stk(lr4)
 
90
c     t1 arg 5 - iskip
 
91
c     ---------------------------
 
92
      kt1=topk-rhs+5-iskip 
 
93
      if(.not.getrmat(fname,topk,kt1,m5,n5,lt5))return
 
94
      nn=m5*n5
 
95
c     checking variable rtol (number 6 - iskip)
 
96
c     -----------------------------------
 
97
      itype = gettype(topk-rhs+6-iskip)
 
98
      if ( itype .ne. 1) then
 
99
         if (.not.cremat(fname,topw,0,1,1,latol,lc)) return
 
100
         topw=topw+1
 
101
         if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
 
102
         topw=topw+1
 
103
         stk(latol)=atol
 
104
         stk(lrtol)=rtol
 
105
         na=1
 
106
         nr=1
 
107
         iskip = iskip+2
 
108
         goto 11105
 
109
      endif
 
110
      kr=top-rhs+6-iskip
 
111
      if(.not.getrvect(fname,topk,kr,m6,n6,lrtol))return
 
112
      nr = m6*n6
 
113
      if(nr.ne.1.and.nr.ne.neq) then
 
114
         err= 6-iskip 
 
115
         call error(89)
 
116
         return
 
117
      endif
 
118
 
 
119
c     checking variable rtol (number 7 - iskip )
 
120
c     --------------------------------
 
121
      itype = gettype(topk-rhs+7-iskip)
 
122
      if (itype .ne. 1) then
 
123
         if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
 
124
         topw=topw+1
 
125
         stk(lrtol)=rtol
 
126
         nr=1
 
127
         iskip = iskip +1
 
128
         goto 11105
 
129
      endif
 
130
      kr= top-rhs+7-iskip
 
131
      if(.not.getrvect(fname,topk,kr,m7,n7,latol))return
 
132
      na = m7*n7
 
133
      if(na.ne.1.and.na.ne.neq) then
 
134
         err= 7-iskip 
 
135
         call error(89)
 
136
         return
 
137
      endif
 
138
c     ----------------------------------
 
139
11105 if(nr.eq.1.and.na.eq.1) itol=1
 
140
      if(nr.eq.1.and.na.gt.1) itol=2
 
141
      if(nr.gt.1.and.na.eq.1) itol=3
 
142
      if(nr.gt.1.and.na.gt.1) itol=4
 
143
c     les externaux : res,adda et jac 
 
144
c     -----------------------------------
 
145
c     checking variable res (number 8 - iskip )
 
146
      kres=topk-rhs+8-iskip
 
147
      typer=.false.
 
148
      if (.not.getexternal(fname,topk,kres,namres,typer,
 
149
     $     setfres)) return
 
150
 
 
151
c     checking variable number 9 - iskip
 
152
c     -----------------------------
 
153
      kadd=topk-rhs+9-iskip
 
154
      if (.not.getexternal(fname,topk,kadd,namadd,typea,
 
155
     $     setfadda)) return
 
156
      if ( typea.neqv.typer) then 
 
157
         buf = fname // ': res and adda must have same type '
 
158
         call error(999)
 
159
      endif
 
160
c     checking variable number 10 - iskip
 
161
c     -----------------------------
 
162
      achaud=gettype(topk).eq.1
 
163
      kjac=topk-rhs+10-iskip
 
164
      if ( kjac.eq.topk.or.(achaud.and.kjac.eq.topk-2)) then 
 
165
         if (.not.getexternal(fname,topk,kjac,namjac,typej,
 
166
     $        setfj2)) return
 
167
         mf = mf+1
 
168
         jaco=.true.
 
169
      else
 
170
         jaco=.false.
 
171
         typej=.false.
 
172
         mf = mf + 2
 
173
      endif
 
174
 
 
175
c     other parameters 
 
176
c     -----------------
 
177
      itask=1
 
178
      istate=1
 
179
      iopt=0
 
180
c     hot restart case 
 
181
c     hot restart is detected when last argument is a matrix
 
182
c     ---------------------
 
183
C     space for result 
 
184
      if(.not.cremat(fname,topw,0,neq,nn,lres,lc)) return
 
185
      kresu=topw
 
186
      topw=topw+1
 
187
      nsizd=219
 
188
      nsizi=41
 
189
      if(achaud) then
 
190
c        iwork 
 
191
         if (.not.vcopyobj(fname,topk,topw)) return
 
192
         topw=topw+1
 
193
         if(.not.getrmat(fname,topk,topw-1,ml,nl,lci))return
 
194
         liwp= ml*nl
 
195
         ilc=iadr(lci)
 
196
         do 400 k= liwp -nsizi+1 , liwp
 
197
            write(06,*) k,'avant sauv iw',stk(lci+k-1)
 
198
 400     continue
 
199
c        rwork 
 
200
         if (.not.vcopyobj(fname,topk-1,topw)) return
 
201
         topw=topw+1
 
202
         if(.not.getrmat(fname,topk,topw-1,ml,nl,lrwp))return
 
203
         ilrw = ml*nl
 
204
         istate=2
 
205
c        restauration des commons
 
206
         lsavs=lrwp+ilrw- nsizd
 
207
         lsavi=lci+liwp- nsizi
 
208
         liwp1=liwp- nsizi
 
209
         call rscom1(stk(lsavs),stk(lsavi))
 
210
c        restauration du tableau entier
 
211
c        the end was used to restore the common 
 
212
         do 40 k=1,liwp
 
213
            istk(ilc+k-1)=int(stk(lci+k-1))
 
214
 40      continue
 
215
         
 
216
      else
 
217
c         ----create Work space 
 
218
         ilrw=0
 
219
         if(mf.gt.10) ilrw=22+16*neq+neq*neq
 
220
         if(mf.gt.20) ilrw=22+9*neq+neq*neq
 
221
         liwp=20+neq
 
222
         if(lhs.gt.1) then
 
223
            ilrw=ilrw+nsizd
 
224
            liwp=liwp+nsizi
 
225
         endif
 
226
         if(.not.cremat(fname,topw,0,1,liwp,li,lc)) return
 
227
         topw=topw+1
 
228
         ilc=iadr(li)
 
229
         do 1 k=1,liwp
 
230
            istk(ilc+k-1) =0
 
231
 1       continue
 
232
         if(.not.cremat(fname,topw,0,1,ilrw,lrwp,lc)) return
 
233
         do 11 k=1,liwp
 
234
            stk(lrwp+k-1) =0
 
235
 11      continue
 
236
         topw=topw+1
 
237
      endif
 
238
      if(jaco) then
 
239
         top=topw
 
240
         lw=lstk(top)
 
241
         ilw1=iadr(lw)
 
242
         istk(ilw1)=3
 
243
         istk(ilw1+1)=ilw1+4
 
244
         istk(ilw1+2)=ilw1+8
 
245
         istk(ilw1+3)=ilw1+12
 
246
         istk(ilw1+4)=kres
 
247
         istk(ilw1+5)=kttop
 
248
         istk(ilw1+6)=kynew
 
249
         istk(ilw1+7)=kydtop
 
250
         istk(ilw1+8)=kadd
 
251
         istk(ilw1+9)=kttop
 
252
         istk(ilw1+10)=kynew
 
253
         istk(ilw1+11)=kydtop
 
254
         istk(ilw1+12)=kjac
 
255
         istk(ilw1+13)=kttop
 
256
         istk(ilw1+14)=kynew
 
257
         istk(ilw1+15)=kydtop
 
258
         lstk(top+1)=sadr(ilw1+17)
 
259
      else
 
260
         top=topw
 
261
         lw=lstk(top)
 
262
         ilw1=iadr(lw)
 
263
         istk(ilw1)=2
 
264
         istk(ilw1+1)=ilw1+3
 
265
         istk(ilw1+2)=ilw1+7
 
266
         istk(ilw1+3)=kres
 
267
         istk(ilw1+4)=kttop
 
268
         istk(ilw1+5)=kynew
 
269
         istk(ilw1+6)=kydtop
 
270
         istk(ilw1+7)=kadd
 
271
         istk(ilw1+8)=kttop
 
272
         istk(ilw1+9)=kynew
 
273
         istk(ilw1+10)=kydtop
 
274
         lstk(top+1)=sadr(ilw1+11)
 
275
      endif
 
276
c     
 
277
      call xsetf(1)
 
278
      call xsetun(wte)
 
279
c     
 
280
c     appel de l'integrateur
 
281
      do 50 k=1,nn
 
282
         t1=stk(lt5 +k-1)
 
283
c        test sur le type des fonctions fournies
 
284
         if(typea) then
 
285
            if(typej) then
 
286
c     f fortran j fortran
 
287
               call lsodi(fres,fadda,fj2,neq,stk(ly),stk(lyd),t0,t1,
 
288
     1              itol,stk(lrtol),stk(latol),itask,istate,iopt,
 
289
     2              stk(lrwp),ilrw,istk(ilc),liwp,mf)
 
290
            else
 
291
c     f fortran j macro 
 
292
               call lsodi(fres,fadda,bj2,neq,stk(ly),stk(lyd),t0,t1,
 
293
     1              itol,stk(lrtol),stk(latol),itask,istate,iopt,
 
294
     2              stk(lrwp),ilrw,istk(ilc),liwp,mf)
 
295
            endif
 
296
         else
 
297
            if(typej) then
 
298
c     f macro j fortran
 
299
               call lsodi(bresid,badd,fj2,neq,stk(ly),stk(lyd),t0,t1,
 
300
     1              itol,stk(lrtol),stk(latol),itask,istate,iopt,
 
301
     2              stk(lrwp),ilrw,istk(ilc),liwp,mf)
 
302
            else
 
303
c     f macro j macro
 
304
               call lsodi(bresid,badd,bj2,neq,stk(ly),stk(lyd),t0,t1,
 
305
     1              itol,stk(lrtol),stk(latol),itask,istate,iopt,
 
306
     2              stk(lrwp),ilrw,istk(ilc),liwp,mf)
 
307
            endif
 
308
         endif
 
309
         if(err.gt.0) return
 
310
         if(istate.lt.0) then
 
311
            call error(24)
 
312
            return
 
313
         endif
 
314
         call unsfdcopy(neq,stk(ly),1,stk(lres+(k-1)*neq),1)
 
315
 50   continue
 
316
      top= topk-rhs+1
 
317
      call copyobj(fname,kresu,topk-rhs+1)
 
318
      if(lhs.eq.1) return
 
319
c     w
 
320
      if (lhs.ne.3) then 
 
321
         buf = fname // ' lhs can only be 1 or 2 '
 
322
         call error(999)
 
323
      endif
 
324
      top=top+1
 
325
      if (.not.cremat(fname,top,0,1,ilrw,lr,lc)) return
 
326
      call unsfdcopy(ilrw-nsizd,stk(lrwp),1,stk(lr),1)
 
327
      lsvs=lr+ilrw-nsizd
 
328
c     iw
 
329
      top=top+1
 
330
      if (.not.cremat(fname,top,0,1,liwp,lr,lc)) return
 
331
      do 60 k=1,liwp-nsizi
 
332
         stk(lr+k-1)=dble(istk(ilc+k-1))
 
333
 60   continue
 
334
      lsvi=lr+liwp-nsizi
 
335
      call svcom1(stk(lsvs),stk(lsvi))
 
336
      return
 
337
      end
 
338
 
 
339
      subroutine dassli(fname)
 
340
      character*(*) fname
 
341
c     ============================================
 
342
      INCLUDE '../stack.h'
 
343
c
 
344
      integer iadr,sadr
 
345
      integer topk,topw, info(15),gettype
 
346
      logical hotstart,getexternal,getrvect,type
 
347
      logical checkrhs,checklhs,getrmat,cremat,getscalar
 
348
      double precision tout,tstop,maxstep,stepin
 
349
      double precision atol,rtol,t0
 
350
      character*(nlgh+1) namer,namej,names
 
351
      character*(nlgh+1) namjac
 
352
      external bresd,bjacd
 
353
      external setfresd,setfjacd
 
354
      common /dassln/ namer,namej,names
 
355
      common/ierode/iero
 
356
      common/cjac/namjac
 
357
c     
 
358
      data atol/1.d-7/,rtol/1.d-9/
 
359
c     
 
360
      iadr(l)=l+l-1
 
361
      sadr(l)=(l/2)+1
 
362
c     
 
363
c     SCILAB function : dassl
 
364
c     --------------------------
 
365
c     [y0 [,hotdata]]=dassl(y0,t0,t1 [,atol,rtol],res [,jac],info..
 
366
c     [,hotdata])
 
367
      iero=0
 
368
      maxord=5
 
369
      lbuf = 1
 
370
      topk=top
 
371
      topw=top+1
 
372
      lw = lstk(topw)
 
373
      l0 = lstk(top+1-rhs)
 
374
      if (.not.checkrhs(fname,4,9)) return
 
375
      if (.not.checklhs(fname,1,2)) return
 
376
c     checking variable y0 (number 1)
 
377
c     -------------------------------
 
378
      ky=top-rhs+1
 
379
      if(.not.getrmat(fname,topk,ky,n1,m1,l1))return
 
380
      neq=n1
 
381
      lydot=l1+n1
 
382
      info(11)=0
 
383
      if (m1 .eq.1) then
 
384
         if (.not.cremat(fname,topw,0,n1,1,lydot,lc)) return
 
385
         topw=topw+1
 
386
         info(11)=1
 
387
         call dset(n1,0.0d0,stk(lydot),1)
 
388
      elseif(m1.ne.2) then
 
389
         err = 1
 
390
         call error(89)
 
391
         return
 
392
      else 
 
393
         il1 = iadr(lstk(top-rhs+1))
 
394
         istk(il1+2)=1
 
395
      endif
 
396
c     checking variable t0 (number 2)
 
397
c     -------------------------------
 
398
      kt0=top-rhs+2
 
399
      if(.not.getscalar(fname,topk,kt0,lr2))return
 
400
      t0=stk(lr2)
 
401
c     checking variable t1 (number 3)
 
402
c     -------------------------------
 
403
      if(.not.getrmat(fname,topk,top-rhs+3,m3,n3,l3))return
 
404
      nt=m3*n3
 
405
c     
 
406
c     checking variable atol (number 4)
 
407
c     -------------------------------
 
408
      iskip=0
 
409
      itype = gettype(top-rhs+4)
 
410
      if ( itype .ne. 1) then
 
411
         if (.not.cremat(fname,topw,0,1,1,latol,lc)) return
 
412
         topw=topw+1
 
413
         if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
 
414
         topw=topw+1
 
415
         stk(latol)=atol
 
416
         stk(lrtol)=rtol
 
417
         info(2)=0
 
418
         iskip=iskip+2
 
419
         goto 105
 
420
      endif
 
421
      if(.not.getrvect(fname,topk,top-rhs+4,m4,n4,latol))return
 
422
      m4 = m4*n4
 
423
c     checking variable rtol (number 5)
 
424
c     --------------------------------
 
425
      itype = gettype(top-rhs+5)
 
426
      if (itype .ne. 1) then
 
427
         if (.not.cremat(fname,topw,0,1,1,lrtol,lc)) return
 
428
         topw=topw+1
 
429
         stk(lrtol)=lrtol
 
430
         info(2)=0
 
431
         iskip=iskip+1
 
432
         goto 105
 
433
      endif
 
434
      if(.not.getrvect(fname,topk,top-rhs+5,m5,n5,lrtol))return
 
435
      m5 = m5*n5
 
436
      if(m5.ne.m4) then
 
437
         call error(60)
 
438
         return
 
439
      endif
 
440
      if(m5.eq.1) then
 
441
         info(2)=0
 
442
      else
 
443
         info(2)=1
 
444
      endif
 
445
c     checking variable res (number 6)
 
446
c     
 
447
 105  kres=top-rhs+6-iskip
 
448
      if (.not.getexternal(fname,topk,kres,namer,type,
 
449
     $     setfresd)) return
 
450
c     checking variable jac (number 7)
 
451
c     
 
452
      kjac=top-rhs+7-iskip
 
453
      if(kjac.gt.top) then
 
454
         iskip=iskip+1
 
455
         info(5)=0
 
456
      else
 
457
         isres=gettype(kjac)
 
458
         if(isres.eq.15) then
 
459
c     .     info or jac ? get list size to decide
 
460
            il6=iadr(lstk(kjac))
 
461
            if (istk(il6).lt.0)  il6=istk(il6+1)
 
462
            if (istk(il6+1).eq.2) then
 
463
c     .        jac
 
464
               isres=13
 
465
            endif
 
466
         endif
 
467
            
 
468
         if((isres.ne.10).and.(isres.ne.11).and.(isres.ne.13)) then
 
469
            iskip=iskip+1
 
470
            info(5)=0
 
471
         else
 
472
            info(5)=1
 
473
            if (.not.getexternal(fname,topk,kjac,namej,type,
 
474
     $           setfjacd)) return
 
475
         endif
 
476
      endif
 
477
c     
 
478
c     checking variable info (number 8)
 
479
c     ---------------------------------
 
480
      kinfo=top-rhs+8-iskip
 
481
      if (kinfo.gt.top) then
 
482
         info(4)=0
 
483
         info(3)=0
 
484
         info(6)=0
 
485
         info(7)=0
 
486
         info(8)=0
 
487
         info(10)=0
 
488
         info(11)=0
 
489
         iskip=iskip+1
 
490
         goto 10
 
491
      endif
 
492
      il8 = iadr(lstk(top-rhs+8-iskip))
 
493
      if (istk(il8) .ne. 15) then
 
494
c     default info values
 
495
         info(4)=0
 
496
         info(3)=0
 
497
         info(6)=0
 
498
         info(7)=0
 
499
         info(8)=0
 
500
         info(10)=0
 
501
         info(11)=0
 
502
         iskip=iskip+1
 
503
         goto 10
 
504
      endif
 
505
      n8=istk(il8+1)
 
506
      l8=sadr(il8+n8+3)
 
507
c     
 
508
c     --   subvariable tstop(info) --
 
509
      il8e1=iadr(l8+istk(il8+1+1)-1)
 
510
      l8e1 = sadr(il8e1+4)
 
511
      m8e1 = istk(il8e1+1)*istk(il8e1+2)
 
512
      if(m8e1.eq.0) then
 
513
         info(4)=0
 
514
      else
 
515
         info(4)=1
 
516
         tstop=stk(l8e1)
 
517
      endif
 
518
      
 
519
c     
 
520
c     --   subvariable imode(info) --
 
521
      il8e2=iadr(l8+istk(il8+1+2)-1)
 
522
      l8e2 = sadr(il8e2+4)
 
523
      info(3)=stk(l8e2)
 
524
      
 
525
c     
 
526
c     --   subvariable band(info) --
 
527
      il8e3=iadr(l8+istk(il8+1+3)-1)
 
528
      m8e3 =istk(il8e3+2)*istk(il8e3+2)
 
529
      l8e3 = sadr(il8e3+4)
 
530
      if(m8e3.eq.0) then
 
531
         info(6)=0
 
532
      elseif(m8e3.eq.2) then
 
533
         info(6)=1
 
534
         ml=stk(l8e3)
 
535
         mu=stk(l8e3+1)
 
536
      else
 
537
         err=8-iskip
 
538
         call error(89)
 
539
         return
 
540
      endif
 
541
c     
 
542
c     --   subvariable maxstep(info) --
 
543
      il8e4=iadr(l8+istk(il8+1+4)-1)
 
544
      m8e4 =istk(il8e4+2)*istk(il8e4+2)
 
545
      l8e4 = sadr(il8e4+4)
 
546
      if(m8e4.eq.0) then
 
547
         info(7)=0
 
548
      else
 
549
         info(7)=1
 
550
         maxstep=stk(l8e4)
 
551
      endif
 
552
      
 
553
c     
 
554
c     --   subvariable stepin(info) --
 
555
      il8e5=iadr(l8+istk(il8+1+5)-1)
 
556
      m8e5 =istk(il8e5+2)*istk(il8e5+2)
 
557
      l8e5 = sadr(il8e5+4)
 
558
      if(m8e5.eq.0) then
 
559
         info(8)=0
 
560
      else
 
561
         info(8)=1
 
562
         stepin=stk(l8e5)
 
563
      endif
 
564
      
 
565
c     
 
566
c     --   subvariable nonneg(info) --
 
567
      il8e6=iadr(l8+istk(il8+1+6)-1)
 
568
      l8e6 = sadr(il8e6+4)
 
569
      info(10)=stk(l8e6)
 
570
c     
 
571
c     --   subvariable isest(info) --
 
572
      il8e7=iadr(l8+istk(il8+1+7)-1)
 
573
      l8e7 = sadr(il8e7+4)
 
574
      isest=stk(l8e7)
 
575
      if(isest.eq.1) info(11)=1
 
576
      
 
577
      
 
578
 10   hotstart=.false.
 
579
      if(rhs.eq.9-iskip) then
 
580
         hotstart=.true.
 
581
c     
 
582
c     checking variable hotdata (number 9)
 
583
c     
 
584
         il9 = iadr(lstk(top-rhs+9-iskip))
 
585
         if (istk(il9) .ne. 1) then
 
586
            err = 9-iskip
 
587
            call error(53)
 
588
            return
 
589
         endif
 
590
         n9 = istk(il9+1)*istk(il9+2)
 
591
         lhot = sadr(il9+4)
 
592
      elseif(rhs.ne.8-iskip) then
 
593
         call error(39)
 
594
         return
 
595
      endif
 
596
c     
 
597
c     --------------------Work Tables 
 
598
      if (.not.cremat(fname,topw,0,1,1,lw15,lc)) return
 
599
      topw=topw+1      
 
600
      if (.not.cremat(fname,topw,0,1,1,lw17,lc)) return
 
601
      topw=topw+1      
 
602
      il17=iadr(lw17)
 
603
      if(info(6).eq.0) then
 
604
C     for the full (dense) JACOBIAN case 
 
605
         lrw = 40+(maxord+4)*neq+neq**2
 
606
      elseif(info(5).eq.1) then
 
607
C     for the banded user-defined JACOBIAN case
 
608
         lrw=40+(maxord+4)*neq+(2*ml+mu+1)*neq
 
609
      elseif(info(5).eq.0) then
 
610
C     for the banded finite-difference-generated JACOBIAN case
 
611
         lrw = 40+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1)
 
612
      endif
 
613
      liw=20+neq
 
614
      if(.not.hotstart) then
 
615
         if (.not.cremat(fname,topw,0,1,lrw,lrwork,lc)) return
 
616
         topw=topw+1
 
617
         if (.not.cremat(fname,topw,0,1,sadr(liw)+1,liwork,lc)) return
 
618
         topw=topw+1
 
619
      else
 
620
         if(lrw+liw.gt.n9) then
 
621
            err=9-iskip
 
622
            call error(89)
 
623
            return
 
624
         endif
 
625
         lrwork=lhot
 
626
         liwork=lhot+lrw
 
627
         call entier(liw,stk(liwork),istk(iadr(liwork)))
 
628
      endif
 
629
c     
 
630
      if(info(4).eq.1) then
 
631
         stk(lrwork)=tstop
 
632
      endif
 
633
      if(info(7).eq.1) then
 
634
         stk(lrwork+1)=maxstep
 
635
      endif
 
636
      if(info(8).eq.1) then
 
637
         stk(lrwork+2)=stepin
 
638
      endif
 
639
      if(info(6).eq.1) then
 
640
         istk(iadr(liwork))=ml
 
641
         istk(iadr(liwork+1))=mu
 
642
      endif
 
643
c     structure d'info pour les externals
 
644
      top=topw
 
645
      lw=lstk(top)
 
646
      ilext=iadr(lw)
 
647
      istk(ilext)=2
 
648
      istk(ilext+1)=ilext+4
 
649
      istk(ilext+2)=ilext+8
 
650
      istk(ilext+3)=ilext+12
 
651
      istk(ilext+4)=kres
 
652
      istk(ilext+5)=neq
 
653
      istk(ilext+6)=kt0
 
654
      istk(ilext+7)=ky
 
655
      istk(ilext+8)=kjac
 
656
      istk(ilext+9)=neq
 
657
      istk(ilext+10)=kt0
 
658
      istk(ilext+11)=ky
 
659
      lw=sadr(ilext)+12
 
660
      lw0=lw
 
661
      ilyr=iadr(lw)
 
662
      istk(ilyr)=1
 
663
      istk(ilyr+1)=2*n1+1
 
664
      istk(ilyr+3)=0
 
665
      lyr=sadr(ilyr+4)
 
666
      lyri=lyr-(2*n1+1)
 
667
      k=0
 
668
      info(1)=0
 
669
      if(hotstart) info(1)=1
 
670
      info(9)=0
 
671
      do 120 i=0,nt-1
 
672
         tout=stk(l3+i)
 
673
c     
 
674
 115     k=k+1
 
675
         lyri=lyri+(2*n1+1)
 
676
         lw=lyri+(2*n1+1)
 
677
         lstk(top+1)=lw
 
678
         margin=(k-1)*(2*n1+1)+4
 
679
         lw1=lw+margin
 
680
         if(lhs.eq.2) lw1=lw1+4+lrw+liw
 
681
         if(lw1-lstk(bot).gt.0) then
 
682
c     not enough memory
 
683
            call msgstxt('Not enough memory to go further')
 
684
            k=k-1
 
685
            goto 125
 
686
         endif
 
687
         if (tout .eq. t0) then
 
688
            stk(lyri)=tout
 
689
            call unsfdcopy(n1,stk(l1),1,stk(lyri+1),1)
 
690
            call unsfdcopy(n1,stk(lydot),1,stk(lyri+n1+1),1)
 
691
            l1=lyri+1
 
692
            lydot=lyri+n1+1
 
693
            t0=tout
 
694
            goto 120            
 
695
         else
 
696
            stk(lyri)=tout
 
697
            call unsfdcopy(n1,stk(l1),1,stk(lyri+1),1)
 
698
            call unsfdcopy(n1,stk(lydot),1,stk(lyri+n1+1),1)
 
699
            l1=lyri+1
 
700
            lydot=lyri+n1+1
 
701
            call ddassl(bresd,n1,t0,stk(l1),stk(lydot),
 
702
     &           stk(lyri),info,stk(lrtol),stk(latol),idid,
 
703
     &           stk(lrwork),lrw,istk(iadr(liwork)),liw,stk(lw15),
 
704
     &           istk(il17),bjacd)
 
705
         endif
 
706
         if(err.gt.0)  return
 
707
         if(idid.eq.1) then
 
708
C     A step was successfully taken in the intermediate-output mode. 
 
709
C     The code has not yet reached TOUT.
 
710
            stk(lyri)=t0
 
711
            info(1)=1
 
712
            goto 115
 
713
            
 
714
         elseif(idid.eq.2) then
 
715
C     The integration to TSTOP was successfully completed (T=TSTOP)
 
716
            goto 125
 
717
            
 
718
         elseif(idid.eq.3) then
 
719
C     The integration to TOUT was successfully completed (T=TOUT) by 
 
720
C     stepping past TOUT. Y,ydot are obtained by interpolation.
 
721
            t0=tout
 
722
            info(1)=1
 
723
            goto 120
 
724
            
 
725
         elseif(idid.eq.-1) then
 
726
C     A large amount of work has been expended (About 500 steps)
 
727
            call msgstxt('to many steps necessary to reached next '//
 
728
     &           'required time discretization point')
 
729
            call msgstxt('Change discretisation of time vector t '//
 
730
     &           'or decrease accuracy')
 
731
            stk(lyri)=t0
 
732
            goto 125
 
733
         elseif(idid.eq.-2) then
 
734
C     The error tolerances are too stringent.
 
735
            t0=tout
 
736
            info(1)=1
 
737
            goto 115
 
738
c     buf='The error tolerances are too stringent'
 
739
c     call error(9982)
 
740
c     return
 
741
         elseif(idid.eq.-3) then
 
742
C     The local error test cannot be satisfied because you specified 
 
743
C     a zero component in ATOL and the corresponding computed solution
 
744
C     component is zero. Thus, a pure relative error test is impossible 
 
745
C     for this component.
 
746
            buf='atol and computed test value are zero'
 
747
            call error(9983)
 
748
            return
 
749
         elseif(idid.eq.-6) then
 
750
C     repeated error test failures on the last attempted step.
 
751
            call msgstxt('A singularity in the solution '//
 
752
     &           'may be present')
 
753
            goto 125
 
754
         elseif(idid.eq.-7) then
 
755
C     The corrector could not converge.
 
756
            call msgstxt('May be inaccurate or ill-conditioned '//
 
757
     &           'JACOBIAN')
 
758
            goto 125
 
759
         elseif(idid.eq.-8) then
 
760
C     The matrix of partial derivatives is singular.
 
761
            buf='Singular partial derivatives matrix'//
 
762
     &           ' (may be redundant equations)'
 
763
            call error(9986)
 
764
            return
 
765
         elseif(idid.eq.-9) then
 
766
C     The corrector could not converge. there were repeated error 
 
767
c     test failures in this step.
 
768
            call msgstxt('Either ill-posed problem or '//
 
769
     &           'discontinuity or singularity encountered')
 
770
            goto 125
 
771
         elseif(idid.eq.-10) then
 
772
            call msgstxt('external ''res'' return many times'//
 
773
     &           ' with ires=-1')
 
774
            goto 125
 
775
         elseif(idid.eq.-11) then
 
776
C     IRES equal to -2 was encountered  and control is being returned to the
 
777
C     calling program.
 
778
            buf='error in external ''res'' '
 
779
            call error(9989)
 
780
            return
 
781
         elseif(idid.eq.-12) then
 
782
C     DDASSL failed to compute the initial YPRIME.
 
783
            buf='dassl failed to compute the initial Ydot.'
 
784
            call error(9990)
 
785
            return
 
786
         elseif(idid.eq.-33) then
 
787
C     The code has encountered trouble from which
 
788
C     it cannot recover. A message is printed
 
789
C     explaining the trouble and control is returned
 
790
C     to the calling program. For example, this occurs
 
791
C     when invalid input is detected.
 
792
            call msgstxt('dassl encountered trouble')
 
793
            goto 125
 
794
         endif
 
795
         t0=tout
 
796
         info(1)=1
 
797
 120  continue
 
798
c     
 
799
 125  top=topk-rhs
 
800
      mv=lw0-l0
 
801
c     
 
802
c     Variable de sortie: y0
 
803
c     
 
804
      top=top+1
 
805
      if(k.eq.0) istk(ilyr+1)=0
 
806
      istk(ilyr+2)=k
 
807
      lw=lyr+(2*n1+1)*k
 
808
      lstk(top+1)=lw-mv
 
809
      if(lhs.eq.1) goto 150
 
810
      
 
811
c     
 
812
c     Variable de sortie: rwork
 
813
c     
 
814
      top=top+1
 
815
      ilw=iadr(lw)
 
816
      err=lw+4+lrw+liw-lstk(bot)
 
817
      if (err .gt. 0) then
 
818
         call error(17)
 
819
         return
 
820
      endif
 
821
      istk(ilw)=1
 
822
      istk(ilw+1)=lrw+liw
 
823
      istk(ilw+2)=1
 
824
      istk(ilw+3)=0
 
825
      lw=sadr(ilw+4)
 
826
      call unsfdcopy(lrw,stk(lrwork),1,stk(lw),1)
 
827
      call int2db(liw,istk(iadr(liwork)),1,stk(lw+lrw),1)
 
828
      lw=lw+lrw+liw
 
829
      lstk(top+1)=lw-mv
 
830
c     
 
831
c     Remise en place de la pile
 
832
 150  call unsfdcopy(lw-lw0,stk(lw0),1,stk(l0),1)
 
833
      return
 
834
      end
 
835