← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/cl-f2cl/intrepid

« back to all changes in this revision

Viewing changes to packages/odepack/opkdemo3.f

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2005-03-03 13:53:18 UTC
  • mfrom: (1.1.1 hoary)
  • Revision ID: james.westby@ubuntu.com-20050303135318-wpmxhzrts93qvs2o
Tags: 1.0+cvs.2005.03.03
New CVS release. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/odepack/opkdemo3.f
 
1
c-----------------------------------------------------------------------
 
2
c Demonstration program for the DLSODA package.
 
3
c This is the version of 14 June 2001.
 
4
c
 
5
c This version is in double precision.
 
6
c
 
7
c The package is used to solve two simple problems,
 
8
c one with a full Jacobian, the other with a banded Jacobian,
 
9
c with the 2 appropriate values of jt in each case.
 
10
c If the errors are too large, or other difficulty occurs,
 
11
c a warning message is printed.  All output is on unit lout = 6.
 
12
c-----------------------------------------------------------------------
 
13
      external f1, jac1, f2, jac2
 
14
      integer i, iopar, iopt, iout, istate, itask, itol, iwork,
 
15
     1   jt, leniw, lenrw, liw, lout, lrw, mband, mused,
 
16
     2   ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst
 
17
      double precision atol, dtout, dtout0, dtout1, er, erm, ero, hu,
 
18
     1     rtol, rwork, t, tout, tout1, tsw, y
 
19
      dimension y(25), rwork(522), iwork(45), neq(1), atol(1), rtol(1)
 
20
      data lout/6/, tout1/16.921743d0/, dtout/17.341162d0/
 
21
c
 
22
      nerr = 0
 
23
      itol = 1
 
24
      rtol(1) = 0.0d0
 
25
      atol(1) = 1.0d-8
 
26
      lrw = 522
 
27
      liw = 45
 
28
      iopt = 0
 
29
c
 
30
c First problem
 
31
c
 
32
      neq(1) = 2
 
33
      nout = 4
 
34
      write (lout,110) neq(1),itol,rtol(1),atol(1)
 
35
 110  format(/'Demonstration program for DLSODA package'////
 
36
     1  ' Problem 1:   Van der Pol oscillator:'/
 
37
     2  '              xdotdot - 20*(1 - x**2)*xdot + x = 0, ',
 
38
     3  '   x(0) = 2, xdot(0) = 0'/' neq =',i2/
 
39
     4  ' itol =',i3,'   rtol =',d10.1,'   atol =',d10.1//)
 
40
c
 
41
      do 190 jt = 1,2
 
42
      write (lout,120) jt
 
43
 120  format(//' Solution with jt =',i3//
 
44
     1       '  t               x               xdot       meth',
 
45
     2       '   nq     h           tsw'//)
 
46
      t = 0.0d0
 
47
      y(1) = 2.0d0
 
48
      y(2) = 0.0d0
 
49
      itask = 1
 
50
      istate = 1
 
51
      dtout0 = 0.5d0*tout1
 
52
      dtout1 = 0.5d0*dtout
 
53
      tout = dtout0
 
54
      ero = 0.0d0
 
55
      do 170 iout = 1,nout
 
56
        call dlsoda(f1,neq,y,t,tout,itol,rtol,atol,itask,istate,
 
57
     1              iopt,rwork,lrw,iwork,liw,jac1,jt)
 
58
        hu = rwork(11)
 
59
        tsw = rwork(15)
 
60
        nqu = iwork(14)
 
61
        mused = iwork(19)
 
62
        write (lout,140) t,y(1),y(2),mused,nqu,hu,tsw
 
63
 140    format(d12.5,d16.5,d14.3,2i6,2d13.3)
 
64
        if (istate .lt. 0) go to 175
 
65
        iopar = iout - 2*(iout/2)
 
66
        if (iopar .ne. 0) go to 160
 
67
        er = abs(y(1))
 
68
        ero = max(ero,er)
 
69
        if (er .gt. 1.0d-2) then
 
70
          write (lout,150)
 
71
 150      format(//' Warning: value at root exceeds 1.0d-2'//)
 
72
          nerr = nerr + 1
 
73
        endif
 
74
 160    if (iout .eq. 1) tout = tout + dtout0
 
75
        if (iout .gt. 1) tout = tout + dtout1
 
76
 170    continue
 
77
 175  continue
 
78
      if (istate .lt. 0) nerr = nerr + 1
 
79
      nst = iwork(11)
 
80
      nfe = iwork(12)
 
81
      nje = iwork(13)
 
82
      lenrw = iwork(17)
 
83
      leniw = iwork(18)
 
84
      nfea = nfe
 
85
      if (jt .eq. 2) nfea = nfe - neq(1)*nje
 
86
      write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
 
87
 180  format(//' Final statistics for this run:'/
 
88
     1  ' rwork size =',i4,'   iwork size =',i4/
 
89
     2  ' number of steps =',i5/
 
90
     3  ' number of f-s   =',i5/
 
91
     4  ' (excluding J-s) =',i5/
 
92
     5  ' number of J-s   =',i5/
 
93
     6  ' max. error at root =',d10.2)
 
94
 190  continue
 
95
c
 
96
c Second problem
 
97
c
 
98
      neq(1) = 25
 
99
      ml = 5
 
100
      mu = 0
 
101
      iwork(1) = ml
 
102
      iwork(2) = mu
 
103
      mband = ml + mu + 1
 
104
      atol(1) = 1.0d-6
 
105
      nout = 5
 
106
      write (lout,210) neq(1),ml,mu,itol,rtol(1),atol(1)
 
107
 210  format(///80('-')///
 
108
     1  ' Problem 2: ydot = A * y , where',
 
109
     2  '  A is a banded lower triangular matrix'/
 
110
     2  '            derived from 2-D advection PDE'/
 
111
     3  ' neq =',i3,'   ml =',i2,'   mu =',i2/
 
112
     4  ' itol =',i3,'   rtol =',d10.1,'   atol =',d10.1//)
 
113
      do 290 jt = 4,5
 
114
      write (lout,220) jt
 
115
 220  format(//' Solution with jt =',i3//
 
116
     1       '     t             max.err.     meth   ',
 
117
     2       'nq      h            tsw'//)
 
118
      t = 0.0d0
 
119
      do 230 i = 2,neq(1)
 
120
 230    y(i) = 0.0d0
 
121
      y(1) = 1.0d0
 
122
      itask = 1
 
123
      istate = 1
 
124
      tout = 0.01d0
 
125
      ero = 0.0d0
 
126
      do 270 iout = 1,nout
 
127
        call dlsoda(f2,neq,y,t,tout,itol,rtol,atol,itask,istate,
 
128
     1              iopt,rwork,lrw,iwork,liw,jac2,jt)
 
129
        call edit2(y,t,erm)
 
130
        hu = rwork(11)
 
131
        tsw = rwork(15)
 
132
        nqu = iwork(14)
 
133
        mused = iwork(19)
 
134
        write (lout,240) t,erm,mused,nqu,hu,tsw
 
135
 240    format(d15.5,d14.3,2i6,2d14.3)
 
136
        if (istate .lt. 0) go to 275
 
137
        er = erm/atol(1)
 
138
        ero = max(ero,er)
 
139
        if (er .gt. 1000.0d0) then
 
140
          write (lout,150)
 
141
          nerr = nerr + 1
 
142
        endif
 
143
 270    tout = tout*10.0d0
 
144
 275  continue
 
145
      if (istate .lt. 0) nerr = nerr + 1
 
146
      nst = iwork(11)
 
147
      nfe = iwork(12)
 
148
      nje = iwork(13)
 
149
      lenrw = iwork(17)
 
150
      leniw = iwork(18)
 
151
      nfea = nfe
 
152
      if (jt .eq. 5) nfea = nfe - mband*nje
 
153
      write (lout,280) lenrw,leniw,nst,nfe,nfea,nje,ero
 
154
 280  format(//' Final statistics for this run:'/
 
155
     1  ' rwork size =',i4,'   iwork size =',i4/
 
156
     2  ' number of steps =',i5/
 
157
     3  ' number of f-s   =',i5/
 
158
     4  ' (excluding J-s) =',i5/
 
159
     5  ' number of J-s   =',i5/
 
160
     6  ' error overrun =',d10.2)
 
161
 290  continue
 
162
      write (lout,300) nerr
 
163
 300  format(///' Number of errors encountered =',i3)
 
164
      stop
 
165
      end
 
166
 
 
167
      subroutine f1 (neq, t, y, ydot)
 
168
      integer neq
 
169
      double precision t, y, ydot
 
170
      dimension y(*), ydot(*)
 
171
      dimension neq(*)
 
172
      ydot(1) = y(2)
 
173
      ydot(2) = 20.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1)
 
174
      return
 
175
      end
 
176
 
 
177
      subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd)
 
178
      integer neq, ml, mu, nrowpd
 
179
      double precision t, y, pd
 
180
      dimension y(*), pd(nrowpd,*)
 
181
      dimension neq(*)
 
182
      pd(1,1) = 0.0d0
 
183
      pd(1,2) = 1.0d0
 
184
      pd(2,1) = -40.0d0*y(1)*y(2) - 1.0d0
 
185
      pd(2,2) = 20.0d0*(1.0d0 - y(1)*y(1))
 
186
      return
 
187
      end
 
188
 
 
189
      subroutine f2 (neq, t, y, ydot)
 
190
      integer neq, i, j, k, ng
 
191
      double precision t, y, ydot, alph1, alph2, d
 
192
      dimension y(*), ydot(*)
 
193
      dimension neq(*)
 
194
      data alph1/1.0d0/, alph2/1.0d0/, ng/5/
 
195
      do 10 j = 1,ng
 
196
      do 10 i = 1,ng
 
197
        k = i + (j - 1)*ng
 
198
        d = -2.0d0*y(k)
 
199
        if (i .ne. 1) d = d + y(k-1)*alph1
 
200
        if (j .ne. 1) d = d + y(k-ng)*alph2
 
201
 10     ydot(k) = d
 
202
      return
 
203
      end
 
204
 
 
205
      subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
 
206
      integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng
 
207
      double precision t, y, pd, alph1, alph2
 
208
      dimension y(*), pd(nrowpd,*)
 
209
      dimension neq(*)
 
210
      data alph1/1.0d0/, alph2/1.0d0/, ng/5/
 
211
      mband = ml + mu + 1
 
212
      mu1 = mu + 1
 
213
      mu2 = mu + 2
 
214
      do 10 j = 1,neq(1)
 
215
        pd(mu1,j) = -2.0d0
 
216
        pd(mu2,j) = alph1
 
217
 10     pd(mband,j) = alph2
 
218
      do 20 j = ng,neq(1),ng
 
219
 20     pd(mu2,j) = 0.0d0
 
220
      return
 
221
      end
 
222
 
 
223
      subroutine edit2 (y, t, erm)
 
224
      integer i, j, k, ng
 
225
      double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt
 
226
      dimension y(*)
 
227
      data alph1/1.0d0/, alph2/1.0d0/, ng/5/
 
228
      erm = 0.0d0
 
229
      if (t .eq. 0.0d0) return
 
230
      ex = 0.0d0
 
231
      if (t .le. 30.0d0) ex = exp(-2.0d0*t)
 
232
      a2 = 1.0d0
 
233
      do 60 j = 1,ng
 
234
        a1 = 1.0d0
 
235
        do 50 i = 1,ng
 
236
          k = i + (j - 1)*ng
 
237
          yt = t**(i+j-2)*ex*a1*a2
 
238
          er = abs(y(k)-yt)
 
239
          erm = max(erm,er)
 
240
          a1 = a1*alph1/i
 
241
 50       continue
 
242
        a2 = a2*alph2/j
 
243
 60     continue
 
244
      return
 
245
      end