← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/scicos/cossim.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/scicos/cossim.f
1
 
      subroutine cossim(neq,x,xptr,z,zptr,iz,izptr,told,tf,
2
 
c     Copyright INRIA
3
 
     $     tevts,evtspt,nevts,pointi,inpptr,inplnk,outptr,
4
 
     $     outlnk,lnkptr,clkptr,ordptr,nptr,
5
 
     $     ordclk,nordcl,ztyp,cord,iord,niord,oord,zord,
6
 
     $     critev,rpar,rpptr,ipar,
7
 
     $     ipptr,funptr,funtyp,rhot,ihot,outtb,jroot,w,iwa,ierr) 
8
 
C     
9
 
C     
10
 
C..   Parameters .. 
11
 
c     maximum number of clock output for one block
12
 
      integer nts
13
 
      parameter (nts=100)
14
 
C     
15
 
      integer neq(*)
16
 
C     neq must contain after #states all integer data for simblk and grblk
17
 
      double precision x(*),z(*),told,tf,tevts(*),rpar(*),outtb(*)
18
 
      double precision w(*),rhot(*)
19
 
      integer iwa(*)
20
 
C     X must contain after state values all real data for simblk and grblk
21
 
      integer xptr(*),zptr(*),iz(*),izptr(*),evtspt(nevts),nevts,pointi
22
 
      integer inpptr(*),inplnk(*),outptr(*),outlnk(*),lnkptr(*)
23
 
      integer clkptr(*),ordptr(nptr),nptr,ztyp(*)
24
 
      integer ordclk(nordcl,2),nordcl,cord(*),iord(*),oord(*),zord(*)
25
 
      integer critev(*),rpptr(*),ipar(*),ipptr(*),funptr(*),funtyp(*)
26
 
      integer ihot(*),jroot(*),ierr
27
 
c
28
 
      logical hot,stuck
29
 
      integer i,k,ierr1,iopt,istate,itask,j,jdum,jt,
30
 
     &     ksz,flag,keve,kpo,nord,nclock
31
 
      double precision t
32
 
      double precision tvec(nts)
33
 
c
34
 
      external grblk,simblk
35
 
 
36
 
      integer otimer,ntimer,stimer
37
 
      external stimer
38
 
 
39
 
      integer         nblk,nordptr,nout,ng,nrwp,niwp,ncord,
40
 
     &     noord,nzord
41
 
      common /cossiz/ nblk,nordptr,nout,ng,nrwp,niwp,ncord,
42
 
     &     noord,nzord
43
 
C     
44
 
      integer halt
45
 
      common /coshlt/ halt
46
 
c
47
 
      integer kfun
48
 
      common /curblk/ kfun
49
 
c
50
 
      double precision scale
51
 
      common /rtfactor/ scale  
52
 
c     
53
 
      double precision atol,rtol,ttol,deltat
54
 
      common /costol/ atol,rtol,ttol,deltat
55
 
c
56
 
      save otimer
57
 
      data otimer/0/
58
 
c    
59
 
      ierr = 0
60
 
      hot = .false.
61
 
      stuck=.false.
62
 
      call xscion(inxsci)
63
 
C     initialization
64
 
      call iset(niwp,0,ihot,1)
65
 
      call dset(nrwp,0.0d0,rhot,1)
66
 
      call realtimeinit(told, scale)
67
 
 
68
 
      ntvec=0
69
 
 
70
 
c     initialisation (propagation of constant blocks outputs)
71
 
      if(niord.eq.0) goto 10
72
 
      do 05 jj=1,niord
73
 
         kfun=iord(jj)
74
 
         nclock = iord(jj+niord)
75
 
         flag=1
76
 
         call callf(kfun,nclock,funptr,funtyp,told,x,x,xptr,z,zptr,iz,
77
 
     $        izptr,rpar,rpptr,ipar,ipptr,tvec,ntvec,inpptr,inplnk
78
 
     $        ,outptr,outlnk,lnkptr,outtb,flag) 
79
 
         if (flag .lt. 0) then
80
 
            ierr = 5 - flag
81
 
            return
82
 
         endif
83
 
 05   continue
84
 
C     main loop on time
85
 
 10   continue
86
 
      if (told .ge. tf) return
87
 
      if (inxsci.eq.1) then
88
 
         ntimer=stimer()
89
 
         if (ntimer.ne.otimer) then
90
 
            call sxevents()
91
 
            otimer=ntimer
92
 
            if (halt.ne.0) then
93
 
               halt=0
94
 
               return
95
 
            endif
96
 
         endif
97
 
      endif
98
 
      if (pointi.eq.0) then
99
 
         t = tf
100
 
      else
101
 
         t = tevts(pointi)
102
 
      endif
103
 
      if (abs(t-told) .lt. ttol) then
104
 
         t = told
105
 
C     update output part
106
 
      endif
107
 
      if (told .gt. t) then
108
 
C     !  scheduling problem
109
 
         ierr = 1
110
 
         return
111
 
      endif
112
 
      if (told .ne. t) then
113
 
         if (xptr(nblk+1) .eq. 1) then
114
 
C     .     no continuous state
115
 
            if(told+deltat+ttol.gt.t) then
116
 
               told=t
117
 
            else
118
 
               told=told+deltat
119
 
            endif
120
 
c     .     update outputs of 'c' type blocks
121
 
            if (ncord.eq.0) goto 343
122
 
 
123
 
            call cdoit(neq,x,xptr,z,zptr,iz,izptr,told,tf
124
 
     $           ,tevts,evtspt,nevts,pointi,inpptr,inplnk,outptr
125
 
     $           ,outlnk,lnkptr,clkptr,ordptr,nptr
126
 
     $           ,ordclk,nordcl,cord,iord,niord,oord,zord,critev,
127
 
     $           rpar,rpptr,ipar
128
 
     $           ,ipptr,funptr,funtyp,outtb,w,hot,ierr) 
129
 
            if(ierr.ne.0) return
130
 
 343        continue
131
 
C     
132
 
         else
133
 
C     integrate
134
 
            if (hot) then
135
 
               istate = 2
136
 
            else
137
 
               istate = 1
138
 
            endif
139
 
            itask = 4
140
 
C     Compute tcrit (rhot(1))
141
 
            rhot(1)=tf+ttol
142
 
            kpo=pointi
143
 
 20         if(critev(kpo).eq.1) then
144
 
               rhot(1)=tevts(kpo)
145
 
               goto 30
146
 
            endif
147
 
            kpo=evtspt(kpo)
148
 
            if(kpo.ne.0) goto 20
149
 
 30         continue
150
 
c     
151
 
            
152
 
c     .     form initial zero crossing input signs
153
 
            ig=1
154
 
            if (ng.gt.0) then
155
 
c     .        loop on zero crossing block
156
 
               do 35 kfun=1,nblk
157
 
                  if (ztyp(kfun).eq.1) then
158
 
c     .           loop on block ports
159
 
                     do 34 kport=inpptr(kfun),inpptr(kfun+1)-1
160
 
                        klink=inplnk(kport)
161
 
                        do 33 i=lnkptr(klink),lnkptr(klink+1)-1
162
 
                           if (outtb(i).gt.0.d0) then
163
 
                              jroot(ng+ig) = 1
164
 
                           else
165
 
                              jroot(ng+ig) = 0
166
 
                           endif
167
 
                           ig=ig+1
168
 
 33                     continue
169
 
 34                  continue
170
 
                  endif
171
 
 35            continue
172
 
            endif
173
 
c     
174
 
            iopt = 0
175
 
c
176
 
            jt = 2
177
 
            t=min(told+deltat,min(t,tf+ttol))
178
 
c     
179
 
c
180
 
            call lsodar(simblk,neq,x,told,t,1,rtol,atol,itask,
181
 
     &           istate,iopt,rhot,nrwp,ihot,niwp,jdum,jt,grblk,
182
 
     &           ng,jroot)
183
 
            if (istate .le. 0) then
184
 
               if (istate .eq. -3) then
185
 
                  if(stuck) then
186
 
                     ierr= 2
187
 
                     return
188
 
                  endif
189
 
                  itask = 2
190
 
                  istate = 1
191
 
                  call lsoda(simblk,neq,x,told,t,
192
 
     &                 1,rtol,atol,itask,
193
 
     &                 istate,iopt,rhot,nrwp,ihot,niwp,jdum,jt)
194
 
                  hot = .false.
195
 
                  stuck=.true.
196
 
                  if (istate .gt. 0) goto 38
197
 
               endif
198
 
C     !        integration problem
199
 
               ierr = 100-istate
200
 
               return
201
 
            endif
202
 
            hot = .true.
203
 
            stuck=.false.
204
 
 38         continue
205
 
c     .     update outputs of 'c' type  blocks
206
 
            nclock = 0
207
 
            ntvec=0
208
 
            if (ncord.gt.0) then
209
 
               call cdoit(neq,x,xptr,z,zptr,iz,izptr,told,tf
210
 
     $              ,tevts,evtspt,nevts,pointi,inpptr,inplnk,outptr
211
 
     $              ,outlnk,lnkptr,clkptr,ordptr,nptr
212
 
     $              ,ordclk,nordcl,cord,iord,niord,oord,zord,critev,
213
 
     $              rpar,rpptr,ipar
214
 
     $              ,ipptr,funptr,funtyp,outtb,w,hot,ierr) 
215
 
               if(ierr.ne.0) return
216
 
            endif
217
 
            if (istate .eq. 3) then
218
 
C     .        at a least one root has been found
219
 
               ig = 1
220
 
               do 50 kfun = 1,nblk
221
 
                  if (ztyp(kfun).eq.1) then
222
 
c     .              loop on block input ports
223
 
                     ksz=0
224
 
                     do 42 kport=inpptr(kfun),inpptr(kfun+1)-1
225
 
c     .                 get corresponding link pointer 
226
 
                        klink=inplnk(kport)
227
 
                        ksz=ksz+lnkptr(klink+1)-lnkptr(klink)
228
 
 42                  continue
229
 
c     .           kev is a base 2 coding of reached zero crossing surfaces
230
 
                     kev=0
231
 
                     do 44 j = 1,ksz
232
 
                        kev=2*kev+jroot(ig+ksz-j)
233
 
 44                  continue
234
 
                     jjflg=1
235
 
                     if (kev.eq.0) jjflg=0
236
 
                     do 45 j = 1,ksz
237
 
                        kev=2*kev+jroot(ng+ig+ksz-j)
238
 
 45                  continue
239
 
                     ig=ig+ksz
240
 
                     if (jjflg .ne. 0) then
241
 
                        flag=3
242
 
                        ntvec=clkptr(kfun+1)-clkptr(kfun)
243
 
c     .              call corresponding block to determine output event (kev)
244
 
                        call callf(kfun,kev,funptr,funtyp,
245
 
     $                       told,x,x,xptr,z,
246
 
     $                       zptr,iz,izptr,rpar,rpptr,ipar,ipptr,tvec,
247
 
     $                       ntvec,inpptr,inplnk,outptr,outlnk,lnkptr,
248
 
     $                       outtb,flag) 
249
 
                        if(flag.lt.0) then
250
 
                           ierr=5-flag
251
 
                           return
252
 
                        endif
253
 
c     .              update event agenda
254
 
                        do 47 k=1,clkptr(kfun+1)-clkptr(kfun)
255
 
                           if (tvec(k).ge.told) then
256
 
                              if (critev(k+clkptr(kfun)-1).eq.1)
257
 
     $                             hot=.false.
258
 
                              call addevs(tevts,evtspt,nevts,pointi,
259
 
     &                             tvec(k),k+clkptr(kfun)-1,ierr1)
260
 
                              if (ierr1 .ne. 0) then
261
 
C     .                       nevts too small
262
 
                                 ierr = 3
263
 
                                 return
264
 
                              endif
265
 
                           endif
266
 
 47                     continue
267
 
                     endif
268
 
                  endif
269
 
 50            continue
270
 
            endif
271
 
c     !     save initial signs of zero crossing surface
272
 
c     
273
 
         endif
274
 
         call realtime(told)
275
 
      else
276
 
C     .  t==told
277
 
         call ddoit(neq,x,xptr,z,zptr,iz,izptr,told,tf,
278
 
     $        tevts,evtspt,nevts,pointi,inpptr,inplnk,outptr,
279
 
     $        outlnk,lnkptr,clkptr,ordptr,nptr,
280
 
     $        ordclk,nordcl,cord,iord,niord,oord,zord,critev,
281
 
     $        rpar,rpptr,ipar,
282
 
     $        ipptr,funptr,funtyp,outtb,w,iwa,hot,ierr) 
283
 
         if(ierr.ne.0) return
284
 
C     
285
 
      endif
286
 
C     end of main loop on time
287
 
      goto 10
288
 
      end
289
 
 
290