1 SUBROUTINE doakima(s, f, nin, ier, case)
5 REAL(rprec),
DIMENSION(nin),
INTENT(in) :: s, f
6 INTEGER,
INTENT(in) :: nin
7 CHARACTER(len=*),
INTENT(in) :: case
9 INTEGER,
PARAMETER :: nset1 = 11, nset2 = 20
11 REAL(rprec),
DIMENSION(nin) :: s3
12 REAL(rprec),
DIMENSION(:),
ALLOCATABLE :: d1f, d2f, br
13 REAL(rprec),
DIMENSION(:,:),
ALLOCATABLE :: co
14 REAL,
DIMENSION(:),
ALLOCATABLE :: xplt,yplt
16 INTEGER,
DIMENSION(:),
ALLOCATABLE :: slctn, olli, span
17 INTEGER :: n, i, lastspan, tot, j, count, ier
19 REAL(rprec) :: fvn_spline_eval, fvn_spline_devaldx
20 REAL(rprec) :: fvn_spline_d2evaldx2
22 CHARACTER*99 xlabel, ylabel, slabel*60, mlabel*50
24 print*,
" Akima Spline for ",
case
28 ALLOCATE( br(n),stat = ier)
29 ALLOCATE( co(4,n),stat = ier)
30 ALLOCATE( slctn(n),stat = ier)
31 ALLOCATE(f2(n), s2(n),stat = ier)
33 CALL fvn_akima(n,s,f,br,co)
37 f2(i) = fvn_spline_eval(s(i),n-1,br,co)
41 mlabel =
"Test that spline of ALL points would be successful."
43 CALL graf2(real(s),real(f),real(f2),n,
's',
44 & trim(ylabel),trim(mlabel),
" ")
45 IF (
ALLOCATED(xplt))
DEALLOCATE(xplt,yplt)
46 ALLOCATE(xplt(2),yplt(2))
50 yplt(1) = maxval(f2)/6
52 CALL pgline(2, (/.02,.08/), yplt)
53 CALL pgtext (0.1, yplt(1) ,
"Spline Fit")
58 yplt(1) = maxval(f2)/5
59 CALL pgpt1 (xplt(2),yplt(1),2)
60 CALL pgtext (0.1, yplt(1),
"knot location/values")
62 CALL pgpt1 (xplt(2),yplt(1),2)
64 IF (
ALLOCATED(xplt))
DEALLOCATE(xplt,yplt)
66 slctn(1) = 1; slctn(2) =
SIZE(f)
68 ALLOCATE(span(nin/nset1))
70 span = (/(i, i = 1,nin,nset1)/)
71 lastspan = lastspan+nin/nset1-1
73 IF(slctn(i) == 0)
EXIT
75 slctn(i:i+
SIZE(span)-1) = span
78 IF(slctn(i) > nin)slctn(i) = 0
79 IF(slctn(i) == 0) cycle
86 IF((f(i)-f(i-1)) == 0.)cycle
87 c = (f(i+1)-f(i))/(f(i)-f(i-1))
88 IF(c < 0 ) olli(i) = -1
90 IF (sum(olli) /= 0)
THEN
92 IF(olli(i) == 0) cycle
98 IF(slctn(i) == 0) cycle
106 f2(i) = fvn_spline_devaldx(s(i),n-1,br,co)
112 IF((f(i)-f(i-1)) == 0.)cycle
113 c = (f2(i+1)-f2(i))/(f2(i)-f2(i-1))
114 IF(c < 0 ) olli(i) = -1
116 IF (sum(olli) /= 0)
THEN
118 IF(olli(i) == 0) cycle
119 slctn(lastspan+1) = i
120 lastspan = lastspan+1
124 IF(slctn(i) == 0) cycle
128 ALLOCATE(span(
SIZE(f2)))
129 span = (/(i,i = 1,nin)/)
131 CALL ssort( s2, span, nin)
133 2000
FORMAT(
' Original ',5x,
' Sorted ',/,
' Array ',
136 2001
FORMAT(i3,2x,f5.3,5x, i3,2x, f5.3)
137 IF(
ALLOCATED(olli))
DEALLOCATE(olli)
139 ALLOCATE(olli(nin/nset1))
140 olli = (/(i, i = 1,nin,nset1)/)
141 lastspan = lastspan+nin/nset1
143 IF(slctn(i) == 0)
EXIT
145 slctn(i:i+
SIZE(olli)-1) = olli
147 IF(slctn(i) == 0) cycle
151 slctn(lastspan+1:lastspan+nset1) = span(1:nset1)
152 lastspan = lastspan+nset1
154 IF(slctn(i) == 0) cycle
157 IF(
ALLOCATED(olli))
DEALLOCATE(olli)
158 ALLOCATE(olli(lastspan))
159 olli = (/(i,i = 1,lastspan)/)
160 CALL ssort(dble(slctn(1:lastspan)),olli,lastspan)
164 match = (slctn(olli(i)) == slctn(olli(j)) .AND. i /= j)
165 IF(.not. match) cycle
172 IF(slctn(i) == 0) cycle
175 IF(
ALLOCATED(span))
DEALLOCATE(span)
176 ALLOCATE(span(lastspan))
177 span = slctn(olli(lastspan:1:-1))
179 IF(
ALLOCATED(olli))
DEALLOCATE(olli)
180 ALLOCATE(olli(lastspan))
181 olli = (/(i, i = 1,lastspan)/)
182 CALL ssort(dble(span),olli,lastspan)
184 IF(span(olli(i)) == 0)
EXIT
187 IF(
ALLOCATED(slctn))
DEALLOCATE(slctn)
188 ALLOCATE(slctn(lastspan))
189 slctn = span(olli(lastspan:1:-1))
190 ALLOCATE(xplt(lastspan),yplt(lastspan))
191 xplt = s(slctn); yplt = f(slctn)
192 WRITE(xlabel,fmt =
'(i3.3)')lastspan
195 IF (lastspan > 100) slabel =
"Thinning FAILED!!"
196 CALL graf1pt(xplt,yplt,lastspan,
's',trim(ylabel),
197 5
"Akima spline -- knots thinned to "//trim(xlabel),trim(slabel))
201 IF(
ALLOCATED(co))
DEALLOCATE(co,br,s2,f2)
203 ALLOCATE(s2(n),f2(n),br(n),co(4,n))
206 CALL fvn_akima(n,s2,f2,br,co)
208 ALLOCATE(s2(nin),f2(nin))
212 call random_number(s3)
214 s2(i) = s(i)+s3(i)*ds
215 f2(i) = fvn_spline_eval(s(i),n-1,br,co)
220 IF(
ALLOCATED(xplt) )
DEALLOCATE(xplt,yplt)
221 ALLOCATE(xplt(nin),yplt(nin))
223 CALL pgline(nin,xplt,yplt)
227 CALL pgline (nin,real(s),real(f))
230 CALL pgpt1(0.1,real(maxval(f2))/5,2)
231 CALL pgtext (0.18, real(maxval(f2))/5,
"knot location/values")
233 IF(
ALLOCATED(xplt) )
DEALLOCATE(xplt,yplt)
234 ALLOCATE(xplt(2),yplt(2))
235 yplt(1) = real(maxval(f2)/6)
239 CALL pgline(2,xplt,yplt)
240 CALL pgtext (0.18, real(maxval(f2))/6,
"Akima spline result")
241 DEALLOCATE(xplt,yplt)
243 ALLOCATE(s2(
SIZE(slctn)),f2(
SIZE(slctn)))
244 s2 = s(slctn); f2 = f(slctn); n2 =
SIZE(slctn)
245 IF (lastspan .GT. 100) ier = -lastspan
246 END SUBROUTINE doakima