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*8 sq
      integer m,n,maxtr,maxbin,nm,mb,ier
c  ..array arguments..
      real*8 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*8 f,wi,xi
c  ..local array..
      real*8 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