Line data Source code
1 46 : subroutine derivk(x,y,dydx,ii,nn,iy,idy,karg)
2 : c
3 : c 2k-th order differentiation routine
4 : c ***********************************
5 : c
6 : c sets derivative of y(i,n) with respect to x into dydx(i,n),i=1,ii,n=1
7 : c id and idy are first dimensions of y and dydx respectively.
8 : c if x is found to be non-monotonic ii is set to 0.
9 : c
10 : c if derivk is called with karg .le. 0, it is assumed that a table
11 : c of differentiation coefficients has been set up.
12 : c
13 : c derivk uses as work space
14 : c wrk((2k+1)*(2k+1))
15 : c the size of wrk set in derivk is 100, which is sufficient, say, for
16 : c 8-th order differentiation in 20 dependent variables. if more
17 : c work space is needed it must be set in common/wrklir/
18 : c
19 : c note that this version of derivk also uses storage of size
20 : c
21 : c dc((2*k+1)*nn)
22 : c
23 : c in common /cderst/ to store differentiation coefficients. the
24 : c default size of 500 probably in general has to be increased in
25 : c the calling programme.
26 : c
27 : c original version: 18/3/1986
28 : c
29 : c Double precision version.
30 : c ++++++++++++++++++++++++
31 : c
32 : c Dated 10/3/90.
33 : c
34 : c Note: this double precision version of the routine has same name
35 : c as single precision version, but is distinguished by its file name
36 : c
37 : c Modified 19/8/05 allowing setting of simple derivative in case
38 : c of problems with leq. This is flagged by hardcoded iset_sim = 1.
39 : c With iset_sim = 0 use previous error exit with ii = 0.
40 : c
41 : implicit double precision(a-h,o-z)
42 : dimension x(1),y(iy,1),dydx(idy,1)
43 : common/wrklir/ wrk(100)
44 : common/cderst/ dc(500)
45 : c
46 : c common defining standard input and output
47 : c
48 : common/cstdio/ istdin, istdou, istdpr, istder
49 : c
50 : save /cderst/,k,k2
51 : c
52 46 : iset_sim=1
53 : c
54 : c test for setting up coefficients
55 : c
56 46 : if(karg.le.0) go to 60
57 : c
58 1 : k=karg
59 : c
60 : c check order
61 : c
62 1 : k2=2*k+1
63 1 : if(nn.lt.k2) then
64 0 : if(istdpr.gt.0) write(istdpr,100) k2,nn
65 0 : k=(nn-1)/2
66 0 : k2=2*k+1
67 : end if
68 : c
69 : c first dimension of wrk, storage information
70 : c
71 1 : iw=k2
72 1 : mst=iw*(iw-1)+2
73 : c
74 : c total range
75 : c
76 1 : diff=x(nn)-x(1)
77 1 : if(diff.eq.0) go to 95
78 : c
79 : c start loop for setting coefficients
80 : c
81 1 : 20 jdc=0
82 : c
83 1469 : do 50 n=1,nn
84 1468 : xtr=x(n)
85 : c
86 : c determine i
87 : c
88 1468 : i=min0(max0(0,n-1-k),nn-k2)
89 : c length of interval
90 1468 : dlx=x(i+k2)-x(i+1)
91 1468 : if(diff*dlx.le.0) go to 97
92 : c
93 : c set equations
94 : c
95 1468 : 32 dlxi=1.d0/dlx
96 : c
97 1468 : j1=jdc
98 1468 : m=0
99 8808 : do 40 l=1,k2
100 7340 : l1=i+l
101 7340 : j1=j1+1
102 7340 : xl=(x(l1)-xtr)*dlxi
103 7340 : aa=1.d0
104 7340 : m=m+1
105 7340 : wrk(m)=aa
106 36700 : do 35 ir=2,k2
107 29360 : m=m+1
108 29360 : aa=xl*aa
109 36700 : 35 wrk(m)=aa
110 8808 : 40 dc(j1)=0
111 : c
112 1468 : dc(jdc+2)=dlxi
113 : c
114 : c solve equations
115 : c
116 1468 : call leq(wrk,dc(jdc+1),k2,1,iw,1,det)
117 : c
118 1468 : if(det.eq.0) then
119 0 : i1=i+1
120 0 : ik2=i+k2
121 0 : if(istdpr.gt.0) write(istdpr,140) (l,x(l),l=i1,ik2)
122 0 : if(iset_sim.eq.0) then
123 0 : ii=0
124 0 : go to 98
125 : else
126 0 : if(istdpr.gt.0) write(istdpr,
127 0 : * '(/'' set simple derivative from parabolic fit''/)')
128 0 : do l=1,k2
129 0 : dc(jdc+l)=0.d0
130 : end do
131 0 : dxp=x(n+1)-x(n)
132 0 : dxm=x(n-1)-x(n)
133 0 : dc(jdc+k+1)=-(dxp+dxm)/(dxp*dxm)
134 0 : dc(jdc+k)=dxp/(dxm*(dxp-dxm))
135 0 : dc(jdc+k+2)=-dxm/(dxp*(dxp-dxm))
136 : end if
137 : end if
138 : c
139 1469 : 50 jdc=jdc+k2
140 : c
141 : c end setting coefficients
142 : c
143 1 : if(istdpr.gt.0) write(istdpr,110) jdc
144 : c
145 : c *****************************************
146 : c
147 : c set derivatives
148 : c
149 46 : 60 jdc=0
150 : c
151 67574 : do 70 n=1,nn
152 : c
153 : c determine i
154 : c
155 67528 : i=min0(max0(0,n-1-k),nn-k2)
156 : c
157 135056 : do 65 is=1,ii
158 67528 : j1=jdc
159 67528 : sum=0
160 405168 : do 62 l=1,k2
161 337640 : j1=j1+1
162 405168 : 62 sum=sum+dc(j1)*y(is,i+l)
163 : c
164 135056 : 65 dydx(is,n)=sum
165 : c
166 67574 : 70 jdc=jdc+k2
167 : c
168 46 : return
169 : c
170 : c diagnostics
171 : c
172 0 : 95 if(istdpr.gt.0) write(istdpr,120) x(1)
173 0 : ii=0
174 0 : return
175 0 : 97 i1=i+1
176 0 : ik2=i+k2
177 0 : if(istdpr.gt.0) write(istdpr,130) x(1),x(nn),i1,x(i1),ik2,x(ik2)
178 0 : ii=0
179 0 : return
180 : 98 continue
181 0 : return
182 : 100 format(1x,10('*'),' 2k+1 =',i4,' is greater than nn =',i4,
183 : * ' in derivk.'/11x,'k has been reset to (nn-1)/2')
184 : 110 format(//' storage needed in common/cderst/ in derivk:',
185 : * i6)
186 : 120 format(//1x,10('*'),' range is zero in derivk. x(1)=',
187 : * 1pe15.5,' = x(nn)'//)
188 : 130 format(//1x,10('*'),' independent variable is not monotonic',
189 : * ' in derivk'//' x(1)=',1pe15.5,' x(nn)=',e15.5,
190 : * 2(' x(',i4,')=',e15.5)//)
191 : 140 format(//1x,10('*'),' Singular matrix in derivk'/
192 : * ((' x(',i4,')=',1pe20.12)))
193 : end
|
