subroutine ppdeux (mm,p,n,jj,x,w,trm,nei,fei,maxit,thr,itape,a,xpa
     1,sp,dp)
c
c two-dimensional exploratory projection pursuit (pp).
c
c Friedman, J. H. (1987). Exploratory projection pursuit. J. Amer. Statist.
c    Assoc., 82, 249-266 (March).
c
c version of 3/10/87.
c
c coded and copywrite (c) by:
c
c                             Jerome H. Friedman
c
c                           Department of Statistics
c                                     and
c                        Stanford Linear Accelerator Center
c                             Stanford University
c
c input:
c
c mm = number of pp solutions (pp followed by structure removal).
c p = number of variables (integer).
c n = number of observations.
c jj = order of the legenre polynomial expansion (for transformed
c      density estimates). typical values, jj=2 to jj=6, depending
c      on the values of p and n.
c x(p,n) = data matrix.
c w(n) = weight (mass) for each observation.
c trm = trimming threshold for robust sphering. observations are
c       ignored in the sphering calculation (only) if their
c       mahalanobis distance from the mean is greater than trm.
c       the value trm = 0.0 indicates no robustification.
c       typical values are trm =0.0, or trm = the 1.0 - nout/n point
c       of a chi-squared distribution on p degrees-of-freedom,
c       with nout a small number.
c nei = maximum dimensionality of search space. solutions are
c       constrained to lie in the subspace spanned by the largest
c       nei principal components. typical value, nei = p.
c fei = threshold for data space dimensionality reduction.
c       solutions are constrained to lie in the subspace spanned by
c       those eigenvectors of the covariance matrix whose eigenvalues
c       are larger than fei times the largest eigenvalue. typical
c       value, fei = 1.e-4.
c maxit = maximum number of iterations in the gradient directed part
c         of the numerical optimization. typical value, maxit = 10.
c thr = convergence threshold for gradient directed part of the
c       numerical optimization. maximization 'converges' when
c       improvement from the previous iteration is less than thr
c       times its value at the current iteration. typical value,
c       thr = 0.01.
c itape = fortran unit number to which output is directed.
c         typical value, itape = 6.
c
c output:
c
c a(p,2,mm) = coefficients of the two linear combinations defining
c             each of the mm pp solutions.
c xpa(2,n,mm) = adjusted data plots for each of the mm pp solutions.
c               data distribution (abscissa and ordinate for each
c               observation ) projected onto each pp solution
c               with the structure associated with all previous
c               solutions removed.
c
c workspace:
c
c sp(n*(p+2*jj+4)+jj*(jj+2)+3*p) : single precision.
c dp(2*p**2+7*p) : double precision.
c
c note :
c
c the typical values noted above are suggestions based on
c limited experience.
c
c this is a fairly computationally intensive procedure, depending
c on problem size and input parameter values. a vax or vax equi-
c valent (1 mip) computer is the smallest on which it would
c reasonably run.
c
      integer p
      real x(p,n),w(n),a(p,2,mm),xpa(2,n,mm),sp(1)
      double precision dp(1),r,s,t,u,v
      common /debug/ ipt
      j1=1
      j2=j1+p
      j3=j2+p*n
      j4=j3+jj*2*n
      j5=j4+2*jj
      j6=j5+2*p
      j7=j6+2*n
      j8=j7+2*n
      k1=1
      k2=k1+p**2
      k3=k2+p
      k4=k3+2*p
      k5=k4+p*2
      k6=k5+p**2
      k7=k6+p
      iter=0
      nst=iter
      l=j2-1
 10   iter=iter+1
      call sphere (p,n,x,w,sw,fei,nei,itape,sp(j1),dp(k2),dp(k1),nj,sp(j
     12),dp(k6))
      if (iter.ge.maxit.or.trm.le.0.0) go to 40
      ns=0
      do 30 j=1,n
      if (w(j).le.0.0) go to 30
      s=0.0
      do 20 i=1,nj
      s=s+sp(p*(j-1)+i+l)**2
 20   continue
      if (s.le.trm) go to 30
      w(j)=0.0
      ns=ns+1
 30   continue
      if (ns.eq.0) go to 40
      nst=nst+ns
      go to 10
 40   write (itape,720) nst,iter,nj
      if (nj.lt.2) return
      do 710 im=1,mm
      write (itape,750) im
      hi=0.0
      l=j2-1
      l4=k4-1
      njm1=nj-1
      do 70 i=1,njm1
      ip1=i+1
      do 60 k=ip1,nj
      do 50 j=1,n
      xpa(1,j,im)=sp(p*(j-1)+i+l)
      xpa(2,j,im)=sp(p*(j-1)+k+l)
 50   continue
      call legndr (jj,p,nj,n,sp(j2),w,a(1,1,im),-1,sw,xpa(1,1,im),sp(j7)
     1,hj,dp(k3),sp(j3),sp(j4),sp(j8))
      if (hj.ge.hi) go to 60
      hi=hj
      it=i
      kt=k
 60   continue
 70   continue
      do 80 j=1,n
      xpa(1,j,im)=sp(p*(j-1)+it+l)
      xpa(2,j,im)=sp(p*(j-1)+kt+l)
 80   continue
      do 90 i=1,nj
      a(i,1,im)=0.0
      a(i,2,im)=0.0
 90   continue
      a(it,1,im)=1.0
      a(kt,2,im)=1.0
      if (ipt.lt.3) go to 110
      write (6,100) hi,it,kt
 100  format (    5h0hi =g12.4,    8h comps =2i5/)
 110  l6=j6-1
      do 120 i=1,nj
      dp(i+l4)=1.0
      dp(i+l4+nj)=1.0
 120  continue
 130  hi0=hi
      do 210 kk=1,2
      ko=2
      if (kk.eq.2) ko=1
      do 140 j=1,n
      sp(l6+2*(j-1)+ko)=xpa(ko,j,im)
 140  continue
      do 200 i=1,nj
      if (a(i,kk,im).ge.1..or.a(i,ko,im).ge.1.) go to 200
      sgn=1.
      ang=1./sqrt(2.*(1.+a(i,kk,im))-a(i,ko,im)**2)
      fang=ang
      do 150 j=1,n
      sp(l6+2*(j-1)+kk)=ang*(xpa(kk,j,im)+sp(p*(j-1)+i+l)-a(i,ko,im)*xpa
     1(ko,j,im))
 150  continue
      call legndr (jj,p,nj,n,sp(j2),w,a(1,1,im),-1,sw,sp(j6),sp(j7),hk,d
     1p(k3),sp(j3),sp(j4),sp(j8))
      bang=1./sqrt(2.*(1.-a(i,kk,im))-a(i,ko,im)**2)
      do 160 j=1,n
      sp(l6+2*(j-1)+kk)=bang*(xpa(kk,j,im)-sp(p*(j-1)+i+l)+a(i,ko,im)*xp
     1a(ko,j,im))
 160  continue
      call legndr (jj,p,nj,n,sp(j2),w,a(1,1,im),-1,sw,sp(j6),sp(j7),hj,d
     1p(k3),sp(j3),sp(j4),sp(j8))
      if (hj.ge.hk) go to 170
      hk=hj
      sgn=-1.0
      fang=bang
 170  if (hk.ge.hi) go to 200
      hi=hk
      a(i,kk,im)=a(i,kk,im)+sgn
      do 180 j=1,nj
      a(j,kk,im)=fang*(a(j,kk,im)-sgn*a(i,ko,im)*a(j,ko,im))
 180  continue
      do 190 j=1,n
      xpa(kk,j,im)=fang*(xpa(kk,j,im)+sgn*sp(p*(j-1)+i+l)-sgn*a(i,ko,im)
     1*xpa(ko,j,im))
 190  continue
 200  continue
 210  continue
      if (ipt.lt.3) go to 250
      write (6,220) hi
 220  format (    5h0hi =g12.4/)
      write (6,230) (a(i,1,im),i=1,nj)
 230  format (    4h a =2(/    1h 5g12.4))
      write (6,240) (a(i,2,im),i=1,nj)
 240  format (    4h b =2(/    1h 5g12.4))
 250  if ((hi-hi0)/hi.gt.0.1*thr) go to 130
      iter=0
      iexit=iter
      l6=k6-1
      l7=k7-1
      l3=k3-1
      g0=0.0
      r0=1.0
 260  iter=iter+1
      call legndr (jj,p,nj,n,sp(j2),w,a(1,1,im),1,sw,xpa(1,1,im),sp(j7),
     1hi,dp(k3),sp(j3),sp(j4),sp(j8))
      if (ipt.lt.2) go to 280
      write (6,270) r0,hi
 270  format (    7h r,hj =3g12.4)
 280  gi=-hi
      if (iter.gt.maxit.or.(iter.gt.1.and.abs(gi-g0)/gi.lt.thr)) iexit=1
      if ((iexit.le.0).and.(ipt.lt.1)) go to 330
      k=j5-1
      l=k1-1
      m=k2-1
      ggi=gi*0.25
      write (itape,730) ggi
      do 320 kk=1,2
      t=0.0
      do 300 j=1,p
      s=0.0
      do 290 i=1,nj
      s=s+dp(p*(i-1)+j+l)*a(i,kk,im)/dp(m+i)
 290  continue
      sp(k+j)=s
      t=t+s**2
 300  continue
      t=1./dsqrt(t)
      do 310 j=1,p
      sp(k+j)=sp(j+k)*t
 310  continue
      write (itape,740) kk,(sp(k+j),j=1,p)
 320  continue
 330  if (iexit.gt.0) go to 640
      if (ipt.lt.2) go to 360
      write (6,340) (dp(i+k3-1),i=1,nj)
 340  format (    9h der(1) =5(/    1h 5g12.4))
      write (6,350) (dp(i+k3-1+nj),i=1,nj)
 350  format (    9h der(2) =5(/    1h 5g12.4))
 360  g0=gi
      k=k3-1
      l=j5-1
      nj2=2*nj
      do 370 i=1,nj2
      dp(i+k)=dp(i+k)*dp(i+l4)
 370  continue
      if (ipt.lt.2) go to 400
      write (6,380) (dp(i+k),i=1,nj)
 380  format (    9h stp(1) =5(/    1h 5g12.4))
      write (6,390) (dp(i+k+nj),i=1,nj)
 390  format (    9h stp(2) =5(/    1h 5g12.4))
 400  s=0.d0
      t=s
      u=t
      do 410 i=1,nj
      s=s+a(i,1,im)*dp(k+i)
      t=t+dp(k+i)**2
      u=u+a(i,2,im)*dp(k+i)
 410  continue
      r=r0
      ni=0
 420  ni=ni+1
      fn=1./dsqrt(1.+2.*r*s+(t-u**2)*r**2)
      do 430 i=1,nj
      sp(l+i)=fn*(a(i,1,im)+r*(dp(k+i)-u*a(i,2,im)))
 430  continue
      v=0.0
      do 440 i=1,nj
      v=v+sp(l+i)*dp(k+i+nj)
 440  continue
      do 450 i=1,nj
      sp(l+i+p)=a(i,2,im)+r*(dp(k+i+nj)-v*sp(l+i))
 450  continue
      v=0.0
      do 460 i=1,nj
      v=v+sp(l+i+p)**2
 460  continue
      v=1./dsqrt(v)
      do 470 i=1,nj
      sp(l+i+p)=sp(l+i+p)*v
 470  continue
      call legndr (jj,p,nj,n,sp(j2),w,sp(j5),0,sw,xpa(1,1,im),sp(j7),hj,
     1dp(k3),sp(j3),sp(j4),sp(j8))
      if (ipt.lt.2) go to 490
      write (6,480) r,hj
 480  format (    7h r,hj =2g12.4)
 490  if (ni.ne.1) go to 500
      gink=0.5
      if (hj.lt.hi) gink=1.5
      hg=hi
      rm=0.
 500  if ((r.ge.0.001).and.(r.le.60.0)) go to 510
      iex=1
      go to 550
 510  if (gink.ge.1.0) go to 520
      if (hj.ge.hg) go to 540
      iex=2
      go to 550
 520  if (hj.lt.hi) go to 530
      iex=3
      go to 550
 530  if (ni.le.1) go to 540
      hg=hi
      rm=r0
 540  hi=hj
      r0=r
      r=r*gink
      go to 420
 550  if (iex.ne.1) go to 560
      r0=r
      go to 580
 560  if (iex.ne.2) go to 570
      v=hi
      hi=hj
      hj=v
      v=r0
      r0=r
      r=v
 570  hg=(r-r0)*hg
      hi=(rm-r)*hi
      hj=(r0-rm)*hj
      r0=0.5*((r+r0)*hg+(rm+r0)*hj+(rm+r)*hi)/(hg+hj+hi)
 580  fn=1./dsqrt(1.+2.*r0*s+(t-u**2)*r0**2)
      do 590 i=1,nj
      a(i,1,im)=fn*(a(i,1,im)+r0*(dp(k+i)-u*a(i,2,im)))
 590  continue
      v=0.0
      do 600 i=1,nj
      v=v+a(i,1,im)*dp(k+i+nj)
 600  continue
      do 610 i=1,nj
      a(i,2,im)=a(i,2,im)+r0*(dp(k+i+nj)-v*a(i,1,im))
 610  continue
      v=0.0
      do 620 i=1,nj
      v=v+a(i,2,im)**2
 620  continue
      v=1./dsqrt(v)
      do 630 i=1,nj
      a(i,2,im)=a(i,2,im)*v
 630  continue
      go to 260
 640  l3=j3-1
      l2=j2-1
      l6=j6-1
      l7=j7-1
      do 650 j=1,n
      sp(l6+j)=xpa(1,j,im)
      sp(l7+j)=xpa(2,j,im)
 650  continue
      call binorm (n,sp(j6),sp(j7),w,sw,sp(j3))
      call gs2 (p,nj,a(1,1,im),dp(k5))
      l3=k3-1
      l5=k5-1
      do 700 l=1,n
      do 680 j=1,nj
      s=0.0
      do 670 i=3,nj
      t=0.0
      do 660 k=1,nj
      t=t+dp(l5+p*(i-1)+k)*sp(l2+p*(l-1)+k)
 660  continue
      s=s+dp(l5+p*(i-1)+j)*t
 670  continue
      dp(l3+j)=dp(l5+j)*sp(l6+l)+dp(l5+j+p)*sp(l7+l)+s
 680  continue
      do 690 j=1,nj
      sp(l2+p*(l-1)+j)=dp(l3+j)
 690  continue
 700  continue
 710  continue
      call xform (mm,p,n,nj,a,xpa,sp(j1),dp(k2),dp(k1),sp(j5))
      return
 720  format(    1h0i4,   16h obs deleted in i2,      9h passes. i2,   1
     18h dim search space.)
 730  format(   19h0projection index =g12.4/)
 740  format(     3h a(i1,    3h) =5(/    1h 5g12.4))
 750  format(   15h0pursuit plane i2/)
      end
      subroutine legndr (jj,pp,p,n,z,w,a,ifl,sw,y,x,hi,d,h,eh,ej)
      integer pp,p
      real z(pp,n),w(n),y(2,n),x(2,n),a(1),h(jj,2,n),eh(jj,2),ej(jj,jj)
      double precision r,s,t,d(1)
      data s2i,const /0.707107,1.59577/
      if (ifl.lt.0) go to 40
      do 30 j=1,n
      s=0.0
      do 10 i=1,p
      s=s+a(i)*z(i,j)
 10   continue
      y(1,j)=s
      s=0.0
      do 20 i=1,p
      s=s+a(i+pp)*z(i,j)
 20   continue
      y(2,j)=s
 30   continue
 40   do 70 j=1,n
      x(1,j)=0.0
      x(2,j)=x(1,j)
      if (w(j).le.0.0) go to 70
      do 60 i=1,2
      if (y(i,j).lt.0.0) go to 50
      x(i,j)=erf(y(i,j)*s2i)
      go to 60
 50   x(i,j)=erfc(-y(i,j)*s2i)-1.0
 60   continue
 70   continue
      do 100 i=1,n
      do 90 k=1,2
      u=x(k,i)
      h(1,k,i)=u
      h(2,k,i)=0.5*(3.0*u**2-1.0)
      do 80 j=3,jj
      h(j,k,i)=((2*j-1)*u*h(j-1,k,i)-(j-1)*h(j-2,k,i))/j
 80   continue
 90   continue
 100  continue
      hi=0.0
      do 130 k=1,2
      do 120 j=1,jj
      s=0.0
      do 110 i=1,n
      s=s+w(i)*h(j,k,i)
 110  continue
      eh(j,k)=s/sw
      hi=hi-(2*j+1)*eh(j,k)**2
 120  continue
 130  continue
      do 160 j=1,jj
      jjmj=jj-j
      do 150 k=1,jjmj
      s=0.0
      do 140 i=1,n
      s=s+w(i)*h(j,1,i)*h(k,2,i)
 140  continue
      ej(j,k)=s/sw
      hi=hi-(2*j+1)*(2*k+1)*ej(j,k)**2
 150  continue
 160  continue
      if (ifl.le.0) return
      l1=1
      l2=2*p
      do 260 ll=l1,l2
      l=mod(ll,p)
      if (l.eq.0) l=p
      if (ll.gt.p) go to 170
      m=1
      mo=2
      go to 180
 170  m=2
      mo=1
 180  r=0.0
      do 250 i=1,n
      s=0.0
      do 240 j=1,jj
      t=0.0
      jjmj=jj-j
      do 210 k=1,jj
      if (m.ne.1) go to 190
      if (k.gt.jjmj) go to 210
      ek=ej(j,k)
      go to 200
 190  if (j.gt.jj-k) go to 210
      ek=ej(k,j)
 200  t=t+(2*k+1)*ek*h(k,mo,i)
 210  continue
      u=x(m,i)
      if (j.ne.1) go to 220
      dm=1.0
      dr=dm
      go to 230
 220  dr=u*dm+j*h(j-1,m,i)
      dm=dr
 230  s=s+(2*j+1)*dr*(eh(j,m)+t)
 240  continue
      if (w(i).gt.0.0) r=r+w(i)*exp(-0.5*y(m,i)**2)*s*(z(l,i)-y(1,i)*a(l
     1)-y(2,i)*a(l+pp))
 250  continue
      d(ll)=const*r/sw
 260  continue
      return
      end
      subroutine binorm (n,x,y,w,sw,sc)
      real x(n),y(n),w(n),sc(n,4),sine(4),csin(4)
      data sine /0.0,0.707107,0.382683,0.92388/
      data csin /1.0,0.707107,0.92388,0.382683/
      common /debug/ ipt
      data maxit /20/
      nit=0
      sm=9.9e30
 10   s=0.0
      nit=nit+1
      do 80 ia=1,4
      do 20 j=1,n
      sc(j,1)=x(j)*csin(ia)+y(j)*sine(ia)
      sc(j,2)=sc(j,1)
 20   continue
      do 30 j=1,n
      sc(j,3)=j+.1
 30   continue
      call sort (sc,sc(1,3),1,n)
      call normin (n,w,sw,sc(1,3),sc)
      t=0.
      do 40 j=1,n
      t=t+w(j)*abs(sc(j,1)-sc(j,2))
 40   continue
      s=amax1(s,t/sw)
      do 50 j=1,n
      sc(j,2)=y(j)*csin(ia)-x(j)*sine(ia)
      sc(j,3)=sc(j,2)
 50   continue
      do 60 j=1,n
      sc(j,4)=j+.1
 60   continue
      call sort (sc(1,2),sc(1,4),1,n)
      call normin (n,w,sw,sc(1,4),sc(1,2))
      t=0.0
      do 70 j=1,n
      t=t+w(j)*abs(sc(j,2)-sc(j,3))
      x(j)=sc(j,1)*csin(ia)-sc(j,2)*sine(ia)
      y(j)=sc(j,1)*sine(ia)+sc(j,2)*csin(ia)
 70   continue
      s=amax1(s,t/sw)
 80   continue
      if (s.ge.sm) go to 90
      sm=s
      more=0
      go to 100
 90   if (nit.le.6) go to 100
      more=more+1
      if (more.ge.4) go to 130
 100  if (ipt.lt.3) go to 120
      write (6,110) nit,s,sm
 110  format (    7h binorm i4,2g12.4)
 120  if (nit.lt.maxit) go to 10
 130  return
      end
      subroutine normin (n,w,sw,x,y)
      real x(n),y(n),w(n)
      l=1
 10   if (w(int(x(l))).gt.0.0) go to 20
      l=l+1
      go to 10
 20   k=n
 30   if (w(int(x(k))).gt.0.0) go to 40
      k=k-1
      go to 30
 40   s=-0.5*amin1(w(int(x(l))),w(int(x(k))))/sw
      do 50 j=l,k
      s=s+w(int(x(j)))/sw
      y(int(x(j)))=phinv(s)
 50   continue
      if (l.le.1) go to 70
      do 60 j=1,l
      y(int(x(j)))=y(int(x(l)))
 60   continue
 70   if (k.ge.n) go to 90
      do 80 j=k,n
      y(int(x(j)))=y(int(x(k)))
 80   continue
 90   return
      end
      function phinv (s)
      double precision p(5),q(5),u
      data p /-0.322232431088d0,-1.d0,-0.342242088547d0,  -0.02042312102
     145d0,-0.453642210148d-4/
      data q /0.0993484626060d0,0.588581570495d0,0.531103462366d0,  0.10
     13537752850d0,0.38560700634d-2/
      u=s
      if (s.ne.0.5) go to 10
      phinv=0.0
      return
 10   if (s.le.0.5) go to 20
      u=1.0-u
 20   u=dsqrt(-2.0*dlog(u))
      u=u+((((u*p(5)+p(4))*u+p(3))*u+p(2))*u+p(1))/((((u*q(5)+q(4))*u+q(
     13))*u+q(2))*u+q(1))
      phinv=u
      if (s.lt.0.5) phinv=-u
      return
      end
      subroutine gs2 (pp,p,a,u)
      integer pp,p
      real a(pp,2)
      double precision u(pp,pp),s,t
      do 10 i=1,p
      u(i,1)=a(i,1)
      u(i,2)=a(i,2)
 10   continue
      do 50 k=3,p
      km1=k-1
      t=0.0
      do 30 j=1,p
      s=0.0
      if (j.eq.k) s=1.0
      do 20 l=1,km1
      s=s-u(k,l)*u(j,l)
 20   continue
      t=t+s**2
      u(j,k)=s
 30   continue
      t=1.0/dsqrt(t)
      do 40 j=1,p
      u(j,k)=u(j,k)*t
 40   continue
 50   continue
      return
      end
      subroutine xform (mm,p,n,nj,a,xpa,xm,eiv,eit,scr)
      integer p
      real a(p,2,mm),xpa(2,n,mm),xm(p),scr(p)
      double precision eiv(p),eit(p,p),s,t,u
      do 60 im=1,mm
      do 50 kk=1,2
      t=0.0
      u=t
      do 20 j=1,p
      s=0.0
      do 10 i=1,nj
      s=s+eit(j,i)*a(i,kk,im)/eiv(i)
 10   continue
      scr(j)=s
      t=t+s**2
      u=u+s*xm(j)
 20   continue
      t=1.0/dsqrt(t)
      do 30 j=1,p
      a(j,kk,im)=scr(j)*t
 30   continue
      do 40 j=1,n
      xpa(kk,j,im)=t*(xpa(kk,j,im)+u)
 40   continue
 50   continue
 60   continue
      return
      end
      block data
      common /debug/ ipt
      data ipt /0/
      end
      subroutine sphere (p,n,x,w,sw,efact,nei,itape,xm,eiv,eit,nj,z,scr)
      integer p
      real x(p,n),w(n),z(p,n),xm(p)
      double precision s,machep,tol,eiv(p),eit(p,p),scr(p)
      machep=0.5**(52)
      tol=1.0e-60
      sw=0.0
      do 10 j=1,n
      sw=sw+w(j)
 10   continue
      if (sw.gt.0.0) go to 20
      write (itape,180) sw
      stop
 20   do 40 i=1,p
      s=0.0
      do 30 j=1,n
      s=s+w(j)*x(i,j)
 30   continue
      xm(i)=s/sw
 40   continue
      do 60 j=1,n
      do 50 i=1,p
      z(i,j)=x(i,j)-xm(i)
 50   continue
 60   continue
      do 90 i=1,p
      do 80 j=1,i
      s=0.0
      do 70 k=1,n
      s=s+w(k)*z(i,k)*z(j,k)
 70   continue
      eit(i,j)=s/sw
 80   continue
 90   continue
      call tred2 (p,p,p,tol,eit,eiv,scr,eit)
      call imtql2 (p,p,p,machep,eiv,scr,eit,ierr)
      if (ierr.eq.0) go to 100
      write (itape,190) ierr
      stop
 100  if (eiv(1).gt.0.0) go to 110
      write (itape,200) eiv(1)
      stop
 110  k=1
 120  if (eiv(k).le.0.0) go to 130
      eiv(k)=dsqrt(eiv(k))
      k=k+1
      if (k.le.p) go to 120
 130  nj=1
      do 140 i=2,p
      if (eiv(i).gt.efact*eiv(1)) nj=nj+1
 140  continue
      nj=min0(nj,nei)
      do 170 j=1,n
      do 160 k=1,nj
      s=0.0
      do 150 i=1,p
      s=s+eit(i,k)*(x(i,j)-xm(i))
 150  continue
      z(k,j)=s/eiv(k)
 160  continue
 170  continue
      return
 180  format(   17h0sum of weights =g12.4,   18h  is not positive.)
 190  format(   33h0matrix error in sphereing. err =i5)
 200  format(   21h0largest eigenvalue =g12.4,   18h  is not positive.)
      end
      subroutine imtql2 (nt,nm,n,machep,d,e,z,error)
      implicit real*8 (a-h,o-z)
      real*8 machep,d(n),e(n),z(nt,nt)
      integer error
      error=0
      if (n.eq.1) go to 140
      do 10 i=2,n
 10   e(i-1)=e(i)
      e(n)=0.0d0
      do 90 l=1,n
      j=0
 20   do 30 m=l,n
      if (m.eq.n) go to 40
      if (dabs(e(m)).le.machep*(dabs(d(m))+dabs(d(m+1)))) go to 40
 30   continue
 40   p=d(l)
      if (m.eq.l) go to 90
      if (j.eq.30) go to 130
      j=j+1
      g=(d(l+1)-p)/(2.0d0*e(l))
      r=dsqrt(1.0d0+g*g)
      g=d(m)-p+e(l)/(g+dsign(r,g))
      s=1.0d0
      c=1.0d0
      p=d(m)
      mml=m-l
      do 80 ii=1,mml
      i=m-ii
      f=s*e(i)
      b=c*e(i)
      if (dabs(f).lt.dabs(g)) go to 50
      c=g/f
      r=dsqrt(c*c+1.0d0)
      e(i+1)=f*r
      s=1.0d0/r
      c=c/r
      go to 60
 50   c=f/g
      r=dsqrt(c*c+1.0d0)
      e(i+1)=g*r
      s=c/r
      c=1.0d0/r
 60   f=c*d(i)-s*b
      g=c*b-s*p
      r=d(i)+p
      p=c*f-s*g
      g=s*f+c*g
      d(i+1)=r-p
      do 70 ia=1,n
      f=z(ia,i+1)
      z(ia,i+1)=s*z(ia,i)+c*f
      z(ia,i)=c*z(ia,i)-s*f
 70   continue
 80   continue
      d(l)=p
      e(l)=g
      e(m)=0.0d0
      go to 20
 90   continue
      nm1=n-1
      do 120 i=1,nm1
      k=i
      p=d(i)
      ip1=i+1
      do 100 j=ip1,n
      if (d(j).le.p) go to 100
      k=j
      p=d(j)
 100  continue
      if (k.eq.i) go to 120
      d(k)=d(i)
      d(i)=p
      do 110 j=1,n
      p=z(j,i)
      z(j,i)=z(j,k)
      z(j,k)=p
 110  continue
 120  continue
      go to 140
 130  error=l
 140  return
      end
      subroutine tred2 (nt,nm,n,tol,a,d,e,z)
      implicit real*8  (a-h,o-z)
      real*8 a(nt,nt),d(n),e(n),z(nt,nt)
      do 10 i=1,n
      do 10 j=1,i
      z(i,j)=a(i,j)
 10   continue
      if (n.eq.1) go to 120
      do 110 ii=2,n
      i=n+2-ii
      l=i-2
      f=z(i,i-1)
      g=0.0d0
      if (l.lt.1) go to 30
      do 20 k=1,l
 20   g=g+z(i,k)*z(i,k)
 30   h=g+f*f
      if (g.gt.tol) go to 40
      e(i)=f
      h=0.0d0
      go to 100
 40   l=l+1
      g=-dsign(dsqrt(h),f)
      e(i)=g
      h=h-f*g
      z(i,i-1)=f-g
      f=0.0d0
      do 80 j=1,l
      z(j,i)=z(i,j)/h
      g=0.0d0
      do 50 k=1,j
 50   g=g+z(j,k)*z(i,k)
      jp1=j+1
      if (l.lt.jp1) go to 70
      do 60 k=jp1,l
 60   g=g+z(k,j)*z(i,k)
 70   e(j)=g/h
      f=f+g*z(j,i)
 80   continue
      hh=f/(h+h)
      do 90 j=1,l
      f=z(i,j)
      g=e(j)-hh*f
      e(j)=g
      do 90 k=1,j
      z(j,k)=z(j,k)-f*e(k)-g*z(i,k)
 90   continue
 100  d(i)=h
 110  continue
 120  d(1)=0.0d0
      e(1)=0.0d0
      do 170 i=1,n
      l=i-1
      if (d(i).eq.0.0d0) go to 150
      do 140 j=1,l
      g=0.0d0
      do 130 k=1,l
 130  g=g+z(i,k)*z(k,j)
      do 140 k=1,l
      z(k,j)=z(k,j)-g*z(k,i)
 140  continue
 150  d(i)=z(i,i)
      z(i,i)=1.0d0
      if (l.lt.1) go to 170
      do 160 j=1,l
      z(i,j)=0.0d0
      z(j,i)=0.0d0
 160  continue
 170  continue
      return
      end
      subroutine sort (v,a,ii,jj)
c
c     puts into a the permutation vector which sorts v into
c     increasing order.  only elements from ii to jj are considered.
c     arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements
c
c     this is a modification of cacm algorithm #347 by r. c. singleton,
c     which is a modified hoare quicksort.
c
      dimension a(jj),v(1),iu(20),il(20)
      integer t,tt
      integer a
      real v
      m=1
      i=ii
      j=jj
 10   if (i.ge.j) go to 80
 20   k=i
      ij=(j+i)/2
      t=a(ij)
      vt=v(ij)
      if (v(i).le.vt) go to 30
      a(ij)=a(i)
      a(i)=t
      t=a(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
 30   l=j
      if (v(j).ge.vt) go to 50
      a(ij)=a(j)
      a(j)=t
      t=a(ij)
      v(ij)=v(j)
      v(j)=vt
      vt=v(ij)
      if (v(i).le.vt) go to 50
      a(ij)=a(i)
      a(i)=t
      t=a(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
      go to 50
 40   a(l)=a(k)
      a(k)=tt
      v(l)=v(k)
      v(k)=vtt
 50   l=l-1
      if (v(l).gt.vt) go to 50
      tt=a(l)
      vtt=v(l)
 60   k=k+1
      if (v(k).lt.vt) go to 60
      if (k.le.l) go to 40
      if (l-i.le.j-k) go to 70
      il(m)=i
      iu(m)=l
      i=k
      m=m+1
      go to 90
 70   il(m)=k
      iu(m)=j
      j=l
      m=m+1
      go to 90
 80   m=m-1
      if (m.eq.0) return
      i=il(m)
      j=iu(m)
 90   if (j-i.gt.10) go to 20
      if (i.eq.ii) go to 10
      i=i-1
 100  i=i+1
      if (i.eq.j) go to 80
      t=a(i+1)
      vt=v(i+1)
      if (v(i).le.vt) go to 100
      k=i
 110  a(k+1)=a(k)
      v(k+1)=v(k)
      k=k-1
      if (vt.lt.v(k)) go to 110
      a(k+1)=t
      v(k+1)=vt
      go to 100
      end
      subroutine varmax (mm,p,n,w,a,xpa,itape)
c
c performes varimax rotations for two - dimensional pp solutions.
c
c the linear combinations and associated adjusted data plots are
c rotated so that the variance of the original (unsphered) variable
c loadings on the vertical solution coordinate (ordinate) is
c maximized. this sometimes helps in interpreting the solutions.
c
c input:
c
c mm, p, n, w, itape = same as in ppdeux (see above).
c a, xpa = output from ppdeux (see above).
c
c output:
c
c a, xpa = corresponding varimax rotated quanities (replaces input).
c
      integer p
      real w(n),a(p,2,mm),xpa(2,n,mm),xm(2),an(2)
      double precision s,t,sw
      sw=0.0
      do 10 j=1,n
      sw=sw+w(j)
 10   continue
      if (itape.gt.0) write (itape,130)
      do 120 im=1,mm
      do 60 k=1,2
      s=0.0
      do 20 j=1,n
      s=s+w(j)*xpa(k,j,im)
 20   continue
      s=s/sw
      xm(k)=s
      t=0.0
      do 30 j=1,n
      t=t+w(j)*(xpa(k,j,im)-s)**2
 30   continue
      t=1./dsqrt(t/sw)
      xm(k)=xm(k)*t
      do 40 i=1,p
      a(i,k,im)=a(i,k,im)*t
 40   continue
      do 50 j=1,n
      xpa(k,j,im)=(xpa(k,j,im)-s)*t
 50   continue
 60   continue
      call rotmin (p,a(1,1,im),a(1,2,im),c,sc)
      do 90 k=1,2
      s=0.0
      do 70 i=1,p
      s=s+a(i,k,im)**2
 70   continue
      s=1./dsqrt(s)
      an(k)=s
      do 80 i=1,p
      a(i,k,im)=s*a(i,k,im)
 80   continue
 90   continue
      x=xm(1)
      y=xm(2)
      xm(1)=c*x+sc*y
      xm(2)=c*y-sc*x
      do 100 j=1,n
      x=xpa(1,j,im)
      y=xpa(2,j,im)
      xpa(1,j,im)=an(1)*(c*x+sc*y+xm(1))
      xpa(2,j,im)=an(2)*(c*y-sc*x+xm(2))
 100  continue
      if (itape.le.0) go to 120
      write (itape,140) im
      do 110 k=1,2
      write (itape,150) k,(a(i,k,im),i=1,p)
 110  continue
 120  continue
      return
 130  format(   27h0varimax rotated solutions:)
 140  format(   15h0pursuit plane i2,    1h:)
 150  format(    3h a(i1,    3h) = 5(/    1h 5g12.4))
      end
      subroutine rotmin (p,a,b,c,s)
      integer p
      real a(p),b(p)
      data omt,big /0.381966,9.9e30/
      b1=1.0
      a1=-1.0
      d1=(b1-a1)/100.0
      fd=big
      c1=a1
      go to 20
 10   c1=c1+(d1)
 20   if ((d1)*((c1)-(b1)).gt.0) go to 40
      s=sqrt(1.0-c1**2)
      t=0.0
      v=t
      do 30 i=1,p
      bc=c1*b(i)-s*a(i)
      t=t+abs(bc)
      v=v+bc**2
 30   continue
      fc=t**2/v
      if (fc.ge.fd) go to 10
      fd=fc
      c=c1
      go to 10
 40   if (c.le.a1.or.c.ge.b1) go to 120
      a1=c-d1
      b1=c+d1
      c1=a1+omt*(b1-a1)
      d1=b1-omt*(b1-a1)
      s=sqrt(1.0-c1**2)
      t=0.0
      v=t
      do 50 i=1,p
      bc=c1*b(i)-s*a(i)
      t=t+abs(bc)
      v=v+bc**2
 50   continue
      fc=t**2/v
      s=sqrt(1.0-d1**2)
      t=0.0
      v=t
      do 60 i=1,p
      bc=d1*b(i)-s*a(i)
      t=t+abs(bc)
      v=v+bc**2
 60   continue
      fd=t**2/v
 70   if (fc.ge.fd) go to 90
      b1=d1
      d1=c1
      fd=fc
      c1=a1+omt*(b1-a1)
      s=sqrt(1.0-c1**2)
      t=0.0
      v=t
      do 80 i=1,p
      bc=c1*b(i)-s*a(i)
      t=t+abs(bc)
      v=v+bc**2
 80   continue
      fc=t**2/v
      go to 110
 90   a1=c1
      c1=d1
      fc=fd
      d1=b1-omt*(b1-a1)
      s=sqrt(1.0-d1**2)
      t=0.0
      v=t
      do 100 i=1,p
      bc=d1*b(i)-s*a(i)
      t=t+abs(bc)
      v=v+bc**2
 100  continue
      fd=t**2/v
 110  if (b1-a1.ge.0.0001) go to 70
      c=0.5*(a1+b1)
 120  s=sqrt(1.0-c**2)
      do 130 i=1,p
      a1=a(i)
      b1=b(i)
      a(i)=c*a1+s*b1
      b(i)=c*b1-s*a1
 130  continue
      return
      end
      subroutine nscore (p,n,w,x,p1,p2,sc)
c
c transforms variable marginal distributions to standard normality.
c variables p1 through p2 (inclusive) are transformed.
c
c input:
c
c p,n,w,x = same as in ppdeux (see above).
c p1, p2 = range of variables to be transformed (integer).
c
c output:
c
c x(p,n) = data matrix with columns p1 through p2 replaced with
c          transformed quanities. other columns are unchanged.
c
c workspace:
c
c sc(n,3) : single precision.
c
c notes:
c
c this routine can be called repeatedly to transform any set of
c variables desired.
c
c ranums(z,l) generates l real random numbers from a uniform
c distribution in the range (0.0,1.0) and stores them in the array
c z(l).
c
      integer p,p1,p2
      real w(n),x(p,n),sc(n,3)
      sw=0.0
      do 10 i=1,n
      sw=sw+w(i)
 10   continue
      do 80 i=p1,p2
      do 20 j=1,n
      sc(j,1)=x(i,j)
      sc(j,2)=j+0.1
 20   continue
      call sort (sc,sc(1,2),1,n)
      j=1
 30   j0=j
      if (j.ge.n) go to 50
 40   if (sc(j+1,1).gt.sc(j,1)) go to 50
      j=j+1
      if (j.lt.n) go to 40
 50   if (j.le.j0) go to 60
      call ranums (sc(j0,3),j-j0+1)
      call sort (sc(1,3),sc(1,2),j0,j)
 60   j=j+1
      if (j.le.n) go to 30
      call normin (n,w,sw,sc(1,2),sc)
      do 70 j=1,n
      x(i,j)=sc(j,1)
 70   continue
 80   continue
      return
      end
      subroutine atoscl (p,n,w,x,z)
c
c centers and scales variables to have zero mean and unit variance.
c
c although not necessary, this can aid in interpreting the linear
c combinations associated with pp solutions. it can also help
c numerical stability in extreme circumstances.
c
c input:
c p, n, w, x = same as in ppdeux (see above).
c
c output:
c
c z(p,n) = corresponding centered and scaled data matrix (to be
c          input to ppdeux).
c
c note:
c
c x and z can be the same array in the call to atoscl.
c
      integer p
      real w(n),x(p,n),z(p,n)
      sw=0.0
      do 10 j=1,n
      sw=sw+w(j)
 10   continue
      do 50 i=1,p
      s=0.0
      t=s
      do 20 j=1,n
      s=s+w(j)*x(i,j)
 20   continue
      s=s/sw
      do 30 j=1,n
      z(i,j)=x(i,j)-s
      t=t+w(j)*z(i,j)**2
 30   continue
      t=1.0/sqrt(t/sw)
      do 40 j=1,n
      z(i,j)=t*z(i,j)
 40   continue
 50   continue
      return
      end
      subroutine eivs (p,dp,itape)
c
c prints out eigenvalues of the data covariance matrix.
c
c input:
c
c p, itape = same as in ppdeux (see above).
c dp = double precision workspace array from ppdeux (see above).
c
      integer p
      double precision dp(1),s
      k=p**2
      do 10 i=1,p
      dp(i+k)=dp(i+k)**2
 10   continue
      write (itape,40)
      write (itape,50) (dp(i+k),i=1,p)
      s=0.0
      do 20 i=1,p
      s=s+dp(i+k)
 20   continue
      dp(k+1)=dp(k+1)/s
      do 30 i=2,p
      dp(i+k)=dp(i+k)/s+dp(i+k-1)
 30   continue
      write (itape,60)
      write (itape,50) (dp(i+k),i=1,p)
      return
 40   format(   31h0covariance matrix eigenvalues:)
 50   format(/(    1h 5g12.4))
 60   format(   21h0cumulative fraction:)
      end
      subroutine seed (iseed)
      real x(1)
      double precision u
      save i
      data i /987654321/
      i=iseed
      return
      entry ranums (x,n)
      do 1 j=1,n
      i=dmod(i*16807.d0,2147483647.d0)
      u=i
      u=u*.465661287d-9
      x(j)=u
    1 continue
      return
      end