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