***  For comments, see sbart.r, written by Finbarr O'Sullivan
***  ratfor output of 14 May 1987 version
***  with bug fix in setreg, 27 Sep 1990
      subroutine sbart(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,
     &ispar,lspar,uspar,tol,isetup,xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,
     &abd,p1ip,p2ip,ld4,ldnk,ier)
      real xs(n),ys(n),ws(n),knot(nk+4),coef(nk),sz(n),lev(n),crit,spar,
     &lspar,uspar,tol,xwy(nk),hs0(nk),hs1(nk),hs2(nk),hs3(nk), sg0(nk),
     &sg1(nk),sg2(nk),sg3(nk),abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk)
      integer n,nk,isetup,icrit,ispar,ld4,ldnk,ier
      realt1,t2,ratio, a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, fu,fv,fw,
     &fx,x,ax,bx
      integer i
      if(.not.(isetup.eq.0))goto 23000
      call sgram(sg0,sg1,sg2,sg3,knot,nk)
      call stxwx(xs,ys,ws,n,knot,nk,xwy,hs0,hs1,hs2,hs3)
      t1=0.
      t2=0.
      do 23002 i=3,nk-3
      t1 = t1 + hs0(i)
23002 continue
      do 23004 i=3,nk-3
      t2 = t2 + sg0(i)
23004 continue
      ratio = t1/t2
      isetup = 1
23000 continue
      if(.not.(ispar.eq.1))goto 23006
      call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
     &xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
      return
23006 continue
      ax=lspar
      bx=uspar
      c = 0.5*(3. - sqrt(5.0))
      eps = 1.00
10    eps = eps/2.00
      tol1 = 1.0 + eps
      if(.not.(tol1 .gt. 1.00))goto 23008
      go to 10
23008 continue
      eps = sqrt(eps)
      a = ax
      b = bx
      v = a + c*(b - a)
      w = v
      x = v
      e = 0.0
      spar = x
      call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
     &xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
      fx = crit
      fv = fx
      fw = fx
20    xm = 0.5*(a + b)
      tol1 = eps*abs(x) + tol/3.0
      tol2 = 2.0*tol1
      if(.not.(abs(x - xm) .le. (tol2 - 0.5*(b - a))))goto 23010
      go to 90
23010 continue
      if(.not.(abs(e) .le. tol1))goto 23012
      go to 40
23012 continue
      r = (x - w)*(fx - fv)
      q = (x - v)*(fx - fw)
      p = (x - v)*q - (x - w)*r
      q = 2.00*(q - r)
      if(.not.(q .gt. 0.0))goto 23014
      p = -p
23014 continue
      q = abs(q)
      r = e
      e = d
30    if(.not.(abs(p) .ge. abs(0.5*q*r)))goto 23016
      go to 40
23016 continue
      if(.not.(p .le. q*(a - x)))goto 23018
      go to 40
23018 continue
      if(.not.(p .ge. q*(b - x)))goto 23020
      go to 40
23020 continue
      d = p/q
      u = x + d
      if(.not.((u - a) .lt. tol2))goto 23022
      d = sign(tol1, xm - x)
23022 continue
      if(.not.((b - u) .lt. tol2))goto 23024
      d = sign(tol1, xm - x)
23024 continue
      go to 50
40    if(.not.(x .ge. xm))goto 23026
      e = a - x
23026 continue
      if(.not.(x .lt. xm))goto 23028
      e = b - x
23028 continue
      d = c*e
50    if(.not.(abs(d) .ge. tol1))goto 23030
      u = x + d
23030 continue
      if(.not.(abs(d) .lt. tol1))goto 23032
      u = x + sign(tol1, d)
23032 continue
      spar = u
      call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
     &xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
      fu = crit
      if(.not.(fu .gt. fx))goto 23034
      go to 60
23034 continue
      if(.not.(u .ge. x))goto 23036
      a = x
23036 continue
      if(.not.(u .lt. x))goto 23038
      b = x
23038 continue
      v = w
      fv = fw
      w = x
      fw = fx
      x = u
      fx = fu
      go to 20
60    if(.not.(u .lt. x))goto 23040
      a = u
23040 continue
      if(.not.(u .ge. x))goto 23042
      b = u
23042 continue
      if(.not.(fu .le. fw))goto 23044
      go to 70
23044 continue
      if(.not.(w .eq. x))goto 23046
      go to 70
23046 continue
      if(.not.(fu .le. fv))goto 23048
      go to 80
23048 continue
      if(.not.(v .eq. x))goto 23050
      go to 80
23050 continue
      if(.not.(v .eq. w))goto 23052
      go to 80
23052 continue
      go to 20
70    v = w
      fv = fw
      w = u
      fw = fu
      go to 20
80    v = u
      fv = fu
      go to 20
90    continue
      spar = x
      crit = fx
      return
23007 continue
      return
      end
      subroutine sgram(sg0,sg1,sg2,sg3,tb,nb)
      real sg0(nb),sg1(nb),sg2(nb),sg3(nb),tb(nb+4),vnikx(4,3),work(16),
     &yw1(4),yw2(4),wpt
      integer nb,ileft,ilo,mflag,i,ii,jj
      do 23054 i=1,nb
      sg0(i)=0.
      sg1(i)=0.
      sg2(i)=0.
      sg3(i)=0.
23054 continue
      ilo = 1
      do 23056 i=1,nb
      call interv(tb(1),(nb+1),tb(i),ileft,mflag)
      call bsplvd (tb,4,tb(i),ileft,work,vnikx,3)
      do 23058 ii=1,4
      yw1(ii) = vnikx(ii,3)
23058 continue
      call bsplvd (tb,4,tb(i+1),ileft,work,vnikx,3)
      do 23060 ii=1,4
      yw2(ii) = vnikx(ii,3) - yw1(ii)
23060 continue
      wpt = tb(i+1) - tb(i)
      if(.not.(ileft.ge.4))goto 23062
      do 23064 ii=1,4
      jj=ii
      sg0(ileft-4+ii) = sg0(ileft-4+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
      jj=ii+1
      if(.not.(jj.le.4))goto 23066
      sg1(ileft+ii-4) = sg1(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23066 continue
      jj=ii+2
      if(.not.(jj.le.4))goto 23068
      sg2(ileft+ii-4) = sg2(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23068 continue
      jj=ii+3
      if(.not.(jj.le.4))goto 23070
      sg3(ileft+ii-4) = sg3(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23070 continue
23064 continue
      goto 23063
23062 continue
      if(.not.(ileft.eq.3))goto 23072
      do 23074 ii=1,3
      jj=ii
      sg0(ileft-3+ii) = sg0(ileft-3+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
      jj=ii+1
      if(.not.(jj.le.3))goto 23076
      sg1(ileft+ii-3) = sg1(ileft+ii-3) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23076 continue
      jj=ii+2
      if(.not.(jj.le.3))goto 23078
      sg2(ileft+ii-3) = sg2(ileft+ii-3) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23078 continue
23074 continue
      goto 23073
23072 continue
      if(.not.(ileft.eq.2))goto 23080
      do 23082 ii=1,2
      jj=ii
      sg0(ileft-2+ii) = sg0(ileft-2+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
      jj=ii+1
      if(.not.(jj.le.2))goto 23084
      sg1(ileft+ii-2) = sg1(ileft+ii-2) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23084 continue
23082 continue
      goto 23081
23080 continue
      if(.not.(ileft.eq.1))goto 23086
      do 23088 ii=1,1
      jj=ii
      sg0(ileft-1+ii) = sg0(ileft-1+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
     &ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23088 continue
23086 continue
23081 continue
23073 continue
23063 continue
23056 continue
      return
      end
      subroutine sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,flag)
      realabd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk),wjm3(3),wjm2(2),wjm1(1)
     &,c0,c1,c2,c3
      integerflag,ld4,nk,ldnk,i,j,k
      wjm3(1)=0.
      wjm3(2)=0.
      wjm3(1)=0.
      wjm2(1)=0.
      wjm2(2)=0.
      wjm1(1)=0.
      do 23090 i=1,nk
      j=nk-i+1
      c0 = 1./abd(4,j)
      if(.not.(j.le.nk-3))goto 23092
      c1 = abd(1,j+3)*c0
      c2 = abd(2,j+2)*c0
      c3 = abd(3,j+1)*c0
      goto 23093
23092 continue
      if(.not.(j.eq.nk-2))goto 23094
      c1 = 0.
      c2 = abd(2,j+2)*c0
      c3 = abd(3,j+1)*c0
      goto 23095
23094 continue
      if(.not.(j.eq.nk-1))goto 23096
      c1 = 0.
      c2 = 0.
      c3 = abd(3,j+1)*c0
      goto 23097
23096 continue
      if(.not.(j.eq.nk))goto 23098
      c1 = 0.
      c2 = 0.
      c3 = 0.
23098 continue
23097 continue
23095 continue
23093 continue
      p1ip(1,j) = 0. - (c1*wjm3(1)+c2*wjm3(2)+c3*wjm3(3))
      p1ip(2,j) = 0. - (c1*wjm3(2)+c2*wjm2(1)+c3*wjm2(2))
      p1ip(3,j) = 0. - (c1*wjm3(3)+c2*wjm2(2)+c3*wjm1(1))
      p1ip(4,j) = c0**2 +c1**2*wjm3(1)+2.*c1*c2*wjm3(2)+2.*c1*c3*wjm3(3)
     & +c2**2*wjm2(1)+2.*c2*c3*wjm2(2) +c3**2*wjm1(1)
      wjm3(1)=wjm2(1)
      wjm3(2)=wjm2(2)
      wjm3(3)=p1ip(2,j)
      wjm2(1)=wjm1(1)
      wjm2(2)=p1ip(3,j)
      wjm1(1)=p1ip(4,j)
23090 continue
      if(.not.(flag.eq.0))goto 23100
      return
23100 continue
      do 23102 i=1,nk
      j=nk-i+1
      k=1
23104 if(.not.(k.le.4.and.j+k-1.le.nk))goto 23106
      p2ip(j,j+k-1) = p1ip(5-k,j)
      k=k+1
      goto 23104
23106 continue
23102 continue
      do 23107 i=1,nk
      j=nk-i+1
      k=j-4
23109 if(.not.(k.ge.1))goto 23111
      c0 = 1./abd(4,k)
      c1 = abd(1,k+3)*c0
      c2 = abd(2,k+2)*c0
      c3 = abd(3,k+1)*c0
      p2ip(k,j) = 0. - ( c1*p2ip(k+3,j) +c2*p2ip(k+2,j) +c3*p2ip(k+1,j)
     &)
      k=k-1
      goto 23109
23111 continue
23107 continue
      return
23101 continue
      end
      subroutine sslvrg(x,y,w,n,knot,nk,coef,sz,lev,crit,icrit,spar,
     &ratio,xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,
     &info)
      real x(n),y(n),w(n),knot(nk+4),coef(nk),sz(n),lev(n),crit,ratio,
     &spar,xwy(nk),hs0(nk),hs1(nk),hs2(nk),hs3(nk), sg0(nk),sg1(nk),sg2(
     &nk),sg3(nk),abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk),lambda,b0,b1,
     &b2,b3,eps,vnikx(4,1),work(16),xv,bvalu2,rss,df
      integer n,nk,icrit,ld4,ldnk,i,icoef,ileft,ilo,info,j,mflag
      ilo = 1
      eps = .1e-10
      lambda = ratio*16.**(-2. + spar*(6.))
      do 23112 i=1,nk
      coef(i) = xwy(i)
23112 continue
      do 23114 i=1,nk
      abd(4,i) = hs0(i)+lambda*sg0(i)
23114 continue
      do 23116 i=1,(nk-1)
      abd(3,i+1) = hs1(i)+lambda*sg1(i)
23116 continue
      do 23118 i=1,(nk-2)
      abd(2,i+2) = hs2(i)+lambda*sg2(i)
23118 continue
      do 23120 i=1,(nk-3)
      abd(1,i+3) = hs3(i)+lambda*sg3(i)
23120 continue
      call spbfa(abd,ld4,nk,3,info)
      if(.not.(info.ne.0))goto 23122
      return
23122 continue
      call spbsl(abd,ld4,nk,3,coef)
      icoef = 1
      do 23124 i=1,n
      xv = x(i)
      sz(i) = bvalu2(knot,coef,nk,4,xv,0)
23124 continue
      if(.not.(icrit.eq.0))goto 23126
      return
23126 continue
      call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0)
      do 23128 i=1,n
      xv = x(i)
      call interv(knot(1),(nk+1),xv,ileft,mflag)
      if(.not.(mflag.eq.-1))goto 23130
      ileft = 4
      xv = knot(4)+eps
23130 continue
      if(.not.(mflag.eq.1))goto 23132
      ileft = nk
      xv = knot(nk+1)-eps
23132 continue
      j=ileft-3
      call bsplvd(knot,4,xv,ileft,work,vnikx,1)
      b0=vnikx(1,1)
      b1=vnikx(2,1)
      b2=vnikx(3,1)
      b3=vnikx(4,1)
      lev(i) = (p1ip(4,j)*b0**2 + 2.*p1ip(3,j)*b0*b1 +2.*p1ip(2,j)*b0*
     &b2 + 2.*p1ip(1,j)*b0*b3 +p1ip(4,j+1)*b1**2 + 2.*p1ip(3,j+1)*b1*b2
     &+2.*p1ip(2,j+1)*b1*b3 +p1ip(4,j+2)*b2**2 + 2.*p1ip(3,j+2)*b2*b3 +
     &p1ip(4,j+3)*b3**2 )*w(i)**2
23128 continue
      if(.not.(icrit.eq.1))goto 23134
      rss = 0.
      df = 0.
      do 23136 i=1,n
      rss = rss + ((y(i)-sz(i))*w(i))**2
23136 continue
      do 23138 i=1,n
      df = df + 1.-lev(i)
23138 continue
      crit = (rss/n)/((df/n)**2)
      goto 23135
23134 continue
      crit = 0.
      do 23140 i=1,n
      crit = crit +(((y(i)-sz(i))*w(i))/(1-lev(i)))**2
23140 continue
23135 continue
      return
23127 continue
      end
      subroutine stxwx(x,z,w,k,xknot,n,y,hs0,hs1,hs2,hs3)
      real z(k),w(k),x(k),xknot(n+4),y(n),hs0(n),hs1(n),hs2(n),hs3(n),
     &eps,vnikx(4,1),work(16)
      integer k,n,j,i,ilo,ileft,mflag
      do 23142 i=1,n
      y(i)=0.
      hs0(i)=0.
      hs1(i)=0.
      hs2(i)=0.
      hs3(i)=0.
23142 continue
      ilo=1
      eps = .1e-9
      do 23144 i=1,k
      call interv(xknot(1),(n+1),x(i),ileft,mflag)
      if(.not.(mflag.eq.-1))goto 23146
      write(6,'("Error in hess ",i2)')mflag
      stop
23146 continue
      if(.not.(mflag.eq. 1))goto 23148
      if(.not.(x(i).le.(xknot(ileft)+eps)))goto 23150
      ileft=ileft-1
      goto 23151
23150 continue
      write(6,'("Error in hess ",i2)')mflag
      stop
23151 continue
23148 continue
      call bsplvd (xknot,4,x(i),ileft,work,vnikx,1)
      j= ileft-4+1
      y(j) = y(j)+w(i)**2*z(i)*vnikx(1,1)
      hs0(j)=hs0(j)+w(i)**2*vnikx(1,1)**2
      hs1(j)=hs1(j)+w(i)**2*vnikx(1,1)*vnikx(2,1)
      hs2(j)=hs2(j)+w(i)**2*vnikx(1,1)*vnikx(3,1)
      hs3(j)=hs3(j)+w(i)**2*vnikx(1,1)*vnikx(4,1)
      j= ileft-4+2
      y(j) = y(j)+w(i)**2*z(i)*vnikx(2,1)
      hs0(j)=hs0(j)+w(i)**2*vnikx(2,1)**2
      hs1(j)=hs1(j)+w(i)**2*vnikx(2,1)*vnikx(3,1)
      hs2(j)=hs2(j)+w(i)**2*vnikx(2,1)*vnikx(4,1)
      j= ileft-4+3
      y(j) = y(j)+w(i)**2*z(i)*vnikx(3,1)
      hs0(j)=hs0(j)+w(i)**2*vnikx(3,1)**2
      hs1(j)=hs1(j)+w(i)**2*vnikx(3,1)*vnikx(4,1)
      j= ileft-4+4
      y(j) = y(j)+w(i)**2*z(i)*vnikx(4,1)
      hs0(j)=hs0(j)+w(i)**2*vnikx(4,1)**2
23144 continue
      return
      end
      subroutine sknotl(x,n,knot,k)
      real x(n),knot(n+6),a1,a2,a3,a4
      integer n,k,ndk,j
      a1 = log(50.)/log(2.)
      a2 = log(100.)/log(2.)
      a3 = log(140.)/log(2.)
      a4 = log(200.)/log(2.)
      if(.not.(n.lt.50))goto 23152
      ndk = n
      goto 23153
23152 continue
      if(.not.(n.ge.50 .and. n.lt.200))goto 23154
      ndk = 2.**(a1+(a2-a1)*(n-50.)/150.)
      goto 23155
23154 continue
      if(.not.(n.ge.200 .and. n.lt.800))goto 23156
      ndk = 2.**(a2+(a3-a2)*(n-200.)/600.)
      goto 23157
23156 continue
      if(.not.(n.ge.800 .and. n.lt.3200))goto 23158
      ndk = 2.**(a3+(a4-a3)*(n-800.)/2400.)
      goto 23159
23158 continue
      if(.not.(n.ge.3200))goto 23160
      ndk = 200. + (n-3200)**.2
23160 continue
23159 continue
23157 continue
23155 continue
23153 continue
      k = ndk + 6
      do 23162 j=1,3
      knot(j) = x(1)
23162 continue
      do 23164 j=1,ndk
      knot(j+3) = x( 1 + (j-1)*(n-1)/(ndk-1) )
23164 continue
      do 23166 j=1,3
      knot(ndk+3+j) = x(n)
23166 continue
      return
      end
      subroutine setreg(x,y,w,n,xw,nx,min,range,knot,nk)
      real x(n),y(n),w(n),xw(n),min,range,knot(n+6)
      integer n,nk,nx
      real eps
      integer ycnt,i,k
      call scopy(n,x,1,xw,1)
      call sortpr(x,n,w)
      call sortpr(xw,n,y)
      range = x(n)-x(1)
      min = x(1)
      eps = .1e-9
      do 23168 i=1,n
      x(i) = (x(i)-min)/range
23168 continue
      call scopy(n,x,1,xw,1)
      nx = 1
      x(nx) = x(1)
      w(nx) = w(1)
      y(nx) = y(1)*w(1)
      i=2
23170 if(.not.(i.le.n))goto 23172
      if(.not.(xw(i).gt.x(nx)+eps))goto 23173
      if(.not.(w(nx).gt.0.0))goto 23175
      y(nx) = y(nx)/w(nx)
23175 continue
      nx = nx + 1
      x(nx) = x(i)
      y(nx) = y(i)*w(i)
      w(nx) = w(i)
      goto 23174
23173 continue
      y(nx) = y(nx)+y(i)*w(i)
      w(nx) = w(nx) + w(i)
23174 continue
      i=i+1
      goto 23170
23172 continue
      if(.not.(w(nx).gt.0.0))goto 23177
      y(nx) = y(nx)/w(nx)
23177 continue
      i=1
23179 if(.not.(i.le.nx))goto 23181
      if(.not.(w(i).gt.0))goto 23182
      w(i)=sqrt(w(i))
23182 continue
      i=i+1
      goto 23179
23181 continue
      call sknotl(x,nx,knot,k)
      nk = k-4
      return
      end