src/library/stats/src/loessf.f

C
C     Copyright (C) 1998--2020  The R Core Team

C  The authors of this software are Cleveland, Grosse, and Shyu.
C  Copyright (c) 1989, 1992 by AT&T.
C  Permission to use, copy, modify, and distribute this software for any
C  purpose without fee is hereby granted, provided that this entire notice
C  is included in all copies of any software which is or includes a copy
C  or modification of this software and in all copies of the supporting
C  documentation for such software.
C  THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
C  WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
C  REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
C  OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

C       altered by B.D. Ripley to
C
C       remove unused variables
C       make phi in ehg139 double precision to match calling sequence
C       pass integer not logical from C
C
C       by M. Maechler, renaming ehg182() to loesswarn()

C     Note that  loesswarn(errormsg_code)  is in ./loessc.c

      subroutine ehg126(d,n,vc,x,v,nvmax)
      integer d,execnt,i,j,k,n,nvmax,vc
      DOUBLE PRECISION machin,alpha,beta,mu,t
      DOUBLE PRECISION v(nvmax,d), x(n,d)

      DOUBLE PRECISION D1MACH
      external D1MACH
      save machin,execnt
      data execnt /0/
c     MachInf -> machin
      execnt=execnt+1
      if(execnt.eq.1)then
c     initialize  d1mach(2) === DBL_MAX:
         machin=D1MACH(2)
      end if
c     fill in vertices for bounding box of $x$
c     lower left, upper right
      do 3 k=1,d
         alpha=machin
         beta=-machin
         do 4 i=1,n
            t=x(i,k)
            alpha=min(alpha,t)
            beta=max(beta,t)
    4    continue
c        expand the box a little
         mu=0.005D0*max(beta-alpha,1.d-10*max(DABS(alpha),DABS(beta))+
     +        1.d-30)
         alpha=alpha-mu
         beta=beta+mu
         v(1,k)=alpha
         v(vc,k)=beta
    3 continue
c     remaining vertices
      do 5 i=2,vc-1
         j=i-1
         do 6 k=1,d
            v(i,k)=v(1+mod(j,2)*(vc-1),k)
c    Integer division would do here
            j=INT(DBLE(j)/2.D0)
    6    continue
    5 continue
      return
      end

      subroutine ehg125(p,nv,v,vhit,nvmax,d,k,t,r,s,f,l,u)
      logical i1,i2,match
      integer d,h,i,i3,j,k,m,mm,nv,nvmax,p,r,s
      integer f(r,0:1,s),l(r,0:1,s),u(r,0:1,s),vhit(nvmax)
      DOUBLE PRECISION t
      DOUBLE PRECISION v(nvmax,d)
      external loesswarn
      h=nv
      do 3 i=1,r
         do 4 j=1,s
            h=h+1
            do 5 i3=1,d
               v(h,i3)=v(f(i,0,j),i3)
    5       continue
            v(h,k)=t
c           check for redundant vertex
            match=.false.
            m=1
c           top of while loop
    6       if(.not.match)then
               i1=(m.le.nv)
            else
               i1=.false.
            end if
            if(.not.(i1))goto 7
               match=(v(m,1).eq.v(h,1))
               mm=2
c              top of while loop
    8          if(match)then
                  i2=(mm.le.d)
               else
                  i2=.false.
               end if
               if(.not.(i2))goto 9
                  match=(v(m,mm).eq.v(h,mm))
                  mm=mm+1
                  goto 8
c              bottom of while loop
    9          m=m+1
               goto 6
c           bottom of while loop
    7       m=m-1
            if(match)then
               h=h-1
            else
               m=h
               if(vhit(1).ge.0)then
                  vhit(m)=p
               end if
            end if
            l(i,0,j)=f(i,0,j)
            l(i,1,j)=m
            u(i,0,j)=m
            u(i,1,j)=f(i,1,j)
    4    continue
    3 continue
      nv=h
      if(.not.(nv.le.nvmax))then
         call loesswarn(180)
      end if
      return
      end

      integer function ehg138(i,z,a,xi,lo,hi,ncmax)
      logical i1
      integer i,j,ncmax
      integer a(ncmax),hi(ncmax),lo(ncmax)
      DOUBLE PRECISION xi(ncmax),z(8)
c     descend tree until leaf or ambiguous
      j=i
c     top of while loop
    3 if(a(j).ne.0)then
         i1=(z(a(j)).ne.xi(j))
      else
         i1=.false.
      end if
      if(.not.(i1))goto 4
         if(z(a(j)).le.xi(j))then
            j=lo(j)
         else
            j=hi(j)
         end if
         goto 3
c     bottom of while loop
    4 ehg138=j
      return
      end

      subroutine ehg106(il,ir,k,nk,p,pi,n)

c Partial sorting of p(1, il:ir) returning the sort indices pi() only
c such that p(1, pi(k)) is correct

c     implicit none
c Arguments
c  Input:
      integer il,ir,k,nk,n
      DOUBLE PRECISION p(nk,n)
c     using only       p(1, pi(*))
c  Output:
      integer pi(n)

c Variables
      DOUBLE PRECISION t
      integer i,ii,j,l,r

c     find the $k$-th smallest of $n$ elements
c     Floyd+Rivest, CACM Mar '75, Algorithm 489
      l=il
      r=ir
c     while (l < r )
    3 if(.not.(l.lt.r))goto 4
c        to avoid recursion, sophisticated partition deleted
c        partition $x sub {l..r}$ about $t$
         t=p(1,pi(k))
         i=l
         j=r
         ii=pi(l)
         pi(l)=pi(k)
         pi(k)=ii
         if(t.lt.p(1,pi(r)))then
            ii=pi(l)
            pi(l)=pi(r)
            pi(r)=ii
         end if
c        top of while loop
    5    if(.not.(i.lt.j))goto 6
            ii=pi(i)
            pi(i)=pi(j)
            pi(j)=ii
            i=i+1
            j=j-1
c           top of while loop
    7       if(.not.(p(1,pi(i)).lt.t))goto 8
               i=i+1
               goto 7
c           bottom of while loop
    8       continue
c           top of while loop
    9       if(.not.(t.lt.p(1,pi(j))))goto 10
               j=j-1
               goto 9
c           bottom of while loop
   10       goto 5
c        bottom of while loop
    6    if(p(1,pi(l)).eq.t)then
            ii=pi(l)
            pi(l)=pi(j)
            pi(j)=ii
         else
            j=j+1
            ii=pi(r)
            pi(r)=pi(j)
            pi(j)=ii
         end if
         if(j.le.k)then
            l=j+1
         end if
         if(k.le.j)then
            r=j-1
         end if
         goto 3
c     bottom of while loop
    4 return
      end


      subroutine ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,
     +     rcond,sing,sigma,u,e,dgamma,qraux,work,tol,dd,tdeg,cdeg,s)
      integer column,d,dd,execnt,i,i3,i9,info,inorm2,j,jj,jpvt,k,kernel,
     +     n,nf,od,sing,tdeg
      integer cdeg(8),psi(n)
      double precision machep,f,i1,i10,i2,i4,i5,i6,i7,i8,rcond,rho,scal,
     +     tol
      double precision g(15),sigma(15),u(15,15),e(15,15),b(nf,k),
     +     colnor(15),dist(n),eta(nf),dgamma(15),q(d),qraux(15),rw(n),
     +     s(0:od),w(nf),work(15),x(n,d),y(n)

      integer idamax
      double precision d1mach, ddot

      external ehg106,loesswarn,ehg184,dqrdc,dqrsl,dsvdc
      external idamax, d1mach, ddot

      save machep,execnt
      data execnt /0/
c     colnorm -> colnor
c     E -> g
c     MachEps -> machep
c     V -> e
c     X -> b
      execnt=execnt+1
      if(execnt.eq.1)then
c     initialize  d1mach(4) === 1 / DBL_EPSILON === 2^52  :
         machep=d1mach(4)
      end if
c     sort by distance
      do 3 i3=1,n
         dist(i3)=0
    3 continue
      do 4 j=1,dd
         i4=q(j)
         do 5 i3=1,n
            dist(i3)=dist(i3)+(x(i3,j)-i4)**2
    5    continue
    4 continue
      call ehg106(1,n,nf,1,dist,psi,n)
      rho=dist(psi(nf))*max(1.d0,f)
      if(rho .le. 0)then
         call loesswarn(120)
      end if
c     compute neighborhood weights
      if(kernel.eq.2)then
         do 6 i=1,nf
            if(dist(psi(i)).lt.rho)then
               i1=dsqrt(rw(psi(i)))
            else
               i1=0
            end if
            w(i)=i1
    6    continue
      else
         do 7 i3=1,nf
            w(i3)=dsqrt(dist(psi(i3))/rho)
    7    continue
         do 8 i3=1,nf
            w(i3)=dsqrt(rw(psi(i3))*(1-w(i3)**3)**3)
    8    continue
      end if
      if(dabs(w(idamax(nf,w,1))).eq.0)then
         call ehg184('at ',q(1),dd,1)
         call ehg184('radius ',rho,1,1)
         if(.not..false.)then
            call loesswarn(121)
         end if
      end if
c     fill design matrix
      column=1
      do 9 i3=1,nf
         b(i3,column)=w(i3)
    9 continue
      if(tdeg.ge.1)then
         do 10 j=1,d
            if(cdeg(j).ge.1)then
               column=column+1
               i5=q(j)
               do 11 i3=1,nf
                  b(i3,column)=w(i3)*(x(psi(i3),j)-i5)
   11          continue
            end if
   10    continue
      end if
      if(tdeg.ge.2)then
         do 12 j=1,d
            if(cdeg(j).ge.1)then
               if(cdeg(j).ge.2)then
                  column=column+1
                  i6=q(j)
                  do 13 i3=1,nf
                     b(i3,column)=w(i3)*(x(psi(i3),j)-i6)**2
   13             continue
               end if
               do 14 jj=j+1,d
                  if(cdeg(jj).ge.1)then
                     column=column+1
                     i7=q(j)
                     i8=q(jj)
                     do 15 i3=1,nf
                        b(i3,column)=w(i3)*(x(psi(i3),j)-i7)*(x(psi(i3),
     +jj)-i8)
   15                continue
                  end if
   14          continue
            end if
   12    continue
         k=column
      end if
      do 16 i3=1,nf
         eta(i3)=w(i3)*y(psi(i3))
   16 continue
c     equilibrate columns
      do 17 j=1,k
         scal=0
         do 18 inorm2=1,nf
            scal=scal+b(inorm2,j)**2
   18    continue
         scal=dsqrt(scal)
         if(0.lt.scal)then
            do 19 i3=1,nf
               b(i3,j)=b(i3,j)/scal
   19       continue
            colnor(j)=scal
         else
            colnor(j)=1
         end if
   17 continue
c     singular value decomposition
      call dqrdc(b,nf,nf,k,qraux,jpvt,work,0)
      call dqrsl(b,nf,nf,k,qraux,eta,work,eta,eta,work,work,1000,info)
      do 20 i9=1,k
         do 21 i3=1,k
            u(i3,i9)=0
   21    continue
   20 continue
      do 22 i=1,k
         do 23 j=i,k
c FIXME: this has i = 3 vs bound 2 in a ggplot2 test
            u(i,j)=b(i,j)
   23    continue
   22 continue
      call dsvdc(u,15,k,k,sigma,g,u,15,e,15,work,21,info)
      if(.not.(info.eq.0))then
         call loesswarn(182)
      end if
      tol=sigma(1)*(100*machep)
      rcond=min(rcond,sigma(k)/sigma(1))
      if(sigma(k).le.tol)then
         sing=sing+1
         if(sing.eq.1)then
            call ehg184('pseudoinverse used at',q(1),d,1)
            call ehg184('neighborhood radius',dsqrt(rho),1,1)
            call ehg184('reciprocal condition number ',rcond,1,1)
         else
            if(sing.eq.2)then
               call ehg184('There are other near singularities as well.'
     +,rho,1,1)
            end if
         end if
      end if
c     compensate for equilibration
      do 24 j=1,k
         i10=colnor(j)
         do 25 i3=1,k
            e(j,i3)=e(j,i3)/i10
   25    continue
   24 continue
c     solve least squares problem
      do 26 j=1,k
         if(tol.lt.sigma(j))then
            i2=ddot(k,u(1,j),1,eta,1)/sigma(j)
         else
            i2=0.d0
         end if
         dgamma(j)=i2
   26 continue
      do 27 j=0,od
c        bug fix 2006-07-04 for k=1, od>1.   (thanks btyner@gmail.com)
         if(j.lt.k)then
            s(j)=ddot(k,e(j+1,1),15,dgamma,1)
         else
            s(j)=0.d0
         end if
   27 continue
      return
      end

      subroutine ehg131(x,y,rw,trl,diagl,kernel,k,n,d,nc,ncmax,vc,nv,
     +     nvmax,nf,f,a,c,hi,lo,pi,psi,v,vhit,vval,xi,dist,eta,b,ntol,
     +     fd,w,vval2,rcond,sing,dd,tdeg,cdeg,lq,lf,setlf)
      logical setlf
      integer identi,d,dd,i1,i2,j,k,kernel,n,nc,ncmax,nf,ntol,nv,
     +     nvmax,sing,tdeg,vc
      integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax),
     +     lo(ncmax),pi(n),psi(n),vhit(nvmax)
      double precision f,fd,rcond,trl
      double precision lf(0:d,nvmax,nf),b(*),delta(8),diagl(n),dist(n),
     +     eta(nf),rw(n),v(nvmax,d),vval(0:d,nvmax),vval2(0:d,nvmax),
     +     w(nf),x(n,d),xi(ncmax),y(n)
      external ehg126,loesswarn,ehg139,ehg124
      double precision dnrm2
      external dnrm2
c     Identity -> identi
c     X -> b
      if(.not.(d.le.8))then
         call loesswarn(101)
      end if
c     build $k$-d tree
      call ehg126(d,n,vc,x,v,nvmax)
      nv=vc
      nc=1
      do 3 j=1,vc
         c(j,nc)=j
         vhit(j)=0
    3 continue
      do 4 i1=1,d
         delta(i1)=v(vc,i1)-v(1,i1)
    4 continue
      fd=fd*dnrm2(d,delta,1)
      do 5 identi=1,n
         pi(identi)=identi
    5 continue
      call ehg124(1,n,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v,vhit,nvmax,
     +ntol,fd,dd)
c     smooth
      if(trl.ne.0)then
         do 6 i2=1,nv
            do 7 i1=0,d
               vval2(i1,i2)=0
    7       continue
    6    continue
      end if
      call ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k,dist,
     +     dist,eta,b,d,w,diagl,vval2,nc,vc,a,xi,lo,hi,c,rcond,
     +     sing,dd,tdeg,cdeg,lq,lf,setlf,vval)
      return
      end

c called from  lowese()  only :
      subroutine ehg133(d,vc,nvmax,ncmax, a,c,hi,lo, v,vval,xi,m,z,s)
      integer           d,vc,nvmax,ncmax,                      m
      integer           a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax)
      double precision v(nvmax,d),vval(0:d,nvmax),xi(ncmax),z(m,d),s(m)
c Var
      double precision delta(8)
      integer i,i1

      double precision ehg128
      external ehg128

      do 3 i=1,m
         do 4 i1=1,d
            delta(i1)=z(i,i1)
    4    continue
         s(i)=ehg128(delta,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval)
    3 continue
      return
      end

      subroutine ehg140(iw,i,j)
      integer i,j
      integer iw(i)
      iw(i)=j
      return
      end

      subroutine ehg141(trl,n,deg,k,d,nsing,dk,delta1,delta2)
      double precision trl,delta1,delta2
      integer n,deg,k,d,nsing,dk

      double precision c(48), c1, c2, c3, c4,  corx,z,zz(1)
      integer i

      external ehg176
      double precision ehg176

c     coef, d, deg, del
      data c / .2971620d0,.3802660d0,.5886043d0,.4263766d0,.3346498d0,
     +.6271053d0,.5241198d0,.3484836d0,.6687687d0,.6338795d0,.4076457d0,
     +.7207693d0,.1611761d0,.3091323d0,.4401023d0,.2939609d0,.3580278d0,
     +.5555741d0,.3972390d0,.4171278d0,.6293196d0,.4675173d0,.4699070d0,
     +.6674802d0,.2848308d0,.2254512d0,.2914126d0,.5393624d0,.2517230d0,
     +.3898970d0,.7603231d0,.2969113d0,.4740130d0,.9664956d0,.3629838d0,
     +.5348889d0,.2075670d0,.2822574d0,.2369957d0,.3911566d0,.2981154d0,
     +.3623232d0,.5508869d0,.3501989d0,.4371032d0,.7002667d0,.4291632d0,
     +.4930370d0 /

      if(deg.eq.0) dk=1
      if(deg.eq.1) dk=d+1
      if(deg.eq.2) dk=int(dble((d+2)*(d+1))/2.d0)
      corx=dsqrt(k/dble(n))
      z=(dsqrt(k/trl)-corx)/(1-corx)
      if(nsing .eq. 0 .and. 1 .lt. z) then
         call ehg184('Chernobyl! trL<k',trl,1,1)
      else if(z .lt. 0) then
         call ehg184('Chernobyl! trL>n',trl,1,1)
      endif
      z=min(1.0d0,max(0.0d0,z))
c R fix
      zz(1)=z
      c4=dexp(ehg176(zz))
      i=1+3*(min(d,4)-1+4*(deg-1))
      if(d.le.4)then
         c1=c(i)
         c2=c(i+1)
         c3=c(i+2)
      else
         c1=c(i)+(d-4)*(c(i)-c(i-3))
         c2=c(i+1)+(d-4)*(c(i+1)-c(i-2))
         c3=c(i+2)+(d-4)*(c(i+2)-c(i-1))
      endif
      delta1=n-trl*dexp(c1*z**c2*(1-z)**c3*c4)
      i=i+24
      if(d.le.4)then
         c1=c(i)
         c2=c(i+1)
         c3=c(i+2)
      else
         c1=c(i)+(d-4)*(c(i)-c(i-3))
         c2=c(i+1)+(d-4)*(c(i+1)-c(i-2))
         c3=c(i+2)+(d-4)*(c(i+2)-c(i-1))
      endif
      delta2=n-trl*dexp(c1*z**c2*(1-z)**c3*c4)
      return
      end

      subroutine lowesc(n,l,ll,trl,delta1,delta2)
      integer i,j,n
      double precision delta1,delta2,trl
      double precision l(n,n),ll(n,n)
      double precision ddot
      external ddot
c     compute $LL~=~(I-L)(I-L)'$
      do 3 i=1,n
         l(i,i)=l(i,i)-1
    3 continue
      do 4 i=1,n
         do 5 j=1,i
            ll(i,j)=ddot(n,l(i,1),n,l(j,1),n)
    5    continue
    4 continue
      do 6 i=1,n
         do 7 j=i+1,n
            ll(i,j)=ll(j,i)
    7    continue
    6 continue
      do 8 i=1,n
         l(i,i)=l(i,i)+1
    8 continue
c     accumulate first two traces
      trl=0
      delta1=0
      do 9 i=1,n
         trl=trl+l(i,i)
         delta1=delta1+ll(i,i)
    9 continue
c     $delta sub 2 = "tr" LL sup 2$
      delta2=0
      do 10 i=1,n
         delta2=delta2+ddot(n,ll(i,1),n,ll(1,i),1)
   10 continue
      return
      end

      subroutine ehg169(d,vc,nc,ncmax,nv,nvmax,v,a,xi,c,hi,lo)
      integer           d,vc,nc,ncmax,nv,nvmax
      integer           a(ncmax), c(vc,ncmax), hi(ncmax), lo(ncmax)
      DOUBLE PRECISION v(nvmax,d),xi(ncmax)

      integer novhit(1),i,j,k,mc,mv,p
      external ehg125,loesswarn
      integer ifloor
      external ifloor

c     as in bbox
c     remaining vertices
      do 3 i=2,vc-1
         j=i-1
         do 4 k=1,d
            v(i,k)=v(1+mod(j,2)*(vc-1),k)
            j=ifloor(DBLE(j)/2.D0)
    4    continue
    3 continue
c     as in ehg131
      mc=1
      mv=vc
      novhit(1)=-1
      do 5 j=1,vc
         c(j,mc)=j
    5 continue
c     as in rbuild
      p=1
c     top of while loop
    6 if(.not.(p.le.nc))goto 7
         if(a(p).ne.0)then
            k=a(p)
c           left son
            mc=mc+1
            lo(p)=mc
c           right son
            mc=mc+1
            hi(p)=mc
            call ehg125(p,mv,v,novhit,nvmax,d,k,xi(p),2**(k-1),2**(d-k),
     +           c(1,p),c(1,lo(p)),c(1,hi(p)))
         end if
         p=p+1
         goto 6
c     bottom of while loop
    7 if(.not.(mc.eq.nc))then
         call loesswarn(193)
      end if
      if(.not.(mv.eq.nv))then
         call loesswarn(193)
      end if
      return
      end

      DOUBLE PRECISION function ehg176(z)
c
      DOUBLE PRECISION z(*)
c
      integer d,vc,nv,nc
      integer a(17), c(2,17)
      integer hi(17), lo(17)
      DOUBLE PRECISION v(10,1)
      DOUBLE PRECISION vval(0:1,10)
      DOUBLE PRECISION xi(17)
      double precision ehg128
      external ehg128

      data d,vc,nv,nc /1,2,10,17/
      data a(1) /1/
      data hi(1),lo(1),xi(1) /3,2,0.3705D0/
      data c(1,1) /1/
      data c(2,1) /2/
      data a(2) /1/
      data hi(2),lo(2),xi(2) /5,4,0.2017D0/
      data c(1,2) /1/
      data c(2,2) /3/
      data a(3) /1/
      data hi(3),lo(3),xi(3) /7,6,0.5591D0/
      data c(1,3) /3/
      data c(2,3) /2/
      data a(4) /1/
      data hi(4),lo(4),xi(4) /9,8,0.1204D0/
      data c(1,4) /1/
      data c(2,4) /4/
      data a(5) /1/
      data hi(5),lo(5),xi(5) /11,10,0.2815D0/
      data c(1,5) /4/
      data c(2,5) /3/
      data a(6) /1/
      data hi(6),lo(6),xi(6) /13,12,0.4536D0/
      data c(1,6) /3/
      data c(2,6) /5/
      data a(7) /1/
      data hi(7),lo(7),xi(7) /15,14,0.7132D0/
      data c(1,7) /5/
      data c(2,7) /2/
      data a(8) /0/
      data c(1,8) /1/
      data c(2,8) /6/
      data a(9) /0/
      data c(1,9) /6/
      data c(2,9) /4/
      data a(10) /0/
      data c(1,10) /4/
      data c(2,10) /7/
      data a(11) /0/
      data c(1,11) /7/
      data c(2,11) /3/
      data a(12) /0/
      data c(1,12) /3/
      data c(2,12) /8/
      data a(13) /0/
      data c(1,13) /8/
      data c(2,13) /5/
      data a(14) /0/
      data c(1,14) /5/
      data c(2,14) /9/
      data a(15) /1/
      data hi(15),lo(15),xi(15) /17,16,0.8751D0/
      data c(1,15) /9/
      data c(2,15) /2/
      data a(16) /0/
      data c(1,16) /9/
      data c(2,16) /10/
      data a(17) /0/
      data c(1,17) /10/
      data c(2,17) /2/
      data vval(0,1) /-9.0572D-2/
      data v(1,1) /-5.D-3/
      data vval(1,1) /4.4844D0/
      data vval(0,2) /-1.0856D-2/
      data v(2,1) /1.005D0/
      data vval(1,2) /-0.7736D0/
      data vval(0,3) /-5.3718D-2/
      data v(3,1) /0.3705D0/
      data vval(1,3) /-0.3495D0/
      data vval(0,4) /2.6152D-2/
      data v(4,1) /0.2017D0/
      data vval(1,4) /-0.7286D0/
      data vval(0,5) /-5.8387D-2/
      data v(5,1) /0.5591D0/
      data vval(1,5) /0.1611D0/
      data vval(0,6) /9.5807D-2/
      data v(6,1) /0.1204D0/
      data vval(1,6) /-0.7978D0/
      data vval(0,7) /-3.1926D-2/
      data v(7,1) /0.2815D0/
      data vval(1,7) /-0.4457D0/
      data vval(0,8) /-6.4170D-2/
      data v(8,1) /0.4536D0/
      data vval(1,8) /3.2813D-2/
      data vval(0,9) /-2.0636D-2/
      data v(9,1) /0.7132D0/
      data vval(1,9) /0.3350D0/
      data vval(0,10) /4.0172D-2/
      data v(10,1) /0.8751D0/
      data vval(1,10) /-4.1032D-2/
      ehg176=ehg128(z,d,nc,vc,a,xi,lo,hi,c,v,nv,vval)
      end

      subroutine lowesa(trl,n,d,tau,nsing,delta1,delta2)
      integer               n,d,tau,nsing
      double precision  trl, delta1,delta2

      integer dka,dkb
      double precision alpha,d1a,d1b,d2a,d2b
      external ehg141

      call ehg141(trl,n,1,tau,d,nsing,dka,d1a,d2a)
      call ehg141(trl,n,2,tau,d,nsing,dkb,d1b,d2b)
      alpha=dble(tau-dka)/dble(dkb-dka)
      delta1=(1-alpha)*d1a+alpha*d1b
      delta2=(1-alpha)*d2a+alpha*d2b
      return
      end

      subroutine ehg191(m,z,l,d,n,nf,nv,ncmax,vc,a,xi,lo,hi,c,v,nvmax,
     +                  vval2,lf,lq)
c Args
      integer m,d,n,nf,nv,ncmax,nvmax,vc
      double precision z(m,d), l(m,n), xi(ncmax), v(nvmax,d),
     +     vval2(0:d,nvmax), lf(0:d,nvmax,nf)
      integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),lo(ncmax),hi(ncmax)
c Var
      integer lq1,i,i1,i2,j,p
      double precision zi(8)
      double precision ehg128
      external ehg128

      do 3 j=1,n
         do 4 i2=1,nv
            do 5 i1=0,d
               vval2(i1,i2)=0
    5       continue
    4    continue
         do 6 i=1,nv
c           linear search for i in Lq
            lq1=lq(i,1)
            lq(i,1)=j
            p=nf
c           top of while loop
    7       if(.not.(lq(i,p).ne.j))goto 8
               p=p-1
               goto 7
c           bottom of while loop
    8       lq(i,1)=lq1
            if(lq(i,p).eq.j)then
               do 9 i1=0,d
                  vval2(i1,i)=lf(i1,i,p)
    9          continue
            end if
    6    continue
         do 10 i=1,m
            do 11 i1=1,d
               zi(i1)=z(i,i1)
   11       continue
            l(i,j)=ehg128(zi,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,vval2)
   10    continue
    3 continue
      return
      end

      subroutine ehg196(tau,d,f,trl)
      integer d,dka,dkb,tau
      double precision alpha,f,trl,trla,trlb
      external ehg197
      call ehg197(1,d,f,dka,trla)
      call ehg197(2,d,f,dkb,trlb)
      alpha=dble(tau-dka)/dble(dkb-dka)
      trl=(1-alpha)*trla+alpha*trlb
      return
      end

      subroutine ehg197(deg,d,f,dk,trl)
      integer deg,d,dk
      double precision f, trl

      double precision g1
      dk = 0
      if(deg.eq.1) dk=d+1
      if(deg.eq.2) dk=int(dble((d+2)*(d+1))/2.d0)
      g1 = (-0.08125d0*d+0.13d0)*d+1.05d0
      trl = dk*(1+max(0.d0,(g1-f)/f))
      return
      end

      subroutine ehg192(y,d,n,nf,nv,nvmax,vval,lf,lq)
      integer d,i,i1,i2,j,n,nf,nv,nvmax
      integer lq(nvmax,nf)
      DOUBLE PRECISION i3
      DOUBLE PRECISION lf(0:d,nvmax,nf),vval(0:d,nvmax),y(n)

      do 3 i2=1,nv
         do 4 i1=0,d
            vval(i1,i2)=0
    4    continue
    3 continue
      do 5 i=1,nv
         do 6 j=1,nf
            i3=y(lq(i,j))
            do 7 i1=0,d
               vval(i1,i)=vval(i1,i)+i3*lf(i1,i,j)
    7       continue
    6    continue
    5 continue
      return
      end

      DOUBLE PRECISION function ehg128(z,d,ncmax,vc,a,xi,lo,hi,c,v,
     +     nvmax,vval)

c     implicit none
c Args
      integer d,ncmax,nvmax,vc
      integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax)
      DOUBLE PRECISION z(d),xi(ncmax),v(nvmax,d), vval(0:d,nvmax)
c Vars
      logical i2,i3,i4,i5,i6,i7,i8,i9,i10
      integer i,i1,i11,i12,ig,ii,j,lg,ll,m,nt,ur
      integer t(20)
      DOUBLE PRECISION ge,gn,gs,gw,gpe,gpn,gps,gpw,h,phi0,phi1,
     +     psi0,psi1,s,sew,sns,v0,v1,xibar
      DOUBLE PRECISION g(0:8,256),g0(0:8),g1(0:8)

      external loesswarn,ehg184
c     locate enclosing cell
      nt=1
      t(nt)=1
      j=1
c     top of while loop
    3 if(.not.(a(j).ne.0))goto 4
         nt=nt+1
         if(z(a(j)).le.xi(j))then
            i1=lo(j)
         else
            i1=hi(j)
         end if
         t(nt)=i1
         if(.not.(nt.lt.20))then
            call loesswarn(181)
         end if
         j=t(nt)
         goto 3
c     bottom of while loop
    4 continue
c     tensor
      do 5 i12=1,vc
         do 6 i11=0,d
            g(i11,i12)=vval(i11,c(i12,j))
    6    continue
    5 continue
      lg=vc
      ll=c(1,j)
      ur=c(vc,j)
      do 7 i=d,1,-1
         h=(z(i)-v(ll,i))/(v(ur,i)-v(ll,i))
         if(h.lt.-.001D0)then
            call ehg184('eval ',z(1),d,1)
            call ehg184('lowerlimit ',v(ll,1),d,nvmax)
         else
            if(1.001D0.lt.h)then
               call ehg184('eval ',z(1),d,1)
               call ehg184('upperlimit ',v(ur,1),d,nvmax)
            end if
         end if
         if(-.001D0.le.h)then
            i2=(h.le.1.001D0)
         else
            i2=.false.
         end if
         if(.not.i2)then
            call loesswarn(122)
         end if
         lg=int(DBLE(lg)/2.D0)
         do 8 ig=1,lg
c           Hermite basis
            phi0=(1-h)**2*(1+2*h)
            phi1=h**2*(3-2*h)
            psi0=h*(1-h)**2
            psi1=h**2*(h-1)
            g(0,ig)=phi0*g(0,ig) + phi1*g(0,ig+lg) +
     +           (psi0*g(i,ig)+psi1*g(i,ig+lg)) * (v(ur,i)-v(ll,i))
            do 9 ii=1,i-1
               g(ii,ig)=phi0*g(ii,ig)+phi1*g(ii,ig+lg)
    9       continue
    8    continue
    7 continue
      s=g(0,1)
c     blending
      if(d.eq.2)then
c        ----- North -----
         v0=v(ll,1)
         v1=v(ur,1)
         do 10 i11=0,d
            g0(i11)=vval(i11,c(3,j))
   10    continue
         do 11 i11=0,d
            g1(i11)=vval(i11,c(4,j))
   11    continue
         xibar=v(ur,2)
         m=nt-1
c        top of while loop
   12    if(m.eq.0)then
            i4=.true.
         else
            if(a(t(m)).eq.2)then
               i3=(xi(t(m)).eq.xibar)
            else
               i3=.false.
            end if
            i4=i3
         end if
         if(.not.(.not.i4))goto 13
            m=m-1
c           voidp junk
            goto 12
c        bottom of while loop
   13    if(m.ge.1)then
            m=hi(t(m))
c           top of while loop
   14       if(.not.(a(m).ne.0))goto 15
               if(z(a(m)).le.xi(m))then
                  m=lo(m)
               else
                  m=hi(m)
               end if
               goto 14
c           bottom of while loop
   15       if(v0.lt.v(c(1,m),1))then
               v0=v(c(1,m),1)
               do 16 i11=0,d
                  g0(i11)=vval(i11,c(1,m))
   16          continue
            end if
            if(v(c(2,m),1).lt.v1)then
               v1=v(c(2,m),1)
               do 17 i11=0,d
                  g1(i11)=vval(i11,c(2,m))
   17          continue
            end if
         end if
         h=(z(1)-v0)/(v1-v0)
c        Hermite basis
         phi0=(1-h)**2*(1+2*h)
         phi1=h**2*(3-2*h)
         psi0=h*(1-h)**2
         psi1=h**2*(h-1)
         gn=phi0*g0(0)+phi1*g1(0)+(psi0*g0(1)+psi1*g1(1))*(v1-v0)
         gpn=phi0*g0(2)+phi1*g1(2)
c        ----- South -----
         v0=v(ll,1)
         v1=v(ur,1)
         do 18 i11=0,d
            g0(i11)=vval(i11,c(1,j))
   18    continue
         do 19 i11=0,d
            g1(i11)=vval(i11,c(2,j))
   19    continue
         xibar=v(ll,2)
         m=nt-1
c        top of while loop
   20    if(m.eq.0)then
            i6=.true.
         else
            if(a(t(m)).eq.2)then
               i5=(xi(t(m)).eq.xibar)
            else
               i5=.false.
            end if
            i6=i5
         end if
         if(.not.(.not.i6))goto 21
            m=m-1
c           voidp junk
            goto 20
c        bottom of while loop
   21    if(m.ge.1)then
            m=lo(t(m))
c           top of while loop
   22       if(.not.(a(m).ne.0))goto 23
               if(z(a(m)).le.xi(m))then
                  m=lo(m)
               else
                  m=hi(m)
               end if
               goto 22
c           bottom of while loop
   23       if(v0.lt.v(c(3,m),1))then
               v0=v(c(3,m),1)
               do 24 i11=0,d
                  g0(i11)=vval(i11,c(3,m))
   24          continue
            end if
            if(v(c(4,m),1).lt.v1)then
               v1=v(c(4,m),1)
               do 25 i11=0,d
                  g1(i11)=vval(i11,c(4,m))
   25          continue
            end if
         end if
         h=(z(1)-v0)/(v1-v0)
c        Hermite basis
         phi0=(1-h)**2*(1+2*h)
         phi1=h**2*(3-2*h)
         psi0=h*(1-h)**2
         psi1=h**2*(h-1)
         gs=phi0*g0(0)+phi1*g1(0)+(psi0*g0(1)+psi1*g1(1))*(v1-v0)
         gps=phi0*g0(2)+phi1*g1(2)
c        ----- East -----
         v0=v(ll,2)
         v1=v(ur,2)
         do 26 i11=0,d
            g0(i11)=vval(i11,c(2,j))
   26    continue
         do 27 i11=0,d
            g1(i11)=vval(i11,c(4,j))
   27    continue
         xibar=v(ur,1)
         m=nt-1
c        top of while loop
   28    if(m.eq.0)then
            i8=.true.
         else
            if(a(t(m)).eq.1)then
               i7=(xi(t(m)).eq.xibar)
            else
               i7=.false.
            end if
            i8=i7
         end if
         if(.not.(.not.i8))goto 29
            m=m-1
c           voidp junk
            goto 28
c        bottom of while loop
   29    if(m.ge.1)then
            m=hi(t(m))
c           top of while loop
   30       if(.not.(a(m).ne.0))goto 31
               if(z(a(m)).le.xi(m))then
                  m=lo(m)
               else
                  m=hi(m)
               end if
               goto 30
c           bottom of while loop
   31       if(v0.lt.v(c(1,m),2))then
               v0=v(c(1,m),2)
               do 32 i11=0,d
                  g0(i11)=vval(i11,c(1,m))
   32          continue
            end if
            if(v(c(3,m),2).lt.v1)then
               v1=v(c(3,m),2)
               do 33 i11=0,d
                  g1(i11)=vval(i11,c(3,m))
   33          continue
            end if
         end if
         h=(z(2)-v0)/(v1-v0)
c        Hermite basis
         phi0=(1-h)**2*(1+2*h)
         phi1=h**2*(3-2*h)
         psi0=h*(1-h)**2
         psi1=h**2*(h-1)
         ge=phi0*g0(0)+phi1*g1(0)+(psi0*g0(2)+psi1*g1(2))*(v1-v0)
         gpe=phi0*g0(1)+phi1*g1(1)
c        ----- West -----
         v0=v(ll,2)
         v1=v(ur,2)
         do 34 i11=0,d
            g0(i11)=vval(i11,c(1,j))
   34    continue
         do 35 i11=0,d
            g1(i11)=vval(i11,c(3,j))
   35    continue
         xibar=v(ll,1)
         m=nt-1
c        top of while loop
   36    if(m.eq.0)then
            i10=.true.
         else
            if(a(t(m)).eq.1)then
               i9=(xi(t(m)).eq.xibar)
            else
               i9=.false.
            end if
            i10=i9
         end if
         if(.not.(.not.i10))goto 37
            m=m-1
c           voidp junk
            goto 36
c        bottom of while loop
   37    if(m.ge.1)then
            m=lo(t(m))
c           top of while loop
   38       if(.not.(a(m).ne.0))goto 39
               if(z(a(m)).le.xi(m))then
                  m=lo(m)
               else
                  m=hi(m)
               end if
               goto 38
c           bottom of while loop
   39       if(v0.lt.v(c(2,m),2))then
               v0=v(c(2,m),2)
               do 40 i11=0,d
                  g0(i11)=vval(i11,c(2,m))
   40          continue
            end if
            if(v(c(4,m),2).lt.v1)then
               v1=v(c(4,m),2)
               do 41 i11=0,d
                  g1(i11)=vval(i11,c(4,m))
   41          continue
            end if
         end if
         h=(z(2)-v0)/(v1-v0)
c        Hermite basis
         phi0=(1-h)**2*(1+2*h)
         phi1=h**2*(3-2*h)
         psi0=h*(1-h)**2
         psi1=h**2*(h-1)
         gw=phi0*g0(0)+phi1*g1(0)+(psi0*g0(2)+psi1*g1(2))*(v1-v0)
         gpw=phi0*g0(1)+phi1*g1(1)
c        NS
         h=(z(2)-v(ll,2))/(v(ur,2)-v(ll,2))
c        Hermite basis
         phi0=(1-h)**2*(1+2*h)
         phi1=h**2*(3-2*h)
         psi0=h*(1-h)**2
         psi1=h**2*(h-1)
         sns=phi0*gs+phi1*gn+(psi0*gps+psi1*gpn)*(v(ur,2)-v(ll,2))
c        EW
         h=(z(1)-v(ll,1))/(v(ur,1)-v(ll,1))
c        Hermite basis
         phi0=(1-h)**2*(1+2*h)
         phi1=h**2*(3-2*h)
         psi0=h*(1-h)**2
         psi1=h**2*(h-1)
         sew=phi0*gw+phi1*ge+(psi0*gpw+psi1*gpe)*(v(ur,1)-v(ll,1))
         s=(sns+sew)-s
      end if
      ehg128=s
      return
      end

      integer function ifloor(x)
      DOUBLE PRECISION x
      ifloor=int(x)
      if(ifloor.gt.x) ifloor=ifloor-1
      end

c DSIGN is unused, causes conflicts on some platforms
c       DOUBLE PRECISION function DSIGN(a1,a2)
c       DOUBLE PRECISION a1, a2
c       DSIGN=DABS(a1)
c       if(a2.ge.0)DSIGN=-DSIGN
c       end


c ehg136()  is the workhorse of lowesf(.)
c     n = number of observations
c     m = number of x values at which to evaluate
c     f = span
c     nf = min(n, floor(f * n))
      subroutine ehg136(u,lm,m,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,
     +     od,o,ihat,w,rcond,sing,dd,tdeg,cdeg,s)
      integer identi,d,dd,i,i1,ihat,info,j,k,kernel,l,lm,m,n,nf,
     +     od,sing,tdeg
      integer cdeg(8),psi(n)
      double precision f,i2,rcond,scale,tol
      double precision o(m,n),sigma(15),e(15,15),g(15,15),b(nf,k),
     $     dist(n),eta(nf),dgamma(15),q(8),qraux(15),rw(n),s(0:od,m),
     $     u(lm,d),w(nf),work(15),x(n,d),y(n)

      external ehg127,loesswarn,dqrsl
      double precision ddot
      external ddot

c     V -> g
c     U -> e
c     Identity -> identi
c     L -> o
c     X -> b
      if(k .gt. nf-1) call loesswarn(104)
      if(k .gt. 15)   call loesswarn(105)
      do 3 identi=1,n
         psi(identi)=identi
    3 continue
      do 4 l=1,m
         do 5 i1=1,d
            q(i1)=u(l,i1)
    5    continue
         call ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,
     +        rcond,sing,sigma,e,g,dgamma,qraux,work,tol,dd,tdeg,cdeg,
     +        s(0,l))
         if(ihat.eq.1)then
c           $L sub {l,l} =
c           V sub {1,:} SIGMA sup {+} U sup T
c           (Q sup T W e sub i )$
            if(.not.(m.eq.n))then
               call loesswarn(123)
            end if
c           find $i$ such that $l = psi sub i$
            i=1
c           top of while loop
    6       if(.not.(l.ne.psi(i)))goto 7
               i=i+1
               if(.not.(i.lt.nf))then
                  call loesswarn(123)
c next line is not in current dloess
                  goto 7
               end if
               goto 6
c           bottom of while loop
    7       do 8 i1=1,nf
               eta(i1)=0
    8       continue
            eta(i)=w(i)
c           $eta = Q sup T W e sub i$
            call dqrsl(b,nf,nf,k,qraux,eta,eta,eta,eta,eta,eta,1000,
     +           info)
c           $gamma = U sup T eta sub {1:k}$
            do 9 i1=1,k
               dgamma(i1)=0
    9       continue
            do 10 j=1,k
               i2=eta(j)
               do 11 i1=1,k
                  dgamma(i1)=dgamma(i1)+i2*e(j,i1)
   11          continue
   10       continue
c           $gamma = SIGMA sup {+} gamma$
            do 12 j=1,k
               if(tol.lt.sigma(j))then
                  dgamma(j)=dgamma(j)/sigma(j)
               else
                  dgamma(j)=0.d0
               end if
   12       continue
c           voidp junk
c           voidp junk
            o(l,1)=ddot(k,g(1,1),15,dgamma,1)
         else
            if(ihat.eq.2)then
c              $L sub {l,:} =
c              V sub {1,:} SIGMA sup {+}
c              ( U sup T Q sup T ) W $
               do 13 i1=1,n
                  o(l,i1)=0
   13          continue
               do 14 j=1,k
                  do 15 i1=1,nf
                     eta(i1)=0
   15             continue
                  do 16 i1=1,k
                     eta(i1)=e(i1,j)
   16             continue
                  call dqrsl(b,nf,nf,k,qraux,eta,eta,work,work,work,work
     +                 ,10000,info)
                  if(tol.lt.sigma(j))then
                     scale=1.d0/sigma(j)
                  else
                     scale=0.d0
                  end if
                  do 17 i1=1,nf
                     eta(i1)=eta(i1)*(scale*w(i1))
   17             continue
                  do 18 i=1,nf
                     o(l,psi(i))=o(l,psi(i))+g(1,j)*eta(i)
   18             continue
   14          continue
            end if
         end if
    4 continue
      return
      end

c called from lowesb() ... compute fit ..?..?...
c somewhat similar to ehg136
      subroutine ehg139(v,nvmax,nv,n,d,nf,f,x,pi,psi,y,rw,trl,kernel,k,
     +     dist,phi,eta,b,od,w,diagl,vval2,ncmax,vc,a,xi,lo,hi,c,
     +     rcond,sing,dd,tdeg,cdeg,lq,lf,setlf,s)
      logical setlf
      integer identi,d,dd,i,i2,i3,i5,i6,ii,ileaf,info,j,k,kernel,
     +     l,n,ncmax,nf,nleaf,nv,nvmax,od,sing,tdeg,vc
      integer lq(nvmax,nf),a(ncmax),c(vc,ncmax),cdeg(8),hi(ncmax),
     +     leaf(256),lo(ncmax),pi(n),psi(n)
      DOUBLE PRECISION f,i1,i4,i7,rcond,scale,term,tol,trl
      DOUBLE PRECISION lf(0:d,nvmax,nf),sigma(15),u(15,15),e(15,15),
     +     b(nf,k),diagl(n),dist(n),eta(nf),DGAMMA(15),q(8),qraux(15),
     +     rw(n),s(0:od,nv),v(nvmax,d),vval2(0:d,nv),w(nf),work(15),
     +     x(n,d),xi(ncmax),y(n),z(8)
      DOUBLE PRECISION phi(n)

      external ehg127,loesswarn,DQRSL,ehg137
      DOUBLE PRECISION ehg128
      external ehg128
      DOUBLE PRECISION DDOT
      external DDOT

c     V -> e
c     Identity -> identi
c     X -> b
c     l2fit with trace(L)
      if(k .gt. nf-1) call loesswarn(104)
      if(k .gt. 15)   call loesswarn(105)
      if(trl.ne.0)then
         do 3 i5=1,n
            diagl(i5)=0
    3    continue
         do 4 i6=1,nv
            do 5 i5=0,d
               vval2(i5,i6)=0
    5       continue
    4    continue
      end if
      do 6 identi=1,n
         psi(identi)=identi
    6 continue
      do 7 l=1,nv
         do 8 i5=1,d
            q(i5)=v(l,i5)
    8    continue
         call ehg127(q,n,d,nf,f,x,psi,y,rw,kernel,k,dist,eta,b,od,w,
     +        rcond,sing,sigma,u,e,DGAMMA,qraux,work,tol,dd,tdeg,cdeg,
     +        s(0,l))
         if(trl.ne.0)then
c           invert $psi$
            do 9 i5=1,n
               phi(i5)=0
    9       continue
            do 10 i=1,nf
               phi(psi(i))=i
   10       continue
            do 11 i5=1,d
               z(i5)=v(l,i5)
   11       continue
            call ehg137(z,leaf,nleaf,d,ncmax,a,xi,
     +           lo,hi)
            do 12 ileaf=1,nleaf
               do 13 ii=lo(leaf(ileaf)),hi(leaf(ileaf))
                  i=int(phi(pi(ii)))
                  if(i.ne.0)then
                     if(.not.(psi(i).eq.pi(ii)))then
                        call loesswarn(194)
                     end if
                     do 14 i5=1,nf
                        eta(i5)=0
   14                continue
                     eta(i)=w(i)
c                    $eta = Q sup T W e sub i$
                     call DQRSL(b,nf,nf,k,qraux,eta,work,eta,eta,work,
     +                    work,1000,info)
                     do 15 j=1,k
                        if(tol.lt.sigma(j))then
                           i4=DDOT(k,u(1,j),1,eta,1)/sigma(j)
                        else
                           i4=0.D0
                        end if
                       DGAMMA(j)=i4
   15                continue
                     do 16 j=1,d+1
c                       bug fix 2006-07-15 for k=1, od>1.   (thanks btyner@gmail.com)
                        if(j.le.k)then
                           vval2(j-1,l)=DDOT(k,e(j,1),15,DGAMMA,1)
                        else
                           vval2(j-1,l)=0.0d0
                        end if
   16                continue
                     do 17 i5=1,d
                        z(i5)=x(pi(ii),i5)
   17                continue
                     term=ehg128(z,d,ncmax,vc,a,xi,lo,hi,c,v,nvmax,
     +                    vval2)
                     diagl(pi(ii))=diagl(pi(ii))+term
                     do 18 i5=0,d
                        vval2(i5,l)=0
   18                continue
                  end if
   13          continue
   12       continue
         end if
         if(setlf)then
c           $Lf sub {:,l,:} = V SIGMA sup {+} U sup T Q sup T W$
            if(.not.(k.ge.d+1))then
               call loesswarn(196)
            end if
            do 19 i5=1,nf
               lq(l,i5)=psi(i5)
   19       continue
            do 20 i6=1,nf
               do 21 i5=0,d
                  lf(i5,l,i6)=0
   21          continue
   20       continue
            do 22 j=1,k
               do 23 i5=1,nf
                  eta(i5)=0
   23          continue
               do 24 i5=1,k
                  eta(i5)=u(i5,j)
   24          continue
               call DQRSL(b,nf,nf,k,qraux,eta,eta,work,work,work,work,
     +              10000,info)
               if(tol.lt.sigma(j))then
                  scale=1.D0/sigma(j)
               else
                  scale=0.D0
               end if
               do 25 i5=1,nf
                  eta(i5)=eta(i5)*(scale*w(i5))
   25          continue
               do 26 i=1,nf
                  i7=eta(i)
                  do 27 i5=0,d
                     if(i5.lt.k)then
                        lf(i5,l,i)=lf(i5,l,i)+e(1+i5,j)*i7
                     else
                        lf(i5,l,i)=0
                     end if
   27             continue
   26          continue
   22       continue
         end if
    7 continue
      if(trl.ne.0)then
         if(n.le.0)then
            trl=0.D0
         else
            i3=n
            i1=diagl(i3)
            do 28 i2=i3-1,1,-1
               i1=diagl(i2)+i1
   28       continue
            trl=i1
         end if
      end if
      return
      end

      subroutine lowesb(xx,yy,ww,diagl,infl,iv,wv)
      integer infl
      integer iv(*)
      DOUBLE PRECISION xx(*),yy(*),ww(*),diagl(*),wv(*)
c Var
      DOUBLE PRECISION trl
      logical setlf

      integer ifloor
      external ifloor
      external ehg131,loesswarn,ehg183

      if(.not.(iv(28).ne.173))then
         call loesswarn(174)
      end if
      if(iv(28).ne.172)then
         if(.not.(iv(28).eq.171))then
            call loesswarn(171)
         end if
      end if
      iv(28)=173
      if(infl.ne.0)then
         trl=1.D0
      else
         trl=0.D0
      end if
      setlf=(iv(27).ne.iv(25))
      call ehg131(xx,yy,ww,trl,diagl,iv(20),iv(29),iv(3),iv(2),iv(5),
     +     iv(17),iv(4),iv(6),iv(14),iv(19),wv(1),iv(iv(7)),iv(iv(8)),
     +     iv(iv(9)),iv(iv(10)),iv(iv(22)),iv(iv(27)),wv(iv(11)),
     +     iv(iv(23)),wv(iv(13)),wv(iv(12)),wv(iv(15)),wv(iv(16)),
     +     wv(iv(18)),ifloor(iv(3)*wv(2)),wv(3),wv(iv(26)),wv(iv(24)),
     +     wv(4),iv(30),iv(33),iv(32),iv(41),iv(iv(25)),wv(iv(34)),
     +     setlf)
      if(iv(14).lt.iv(6)+DBLE(iv(4))/2.D0)then
         call ehg183('k-d tree limited by memory; nvmax=',
     +        iv(14),1,1)
      else
         if(iv(17).lt.iv(5)+2)then
            call ehg183('k-d tree limited by memory. ncmax=',
     +           iv(17),1,1)
         end if
      end if
      return
      end

c lowesd() : Initialize iv(*) and v(1:4)
c ------     called only by loess_workspace()  in ./loessc.c
      subroutine lowesd(iv, liv,lv, v, d,n,f,ideg,nf,nvmax, setlf)
      integer               liv,lv,    d,n,  ideg,nf,nvmax, setlf
c           setlf if(true) need  L [nf x nvmax] matrices
      integer iv(liv)
      double precision f, v(lv)

c  had  nf = min(n,ifloor(n*f))

      integer bound,i,i1,i2,j,ncmax,vc
      external loesswarn
      integer ifloor
      external ifloor
c
c     unnecessary initialization of i1 to keep g77 -Wall happy
c
      i1 = 0
cc    version -> versio
c     if(.not.(versio.eq.106))then
c        call loesswarn(100)
c     end if
      iv(28)=171
      iv(2)=d
      iv(3)=n
      vc=2**d
      iv(4)=vc
      if(.not.(0.lt.f))then
         call loesswarn(120)
      end if

      iv(19)=nf
      iv(20)=1
      if(ideg.eq.0)then
         i1=1
      else
         if(ideg.eq.1)then
            i1=d+1
         else
            if(ideg.eq.2)then
               i1=int(dble((d+2)*(d+1))/2.d0)
            end if
         end if
      end if
      iv(29)=i1
      iv(21)=1
      iv(14)=nvmax
      ncmax=nvmax
      iv(17)=ncmax
      iv(30)=0
      iv(32)=ideg
      if(.not.(ideg.ge.0))then
         call loesswarn(195)
      end if
      if(.not.(ideg.le.2))then
         call loesswarn(195)
      end if
      iv(33)=d
      do 3 i2=41,49
         iv(i2)=ideg
    3 continue
      iv(7)=50
      iv(8)=iv(7)+ncmax
      iv(9)=iv(8)+vc*ncmax
      iv(10)=iv(9)+ncmax
      iv(22)=iv(10)+ncmax
c     initialize permutation
      j=iv(22)-1
      do 4 i=1,n
         iv(j+i)=i
    4 continue
      iv(23)=iv(22)+n
      iv(25)=iv(23)+nvmax
      if(setlf.ne.0)then
         iv(27)=iv(25)+nvmax*nf
      else
         iv(27)=iv(25)
      end if
      bound=iv(27)+n
      if(.not.(bound-1.le.liv))then
         call loesswarn(102)
      end if
      iv(11)=50
      iv(13)=iv(11)+nvmax*d
      iv(12)=iv(13)+(d+1)*nvmax
      iv(15)=iv(12)+ncmax
      iv(16)=iv(15)+n
      iv(18)=iv(16)+nf
      iv(24)=iv(18)+iv(29)*nf
      iv(34)=iv(24)+(d+1)*nvmax
      if(setlf.ne.0)then
         iv(26)=iv(34)+(d+1)*nvmax*nf
      else
         iv(26)=iv(34)
      end if
      bound=iv(26)+nf
      if(.not.(bound-1.le.lv))then
         call loesswarn(103)
      end if
      v(1)=f
      v(2)=0.05d0
      v(3)=0.d0
      v(4)=1.d0
      return
      end

      subroutine lowese(iv,wv,m,z,s)
      integer m
      integer iv(*)
      double precision s(m),wv(*),z(m,1)

      external ehg133,loesswarn

      if(.not.(iv(28).ne.172))then
         call loesswarn(172)
      end if
      if(.not.(iv(28).eq.173))then
         call loesswarn(173)
      end if
      call ehg133(iv(2),iv(4),iv(14),iv(17),iv(iv(7)),iv(iv(
     +8)),iv(iv(9)),iv(iv(10)),wv(iv(11)),wv(iv(13)),wv(iv(12)),m,z,s)
      return
      end

c "direct" (non-"interpolate") fit aka predict() :
      subroutine lowesf(xx,yy,ww,iv,wv,m,z,l,ihat,s)
      integer m,ihat
c     m = number of x values at which to evaluate
      integer iv(*)
      double precision xx(*),yy(*),ww(*),wv(*),z(m,1),l(m,*),s(m)

      logical i1

      external loesswarn,ehg136
      if(171.le.iv(28))then
         i1=(iv(28).le.174)
      else
         i1=.false.
      end if
      if(.not.i1)then
         call loesswarn(171)
      end if
      iv(28)=172
      if(.not.(iv(14).ge.iv(19)))then
         call loesswarn(186)
      end if

c do the work; in ehg136()  give the argument names as they are there:
c          ehg136(u,lm,m, n,    d,    nf,   f,   x,   psi,     y ,rw,
      call ehg136(z,m,m,iv(3),iv(2),iv(19),wv(1),xx,iv(iv(22)),yy,ww,
c          kernel,  k,     dist,       eta,       b,     od,o,ihat,
     +     iv(20),iv(29),wv(iv(15)),wv(iv(16)),wv(iv(18)),0,l,ihat,
c              w,     rcond,sing,    dd,    tdeg,cdeg,  s)
     +     wv(iv(26)),wv(4),iv(30),iv(33),iv(32),iv(41),s)
      return
      end

c Called either from loess_raw() only for case [surf_stat = "interpolate/exact"], or
c               from loess_ise() {used only when 'se = TRUE' and surface = "interpolate"}
      subroutine lowesl(iv,wv,m,z,l)
      integer m
      integer iv(*)
      double precision l(m,*), wv(*), z(m,1)

      external loesswarn,ehg191

      if(.not.(iv(28).ne.172))then
         call loesswarn(172)
      end if
      if(.not.(iv(28).eq.173))then
         call loesswarn(173)
      end if
      if(.not.(iv(26).ne.iv(34)))then
         call loesswarn(175)
      end if
      call ehg191(m,z,l,iv(2),iv(3),iv(19),iv(6),iv(17),iv(4),iv(iv(7)),
     +     wv(iv(12)),iv(iv(10)),iv(iv(9)),iv(iv(8)),wv(iv(11)),iv(14),
     +     wv(iv(24)),wv(iv(34)),iv(iv(25)))
      return
      end

c  Not used
      subroutine lowesr(yy,iv,wv)
      integer iv(*)
      DOUBLE PRECISION yy(*),wv(*)

      external loesswarn,ehg192
      if(.not.(iv(28).ne.172))then
         call loesswarn(172)
      end if
      if(.not.(iv(28).eq.173))then
         call loesswarn(173)
      end if
      call ehg192(yy,iv(2),iv(3),iv(19),iv(6),iv(14),wv(iv(13)),
     +     wv(iv(34)),iv(iv(25)))
      return
      end

c Update Robustness Weights -- called via .Fortran() from R's simpleLoess() in ../R/loess.R
c
      subroutine lowesw(res,n,rw,pi)
c Tranliterated from Devlin's ratfor

c     implicit none
c Args
      integer n
      double precision res(n), rw(n)
      integer pi(n)
c Var
      integer identi,i,i1,nh
      double precision cmad,rsmall

      integer ifloor
      double precision d1mach

      external ehg106
      external ifloor
      external d1mach

c     Identity -> identi

c     find median of absolute residuals
      do 3 i1=1,n
         rw(i1)=dabs(res(i1))
    3 continue
      do 4 identi=1,n
         pi(identi)=identi
    4 continue
      nh=ifloor(dble(n)/2.d0)+1
c     partial sort to find 6*mad
      call ehg106(1,n,nh,1,rw,pi,n)
      if((n-nh)+1.lt.nh)then
         call ehg106(1,nh-1,nh-1,1,rw,pi,n)
         cmad=3*(rw(pi(nh))+rw(pi(nh-1)))
      else
         cmad=6*rw(pi(nh))
      end if
      rsmall=d1mach(1)
      if(cmad.lt.rsmall)then
         do 5 i1=1,n
            rw(i1)=1
    5    continue
      else
         do 6 i=1,n
            if(cmad*0.999d0.lt.rw(i))then
               rw(i)=0
            else
               if(cmad*0.001d0.lt.rw(i))then
                  rw(i)=(1-(rw(i)/cmad)**2)**2
               else
                  rw(i)=1
               end if
            end if
    6    continue
      end if
      return
      end

c Compute Pseudovalues -- called via .Fortran() from R's simpleLoess() in ../R/loess.R
c                         in the case of robustness iterations
      subroutine lowesp(n,y,yhat,pwgts,rwgts,pi,ytilde)
      integer n
      integer pi(n)
      double precision y(n),yhat(n),pwgts(n),rwgts(n),ytilde(n)
c Var
      double precision c,i1,i4,mad
      integer i2,i3,i,m

      external ehg106
      integer ifloor
      external ifloor
c     median absolute deviation (using partial sort):
      do 3 i=1,n
         ytilde(i)=dabs(y(i)-yhat(i))*dsqrt(pwgts(i))
         pi(i) = i
    3 continue
      m=ifloor(dble(n)/2.d0)+1
      call ehg106(1,n,m,1,ytilde,pi,n)
      if((n-m)+1.lt.m)then
         call ehg106(1,m-1,m-1,1,ytilde,pi,n)
         mad=(ytilde(pi(m-1))+ytilde(pi(m)))/2
      else
         mad=ytilde(pi(m))
      end if
c     magic constant
      c=(6*mad)**2/5
      do 5 i=1,n
         ytilde(i)= 1 - ((y(i)-yhat(i))**2 * pwgts(i))/c
    5 continue
      do 6 i=1,n
         ytilde(i)=ytilde(i)*dsqrt(rwgts(i))
    6 continue
      if(n.le.0)then
         i4=0.d0
      else
         i3=n
         i1=ytilde(i3)
         do 7 i2=i3-1,1,-1
            i1=ytilde(i2)+i1
    7    continue
         i4=i1
      end if
      c=n/i4
c     pseudovalues
      do 8 i=1,n
         ytilde(i)=yhat(i) + (c*rwgts(i))*(y(i)-yhat(i))
    8 continue
      return
      end

      subroutine ehg124(ll,uu,d,n,nv,nc,ncmax,vc,x,pi,a,xi,lo,hi,c,v,
     +     vhit,nvmax,fc,fd,dd)

      integer ll,uu,d,n,nv,nc,ncmax,vc,nvmax,fc,dd
      integer a(ncmax),c(vc,ncmax),hi(ncmax),lo(ncmax),pi(n),vhit(nvmax)
      DOUBLE PRECISION fd, v(nvmax,d),x(n,d),xi(ncmax)

      logical i1,i2,leaf
      integer i4,inorm2,k,l,m,p,u, upper, lower, check, offset
      DOUBLE PRECISION diam,diag(8),sigma(8)

      external ehg125,ehg106,ehg129
      integer IDAMAX
      external IDAMAX
      p=1
      l=ll
      u=uu
      lo(p)=l
      hi(p)=u
c     top of while loop
    3 if(.not.(p.le.nc))goto 4
         do 5 i4=1,dd
            diag(i4)=v(c(vc,p),i4)-v(c(1,p),i4)
    5    continue
         diam=0
         do 6 inorm2=1,dd
            diam=diam+diag(inorm2)**2
    6    continue
         diam=DSQRT(diam)
         if((u-l)+1.le.fc)then
            i1=.true.
         else
            i1=(diam.le.fd)
         end if
         if(i1)then
            leaf=.true.
         else
            if(ncmax.lt.nc+2)then
               i2=.true.
            else
               i2=(nvmax.lt.nv+DBLE(vc)/2.D0)
            end if
            leaf=i2
         end if
         if(.not.leaf)then
            call ehg129(l,u,dd,x,pi,n,sigma)
            k=IDAMAX(dd,sigma,1)
            m=int(DBLE(l+u)/2.D0)
            call ehg106(l,u,m,1,x(1,k),pi,n)

c           all ties go with hi son
c           top of while loop
c     bug fix from btyner@gmail.com 2006-07-20
            offset = 0
 7          if(((m+offset).ge.u).or.((m+offset).lt.l))goto 8
            if(offset .lt. 0)then
               lower = l
               check = m + offset
               upper = check
            else
               lower = m + offset + 1
               check = lower
               upper = u
            end if
            call ehg106(lower,upper,check,1,x(1,k),pi,n)
            if(x(pi(m + offset),k).eq.x(pi(m+offset+1),k))then
               offset = -offset
               if(offset .ge. 0) then
                  offset = offset + 1
               end if
               goto 7
            else
               m = m + offset
               goto 8
            end if

c           bottom of while loop
    8       if(v(c(1,p),k).eq.x(pi(m),k))then
               leaf=.true.
            else
               leaf=(v(c(vc,p),k).eq.x(pi(m),k))
            end if
         end if
         if(leaf)then
            a(p)=0
         else
            a(p)=k
            xi(p)=x(pi(m),k)
c           left son
            nc=nc+1
            lo(p)=nc
            lo(nc)=l
            hi(nc)=m
c           right son
            nc=nc+1
            hi(p)=nc
            lo(nc)=m+1
            hi(nc)=u
            call ehg125(p,nv,v,vhit,nvmax,d,k,xi(p),2**(k-1),2**(d-k),
     +                  c(1,p),c(1,lo(p)),c(1,hi(p)))
         end if
         p=p+1
         l=lo(p)
         u=hi(p)
         goto 3
c     bottom of while loop
    4 return
      end

      subroutine ehg129(l,u,d,x,pi,n,sigma)
      integer d,execnt,i,k,l,n,u
      integer pi(n)
      DOUBLE PRECISION machin,alpha,beta,t
      DOUBLE PRECISION sigma(d),x(n,d)
      DOUBLE PRECISION D1MACH
      external D1MACH
      save machin,execnt
      data execnt /0/
c     MachInf -> machin
      execnt=execnt+1
      if(execnt.eq.1)then
c     initialize  d1mach(2) === DBL_MAX:
         machin=D1MACH(2)
      end if
      do 3 k=1,d
         alpha=machin
         beta=-machin
         do 4 i=l,u
            t=x(pi(i),k)
            alpha=min(alpha,x(pi(i),k))
            beta=max(beta,t)
    4    continue
         sigma(k)=beta-alpha
    3 continue
      return
      end

c {called only from ehg127}  purpose...?...
      subroutine ehg137(z,leaf,nleaf,d,ncmax,a,xi,lo,hi)
      integer d,nleaf,ncmax
      integer leaf(256),a(ncmax),hi(ncmax),lo(ncmax),pstack(20)
      DOUBLE PRECISION z(d),xi(ncmax)

      integer p,stackt

      external loesswarn
c     stacktop -> stackt
c     find leaf cells affected by $z$
      stackt=0
      p=1
      nleaf=0
c     top of while loop
    3 if(.not.(0.lt.p))goto 4
         if(a(p).eq.0)then
c           leaf
            nleaf=nleaf+1
            leaf(nleaf)=p
c           Pop
            if(stackt.ge.1)then
               p=pstack(stackt)
            else
               p=0
            end if
            stackt=max(0,stackt-1)
         else
            if(z(a(p)).eq.xi(p))then
c              Push
               stackt=stackt+1
               if(.not.(stackt.le.20))then
                  call loesswarn(187)
               end if
               pstack(stackt)=hi(p)
               p=lo(p)
            else
               if(z(a(p)).le.xi(p))then
                  p=lo(p)
               else
                  p=hi(p)
               end if
            end if
         end if
         goto 3
c     bottom of while loop
    4 if(.not.(nleaf.le.256))then
         call loesswarn(185)
      end if
      return
      end

C-- For Error messaging, call the "a" routines at the bottom of  ./loessc.c  :
      subroutine ehg183(s, i, n, inc)
      character s*(*)
      integer i, n, inc
      call ehg183a(s, len(s), i, n, inc)
      end

      subroutine ehg184(s, x, n, inc)
      character s*(*)
      double precision x
      integer n, inc
      call ehg184a(s, len(s), x, n, inc)
      end