diff options
Diffstat (limited to 'math/deboor/setupq.f')
| -rw-r--r-- | math/deboor/setupq.f | 40 |
1 files changed, 40 insertions, 0 deletions
diff --git a/math/deboor/setupq.f b/math/deboor/setupq.f new file mode 100644 index 00000000..dea21a3b --- /dev/null +++ b/math/deboor/setupq.f @@ -0,0 +1,40 @@ + subroutine setupq ( x, dx, y, npoint, v, qty ) +c from * a practical guide to splines * by c. de boor +c to be called in s m o o t h +c put delx = x(.+1) - x(.) into v(.,4), +c put the three bands of q-transp*d into v(.,1-3), and +c put the three bands of (d*q)-transp*(d*q) at and above the diagonal +c into v(.,5-7) . +c here, q is the tridiagonal matrix of order (npoint-2,npoint) +c with general row 1/delx(i) , -1/delx(i) - 1/delx(i+1) , 1/delx(i+1) +c and d is the diagonal matrix with general row dx(i) . + integer npoint, i,npm1 + real dx(npoint),qty(npoint),v(npoint,7),x(npoint),y(npoint), + * diff, prev + npm1 = npoint - 1 + v(1,4) = x(2) - x(1) + do 11 i=2,npm1 + v(i,4) = x(i+1) - x(i) + v(i,1) = dx(i-1)/v(i-1,4) + v(i,2) = - dx(i)/v(i,4) - dx(i)/v(i-1,4) + 11 v(i,3) = dx(i+1)/v(i,4) + v(npoint,1) = 0. + do 12 i=2,npm1 + 12 v(i,5) = v(i,1)**2 + v(i,2)**2 + v(i,3)**2 + if (npm1 .lt. 3) go to 14 + do 13 i=3,npm1 + 13 v(i-1,6) = v(i-1,2)*v(i,1) + v(i-1,3)*v(i,2) + 14 v(npm1,6) = 0. + if (npm1 .lt. 4) go to 16 + do 15 i=4,npm1 + 15 v(i-2,7) = v(i-2,3)*v(i,1) + 16 v(npm1-1,7) = 0. + v(npm1,7) = 0. +construct q-transp. * y in qty. + prev = (y(2) - y(1))/v(1,4) + do 21 i=2,npm1 + diff = (y(i+1)-y(i))/v(i,4) + qty(i) = diff - prev + 21 prev = diff + return + end |