diff options
Diffstat (limited to 'math/deboor/chol1d.f')
| -rw-r--r-- | math/deboor/chol1d.f | 58 |
1 files changed, 58 insertions, 0 deletions
diff --git a/math/deboor/chol1d.f b/math/deboor/chol1d.f new file mode 100644 index 00000000..fa0434ab --- /dev/null +++ b/math/deboor/chol1d.f @@ -0,0 +1,58 @@ + subroutine chol1d ( p, v, qty, npoint, ncol, u, qu ) +c from * a practical guide to splines * by c. de boor +c from * a practical guide to splines * by c. de boor +c to be called in s m o o t h +constructs the upper three diags. in v(i,j), i=2,,npoint-1, j=1,3, of +c the matrix 6*(1-p)*q-transp.*(d**2)*q + p*r, then computes its +c l*l-transp. decomposition and stores it also in v, then applies +c forward and backsubstitution to the right side q-transp.*y in qty +c to obtain the solution in u . + integer ncol,npoint, i,npm1,npm2 + real p,qty(npoint),qu(npoint),u(npoint),v(npoint,7), prev,ratio + * ,six1mp,twop + npm1 = npoint - 1 +c construct 6*(1-p)*q-transp.*(d**2)*q + p*r + six1mp = 6.*(1.-p) + twop = 2.*p + do 2 i=2,npm1 + v(i,1) = six1mp*v(i,5) + twop*(v(i-1,4)+v(i,4)) + v(i,2) = six1mp*v(i,6) + p*v(i,4) + 2 v(i,3) = six1mp*v(i,7) + npm2 = npoint - 2 + if (npm2 .ge. 2) go to 10 + u(1) = 0. + u(2) = qty(2)/v(2,1) + u(3) = 0. + go to 41 +c factorization + 10 do 20 i=2,npm2 + ratio = v(i,2)/v(i,1) + v(i+1,1) = v(i+1,1) - ratio*v(i,2) + v(i+1,2) = v(i+1,2) - ratio*v(i,3) + v(i,2) = ratio + ratio = v(i,3)/v(i,1) + v(i+2,1) = v(i+2,1) - ratio*v(i,3) + 20 v(i,3) = ratio +c +c forward substitution + u(1) = 0. + v(1,3) = 0. + u(2) = qty(2) + do 30 i=2,npm2 + 30 u(i+1) = qty(i+1) - v(i,2)*u(i) - v(i-1,3)*u(i-1) +c back substitution + u(npoint) = 0. + u(npm1) = u(npm1)/v(npm1,1) + i = npm2 + 40 u(i) = u(i)/v(i,1)-u(i+1)*v(i,2)-u(i+2)*v(i,3) + i = i - 1 + if (i .gt. 1) go to 40 +c construct q*u + 41 prev = 0. + do 50 i=2,npoint + qu(i) = (u(i) - u(i-1))/v(i-1,4) + qu(i-1) = qu(i) - prev + 50 prev = qu(i) + qu(npoint) = -qu(npoint) + return + end |