subroutine fpcosp(m,x,y,w,n,t,e,maxtr,maxbin,c,sq,sx,bind,nm,mb,a,
* b,const,z,zz,u,q,info,up,left,right,jbind,ibind,ier)
c ..
c ..scalar arguments..
real sq
integer m,n,maxtr,maxbin,nm,mb,ier
c ..array arguments..
real x(m),y(m),w(m),t(n),e(n),c(n),sx(m),a(n,4),b(nm,maxbin),
* const(n),z(n),zz(n),u(maxbin),q(m,4)
integer info(maxtr),up(maxtr),left(maxtr),right(maxtr),jbind(mb),
* ibind(mb)
logical bind(n)
c ..local scalars..
integer count,i,i1,j,j1,j2,j3,k,kdim,k1,k2,k3,k4,k5,k6,
* l,lp1,l1,l2,l3,merk,nbind,number,n1,n4,n6
real f,wi,xi
c ..local array..
real h(4)
c ..subroutine references..
c fpbspl,fpadno,fpdeno,fpfrno,fpseno
c ..
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c if we use the b-spline representation of s(x) our approximation c
c problem results in a quadratic programming problem: c
c find the b-spline coefficients c(j),j=1,2,...n-4 such that c
c (1) sumi((wi*(yi-sumj(cj*nj(xi))))**2),i=1,2,...m is minimal c
c (2) sumj(cj*n''j(t(l+3)))*e(l) <= 0, l=1,2,...n-6. c
c to solve this problem we use the theil-van de panne procedure. c
c if the inequality constraints (2) are numbered from 1 to n-6, c
c this algorithm finds a subset of constraints ibind(1)..ibind(nbind) c
c such that the solution of the minimization problem (1) with these c
c constraints in equality form, satisfies all constraints. such a c
c feasible solution is optimal if the lagrange parameters associated c
c with that problem with equality constraints, are all positive. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c determine n6, the number of inequality constraints.
n6 = n-6
c fix the parameters which determine these constraints.
do 10 i=1,n6
const(i) = e(i)*(t(i+4)-t(i+1))/(t(i+5)-t(i+2))
10 continue
c initialize the triply linked tree which is used to find the subset
c of constraints ibind(1),...ibind(nbind).
count = 1
info(1) = 0
left(1) = 0
right(1) = 0
up(1) = 1
merk = 1
c set up the normal equations n'nc=n'y where n denotes the m x (n-4)
c observation matrix with elements ni,j = wi*nj(xi) and y is the
c column vector with elements yi*wi.
c from the properties of the b-splines nj(x),j=1,2,...n-4, it follows
c that n'n is a (n-4) x (n-4) positive definit bandmatrix of
c bandwidth 7. the matrices n'n and n'y are built up in a and z.
n4 = n-4
c initialization
do 20 i=1,n4
z(i) = 0.
do 20 j=1,4
a(i,j) = 0.
20 continue
l = 4
lp1 = l+1
do 70 i=1,m
c fetch the current row of the observation matrix.
xi = x(i)
wi = w(i)**2
c search for knot interval t(l) <= xi < t(l+1)
30 if(xi.lt.t(lp1) .or. l.eq.n4) go to 40
l = lp1
lp1 = l+1
go to 30
c evaluate the four non-zero cubic b-splines nj(xi),j=l-3,...l.
40 call fpbspl(t,n,3,xi,l,h)
c store in q these values h(1),h(2),...h(4).
do 50 j=1,4
q(i,j) = h(j)
50 continue
c add the contribution of the current row of the observation matrix
c n to the normal equations.
l3 = l-3
k1 = 0
do 60 j1 = l3,l
k1 = k1+1
f = h(k1)
z(j1) = z(j1)+f*wi*y(i)
k2 = k1
j2 = 4
do 60 j3 = j1,l
a(j3,j2) = a(j3,j2)+f*wi*h(k2)
k2 = k2+1
j2 = j2-1
60 continue
70 continue
c since n'n is a symmetric matrix it can be factorized as
c (3) n'n = (r1)'(d1)(r1)
c with d1 a diagonal matrix and r1 an (n-4) x (n-4) unit upper
c triangular matrix of bandwidth 4. the matrices r1 and d1 are built
c up in a. at the same time we solve the systems of equations
c (4) (r1)'(z2) = n'y
c (5) (d1) (z1) = (z2)
c the vectors z2 and z1 are kept in zz and z.
do 140 i=1,n4
k1 = 1
if(i.lt.4) k1 = 5-i
k2 = i-4+k1
k3 = k2
do 100 j=k1,4
k4 = j-1
k5 = 4-j+k1
f = a(i,j)
if(k1.gt.k4) go to 90
k6 = k2
do 80 k=k1,k4
f = f-a(i,k)*a(k3,k5)*a(k6,4)
k5 = k5+1
k6 = k6+1
80 continue
90 if(j.eq.4) go to 110
a(i,j) = f/a(k3,4)
k3 = k3+1
100 continue
110 a(i,4) = f
f = z(i)
if(i.eq.1) go to 130
k4 = i
do 120 j=k1,3
k = k1+3-j
k4 = k4-1
f = f-a(i,k)*z(k4)*a(k4,4)
120 continue
130 z(i) = f/a(i,4)
zz(i) = f
140 continue
c start computing the least-squares cubic spline without taking account
c of any constraint.
nbind = 0
n1 = 1
ibind(1) = 0
c main loop for the least-squares problems with different subsets of
c the constraints (2) in equality form. the resulting b-spline coeff.
c c and lagrange parameters u are the solution of the system
c ! n'n b' ! ! c ! ! n'y !
c (6) ! ! ! ! = ! !
c ! b 0 ! ! u ! ! 0 !
c z1 is stored into array c.
150 do 160 i=1,n4
c(i) = z(i)
160 continue
c if there are no equality constraints, compute the coeff. c directly.
if(nbind.eq.0) go to 370
c initialization
kdim = n4+nbind
do 170 i=1,nbind
do 170 j=1,kdim
b(j,i) = 0.
170 continue
c matrix b is built up,expressing that the constraints nrs ibind(1),...
c ibind(nbind) must be satisfied in equality form.
do 180 i=1,nbind
l = ibind(i)
b(l,i) = e(l)
b(l+1,i) = -(e(l)+const(l))
b(l+2,i) = const(l)
180 continue
c find the matrix (b1) as the solution of the system of equations
c (7) (r1)'(d1)(b1) = b'
c (b1) is built up in the upper part of the array b(rows 1,...n-4).
do 220 k1=1,nbind
l = ibind(k1)
do 210 i=l,n4
f = b(i,k1)
if(i.eq.1) go to 200
k2 = 3
if(i.lt.4) k2 = i-1
do 190 k3=1,k2
l1 = i-k3
l2 = 4-k3
f = f-b(l1,k1)*a(i,l2)*a(l1,4)
190 continue
200 b(i,k1) = f/a(i,4)
210 continue
220 continue
c factorization of the symmetric matrix -(b1)'(d1)(b1)
c (8) -(b1)'(d1)(b1) = (r2)'(d2)(r2)
c with (d2) a diagonal matrix and (r2) an nbind x nbind unit upper
c triangular matrix. the matrices r2 and d2 are built up in the lower
c part of the array b (rows n-3,n-2,...n-4+nbind).
do 270 i=1,nbind
i1 = i-1
do 260 j=i,nbind
f = 0.
do 230 k=1,n4
f = f+b(k,i)*b(k,j)*a(k,4)
230 continue
k1 = n4+1
if(i1.eq.0) go to 250
do 240 k=1,i1
f = f+b(k1,i)*b(k1,j)*b(k1,k)
k1 = k1+1
240 continue
250 b(k1,j) = -f
if(j.eq.i) go to 260
b(k1,j) = b(k1,j)/b(k1,i)
260 continue
270 continue
c according to (3),(7) and (8) the system of equations (6) becomes
c ! (r1)' 0 ! ! (d1) 0 ! ! (r1) (b1) ! ! c ! ! n'y !
c (9) ! ! ! ! ! ! ! ! = ! !
c ! (b1)' (r2)'! ! 0 (d2) ! ! 0 (r2) ! ! u ! ! 0 !
c backward substitution to obtain the b-spline coefficients c(j),j=1,..
c n-4 and the lagrange parameters u(j),j=1,2,...nbind.
c first step of the backward substitution: solve the system
c ! (r1)'(d1) 0 ! ! (c1) ! ! n'y !
c (10) ! ! ! ! = ! !
c ! (b1)'(d1) (r2)'(d2) ! ! (u1) ! ! 0 !
c from (4) and (5) we know that this is equivalent to
c (11) (c1) = (z1)
c (12) (r2)'(d2)(u1) = -(b1)'(z2)
do 310 i=1,nbind
f = 0.
do 280 j=1,n4
f = f+b(j,i)*zz(j)
280 continue
i1 = i-1
k1 = n4+1
if(i1.eq.0) go to 300
do 290 j=1,i1
f = f+u(j)*b(k1,i)*b(k1,j)
k1 = k1+1
290 continue
300 u(i) = -f/b(k1,i)
310 continue
c second step of the backward substitution: solve the system
c ! (r1) (b1) ! ! c ! ! c1 !
c (13) ! ! ! ! = ! !
c ! 0 (r2) ! ! u ! ! u1 !
k1 = nbind
k2 = kdim
c find the lagrange parameters u.
do 340 i=1,nbind
f = u(k1)
if(i.eq.1) go to 330
k3 = k1+1
do 320 j=k3,nbind
f = f-u(j)*b(k2,j)
320 continue
330 u(k1) = f
k1 = k1-1
k2 = k2-1
340 continue
c find the b-spline coefficients c.
do 360 i=1,n4
f = c(i)
do 350 j=1,nbind
f = f-u(j)*b(i,j)
350 continue
c(i) = f
360 continue
370 k1 = n4
do 390 i=2,n4
k1 = k1-1
f = c(k1)
k2 = 1
if(i.lt.5) k2 = 5-i
k3 = k1
l = 3
do 380 j=k2,3
k3 = k3+1
f = f-a(k3,l)*c(k3)
l = l-1
380 continue
c(k1) = f
390 continue
c test whether the solution of the least-squares problem with the
c constraints ibind(1),...ibind(nbind) in equality form, satisfies
c all of the constraints (2).
k = 1
c number counts the number of violated inequality constraints.
number = 0
do 440 j=1,n6
l = ibind(k)
k = k+1
if(j.eq.l) go to 440
k = k-1
c test whether constraint j is satisfied
f = e(j)*(c(j)-c(j+1))+const(j)*(c(j+2)-c(j+1))
if(f.le.0.) go to 440
c if constraint j is not satisfied, add a branch of length nbind+1
c to the tree. the nodes of this branch contain in their information
c field the number of the constraints ibind(1),...ibind(nbind) and j,
c arranged in increasing order.
number = number+1
k1 = k-1
if(k1.eq.0) go to 410
do 400 i=1,k1
jbind(i) = ibind(i)
400 continue
410 jbind(k) = j
if(l.eq.0) go to 430
do 420 i=k,nbind
jbind(i+1) = ibind(i)
420 continue
430 call fpadno(maxtr,up,left,right,info,count,merk,jbind,n1,ier)
c test whether the storage space which is required for the tree,exceeds
c the available storage space.
if(ier.ne.0) go to 560
440 continue
c test whether the solution of the least-squares problem with equality
c constraints is a feasible solution.
if(number.eq.0) go to 470
c test whether there are still cases with nbind constraints in
c equality form to be considered.
450 if(merk.gt.1) go to 460
nbind = n1
c test whether the number of knots where s''(x)=0 exceeds maxbin.
if(nbind.gt.maxbin) go to 550
n1 = n1+1
ibind(n1) = 0
c search which cases with nbind constraints in equality form
c are going to be considered.
call fpdeno(maxtr,up,left,right,nbind,merk)
c test whether the quadratic programming problem has a solution.
if(merk.eq.1) go to 570
c find a new case with nbind constraints in equality form.
460 call fpseno(maxtr,up,left,right,info,merk,ibind,nbind)
go to 150
c test whether the feasible solution is optimal.
470 ier = 0
do 480 i=1,n6
bind(i) = .false.
480 continue
if(nbind.eq.0) go to 500
do 490 i=1,nbind
if(u(i).le.0.) go to 450
j = ibind(i)
bind(j) = .true.
490 continue
c evaluate s(x) at the data points x(i) and calculate the weighted
c sum of squared residual right hand sides sq.
500 sq = 0.
l = 4
lp1 = 5
do 530 i=1,m
510 if(x(i).lt.t(lp1) .or. l.eq.n4) go to 520
l = lp1
lp1 = l+1
go to 510
520 sx(i) = c(l-3)*q(i,1)+c(l-2)*q(i,2)+c(l-1)*q(i,3)+c(l)*q(i,4)
sq = sq+(w(i)*(y(i)-sx(i)))**2
530 continue
go to 600
c error codes and messages.
550 ier = 1
go to 600
560 ier = 2
go to 600
570 ier = 3
600 return
end