	 program SFPDv1
	 
	 
	
c	copyright 2002 by Leo Breiman

c       This program is free software; you can redistribute it and/or
c       modify it under the terms of the GNU General Public License
c       as published by the Free Software Foundation; either version 2
c       of the License, or (at your option) any later version.

c       This program is distributed in the hope that it will be useful,
c       but WITHOUT ANY WARRANTY; without even the implied warranty of
c       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c       GNU General Public License for more details.

	

        implicit double precision (a-h,o-z)

	parameter(ns=136,mdim=6,jbt=100,
     1	ndsize=1,nrnodes=2*ns+1,ntsm=100,
     1	look=10,itime=0,isurv=0,indsurv=0,
     1	icor=1,ireg=0,
     1	iscale=0,mdimsc=3)
    
     	double precision td(ns),tm(nrnodes),ut(ns),xt(ns),
     1	tmin(nrnodes),tmax(nrnodes),xsplit(nrnodes),v(ns),
     1	xb(mdim,ns),tb(ns),bestcrit(nrnodes),sgn(mdim),
     1	eml(ns),trtp(ns),tp(ns),tx(ns),avb(ntsm,mdim),
     1	erg(ntsm),vt(ntsm),gtime(ntsm,ns),sot(mdim,ns),
     1	zt(mdim),rsnodecost(nrnodes),avt(ntsm),ztt(ns),
     1	epred(mdim),ftime(ntsm,ns),otd(ns),avrsq(ntsm),
     1	agtime(ntsm,ns),gtime1(ntsm,ns),wtt(ns),ytt(ns),
     1	aftime(ntsm,ns),azip(ns,ntsm),gin(ns),xc(mdim),
     1	b(mdim),cov(mdim,mdim),y(ns),xy(mdim),xtt(ns),
     1	x(mdim,ns),acor(ntsm,mdim),yt(ns),ty(ns),atb(mdim),
     1	avcor(mdim)
     	
c	used in scaling 
     	double precision prox(ns,ns),zy(ns),u(ns),dl(mdimsc),
     1	xsc(ns,mdimsc),red(ns)
     
	integer ce(ns),nd(nrnodes),jdex(ns),ncase(ns),isg(mdim),
     1	nodestatus(nrnodes),nodepop(nrnodes),msplit(nrnodes),
     1	nodestart(nrnodes),treemap(2,nrnodes),parent(nrnodes),
     1	nalive(nrnodes),nout(ns),cb(ns),jin(ns),nvec(nrnodes),
     1	jbig(ns,nrnodes),treeseq(nrnodes,nrnodes),
     1	cat(0:mdim),noc(nrnodes,ns),ncx(ns),ivncx(ns),
     1	nox(mdim,ns),ivnox(mdim,ns)
     
     	if(isurv.eq.1) open(7,file='survout2',status='new')
	if(indsurv.eq.1) open(8,file='indsurvout',status='new')
	if(icor.ge.1) open(9,file='corrout',status='new') 
	if(ireg.ge.1) open(1,file='regout',status='new') 
	if(iscale.eq.1) open(2,file='scaleout',status='new') 
     	
	open(9, file='vet-lung',status='old')
	do n=1,ns
	read(9,*) x(1,n),x(2,n),td(n),ce(n),(x(m,n),m=3,6)
	if(dnint(x(6,n)).eq.0) x(6,n)=1
	if(dnint(x(6,n)).eq.10) x(6,n)=2
	end do
	cat(1)=2
	cat(2)=4
	cat(3)=1
	cat(4)=1
	cat(5)=1
	cat(6)=2
	close(9)
	
	
	call sortdie(x,td,ce,xt,ncase,mdim,ns)
	call reducetime(ce,td,otd,vt,ntsm,ns,ncase)

	

	do n=1,271
	zzt=rand(1)
	end do
	
	call izerv(nout,ns)
	call zermr(gtime,ntsm,ns)
	call zermr(gtime1,ntsm,ns)
	if(iscale.eq.1) call zermr(prox,ns,ns)
	call zervr(avrsq,ntsm)
	call zermr(avb,ntsm,mdim)
	
	
	
c	start iteration $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
	
	do jb=1,jbt
	
	call izerv(jin,ns)
	do n=1,ns
	k=int(rand(1)*ns)+1
	do m=1,mdim
	xb(m,n)=x(m,k)
	end do
	tb(n)=td(k)
	cb(n)=ce(k)
	jin(k)=1
	end do
	
	do n=1,ns
	if(jin(n).eq.0) nout(n)=nout(n)+1
	end do
	
	call buildtree(xb,tm,tb,ut,xt,tmin,tmax,xsplit,
     1	cb,nd,jdex,ncase,nodestatus,nodepop,msplit,nodestart,
     1	v,treemap,parent,ns,mdim,nrnodes,bestcrit,
     1	ndsize,jbig,cat,rsnodecost,itime,ndbigtree)
  

     	call where(x,xsplit,nalive,msplit,treemap,
     1	nodestatus,ns,mdim,nrnodes,ndbigtree,cat,noc)
     
	
	if(iscale.eq.1) then
        call proxim(prox,nvec,noc,ns,nrnodes)
        end if
	
	do k=1,ntsm
	tip=vt(k)
	call surv(eml,nodestatus,jin,noc,tmax,tmin,
     1	nd,tm,td,ndbigtree,ns,nrnodes,tip)
	do n=1,ns
	if(jin(n).eq.0) then
	gtime(k,n)=gtime(k,n)+exp(-eml(n))
	gtime1(k,n)=gtime1(k,n)-eml(n)
	agtime(k,n)=-gtime1(k,n)/nout(n)
	end if
	end do
	end do !k
	
	if(mod(jb,look).eq.0) then
	call errorpred(agtime,vt,ns,ntsm,azip,gin,td,
     &	ce,xt,ut,std)
     	write(*,'(i5,4f10.3)') jb,100*(1-std)
	end if
	
	end do !jbt
	
c	end iteration $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
	
	do k=1,ntsm
	avt(k)=0
	pro=(1-1.0*(dble(k)/ntsm))
	if(itime.eq.0) pro=1
	do n=1,ns
	gtime(k,n)=gtime(k,n)/nout(n)
	gtime1(k,n)=exp(gtime1(k,n)/nout(n))
	gtime(k,n)=pro*gtime(k,n)+(1-pro)*gtime1(k,n)
	agtime(k,n)=-dlog(gtime(k,n))
	avt(k)=avt(k)+gtime(k,n)
	end do
	avt(k)=avt(k)/ns
	end do
	
	if(isurv.eq.1) then
	do k=1,ntsm
	write(7,'(2f10.3)') vt(k),avt(k)
	end do
	stop
	end if
	
	if(indsurv.eq.1) then
	do n=1,ns
	write(8,'(200f10.3)') (gtime(k,n),k=1,ntsm),ce(n),
     1	td(n),(x(m,n),m=1,mdim)
	end do
	close(8)
	end if
	
	
	
	if(ireg.ge.1) then
	call regress(agtime,vt,ntsm,mdim,ns,cov,b,y,x,sgn,xy,
     1	xt,yt,xtt,wtt,ytt,ztt,ncx,ivncx,sot,nox,ivnox,
     1	cat,ut,ireg,avrsq,avb,atb,trsq)
     	write(1,'(20f10.3)')  100.0, trsq, (atb(m),m=1,mdim)
	write(1,*) ""
     	do k=1,ntsm-1
	write(1,'(20f10.3)') vt(k),avrsq(k),(avb(k,m),m=1,mdim)
	end do
	close(1)
	pause
	end if
	
	if(icor.ge.1) then
	call varnorm(x,mdim,ns)
	call corr(agtime,x,vt,ntsm,mdim,ns,ncx,xt,
     1	yt,acor,ty,ytt,zt,ut,cat,icor,sot,avcor)
     	write(9,'(20f10.3)')100.0, (avcor(m),m=1,mdim)
	write(9,*) ""
     	do k=1,ntsm
	write(9,'(20f10.3)')vt(k),(acor(k,m),m=1,mdim)
	end do
	close(9)
	pause
	end if
	
	if(iscale.eq.1) then
         do n=1,ns
         do k=1,ns
	if(n.ne.k) prox(n,k)=prox(n,k)/(dsqrt(prox(n,n)*prox(k,k)))
         end do
	end do
	do n=1,ns
	prox(n,n)=1
	end do
	end if
	
	if(iscale.eq.1) then
	call scale(prox,xsc,zy,u,dl,ns,mdimsc,red)
	do n=1,mdimsc
	write(*,'(i5,f10.3)') n,dl(n)
	end do
	do n=1,ns
	write(2,'(2i5,15f10.3)') n,ce(n),td(n), 
     1	(xsc(n,k),k=1,mdimsc), (x(m,n),m=1,mdim)
     	end do
	close(2)
	end if
	end 

c	END MAIN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
 	subroutine buildtree(x,tm,td,ut,xt,tmin,tmax,xsplit,
     1	ce,nd,jdex,ncase,nodestatus,nodepop,msplit,nodestart,
     1	v,treemap,parent,ns,mdim,nrnodes,bestcrit,
     1	ndsize,jbig,cat,rsnodecost,itime,ndbigtree)
     	implicit double precision (a-h,o-z)
     
	
	double precision x(mdim,ns),td(ns),tm(nrnodes),th,
     1	tmin(nrnodes),tmax(nrnodes), rsnodecost(nrnodes),
     1	xsplit(nrnodes),v(ns),bestcrit(nrnodes),ut(ns),xt(ns)

     
     
	integer ce(ns),nd(nrnodes),jdex(ns),cat(0:mdim),ncase(ns),
     1	nodestatus(nrnodes),nodepop(nrnodes),msplit(nrnodes),
     1	nodestart(nrnodes),treemap(2,nrnodes),parent(nrnodes),
     1	jbig(ns,nrnodes)
     
     	if(itime.eq.1) then
	th=.75
	else
	th=.5
	end if
     
     	call izerv(nodestatus,nrnodes)
	call izerv(msplit,nrnodes)
	call izerv(nd,nrnodes)
	call zervr(tm,nrnodes)

	nd(1)=0
	tm(1)=0
	parent(1)=0
	
	do n=1,ns
	nd(1)=nd(1)+ce(n)
	tm(1)=tm(1)+td(n)
	end do
	
	
	
	do n=1,ns
	ut(n)=0
	jbig(n,1)=n
	end do
	
	ncur=1
	nodestart(1)=1
	nodepop(1)=ns
	nodestatus(1)=2
	rsnodecost(1)=-nd(1)*dlog(1.0*nd(1)/tm(1))
	tmin(1)=0
	tmax(1)=0
	do n=1,ns
	tmax(1)=dmax1(tmax(1),td(n))
	end do
	range=tmax(1)-tmin(1)
	
c	start main loop
	
	do kbuild=1,nrnodes
	if (kbuild.gt.ncur) goto 50
	if(ncur+2.gt.nrnodes) goto 50
	if (nodestatus(kbuild).ne.2) goto 30
	
	
c		initialize for next call to findbestsplit

	ndstart=nodestart(kbuild)
	ndend=ndstart+nodepop(kbuild)-1
	ndnd=nd(kbuild)
	
	tmnd=tm(kbuild)
	tmint=tmin(kbuild)
	tmaxt=tmax(kbuild)
	
	do n=ndstart,ndend
	jdex(n)=jbig(n,kbuild)
	end do
	jstat=0
	
	call findbestsplit(x,ce,td,jdex,ndnd,tmnd,tmint,
     1	ns,mdim,v,xt,ut,ncase,ubest,mbest,ndrt,ndlt,tmlt,
     1	tmrt,ndendl,ndstart,ndend,jstat,critinc,kbuild,tmaxt,
     1	cat,th,range)
     	
     	
     	if(jstat.eq.1) then
	nodestatus(kbuild)=-1
	goto 30
	end if
	
	msplit(kbuild)=mbest
	
	xsplit(kbuild)=ubest
	
	
	
	nd(ncur+1)=ndlt
	nd(ncur+2)=ndrt
	tm(ncur+1)=tmlt
	tm(ncur+2)=tmrt
	
	
	rsnodecost(ncur+1)=-nd(ncur+1)*
     1	dlog(1.0*nd(ncur+1)/tm(ncur+1))
	rsnodecost(ncur+2)=-nd(ncur+2)*
     1	dlog(1.0*nd(ncur+2)/tm(ncur+2))
	
	if(msplit(kbuild).ge.1) then
	tmin(ncur+1)=tmin(kbuild)
	tmin(ncur+2)=tmin(kbuild)
	tmax(ncur+1)=tmax(kbuild)
	tmax(ncur+2)=tmax(kbuild)
	else
	tmin(ncur+1)=tmin(kbuild)
	tmax(ncur+1)=ubest
	tmin(ncur+2)=ubest
	tmax(ncur+2)=tmax(kbuild)
	end if
	
	
	
	bestcrit(kbuild)=critinc
	
c	leftnode no.= ncur+1, rightnode no. = ncur+2
c	find populations in both nodes.
	
	if(msplit(kbuild).ge.1) then
	nodepop(ncur+1)=ndendl-ndstart+1
	nodepop(ncur+2)=ndend-ndendl
	nodestart(ncur+1)=ndstart
	nodestart(ncur+2)=ndendl+1
	else
	nodepop(ncur+1)=nodepop(kbuild)
	nodepop(ncur+2)=nodepop(kbuild)
	nodestart(ncur+1)=nodestart(kbuild)
	nodestart(ncur+2)=nodestart(kbuild)
	end if
	
	
	do n=ndstart,ndend
	jbig(n,ncur+1)=jdex(n)
	end do
	do n=ndstart,ndend
	jbig(n,ncur+2)=jdex(n)
	end do
	

c	check on nodestatus
	
	nodestatus(ncur+1)=2
	nodestatus(ncur+2)=2
	
	if(ndlt.le.1) nodestatus(ncur+1)=-1
	if(ndrt.le.1) nodestatus(ncur+2)=-1
	if (nodepop(ncur+1).le.ndsize) nodestatus(ncur+1)=-1
	
	treemap(1,kbuild)=ncur+1
	treemap(2,kbuild)=ncur+2
	parent(ncur+1)=kbuild
	parent(ncur+2)=kbuild
	nodestatus(kbuild)=1
	
	ncur=ncur+2
	
	if (ncur.ge.nrnodes) goto 50

	
30	continue
	end do !kbuild
	
50	continue

c	write(*,*) "end"
c	pause
	
	ndbigtree=nrnodes
	do k=nrnodes,1,-1
	if (nodestatus(k).eq.0) ndbigtree=ndbigtree-1
	if (nodestatus(k).eq.2) nodestatus(k)=-1
	end do
	end
	
	subroutine findbestsplit(x,ce,td,jdex,ndnd,tmnd,tmint,
     1	ns,mdim,v,xt,ut,ncase,ubest,mbest,ndrt,ndlt,tmlt,
     1	tmrt,ndendl,ndstart,ndend,jstat,critinc,kbuild,tmaxt,
     1	cat,th,range)
     	implicit double precision (a-h,o-z)
	
	double precision x(mdim,ns),xt(ns),ut(ns),td(ns),v(ns)
    	integer ce(ns),jdex(ns),ncase(ns),cat(0:mdim)
    
     	parameter(nbig=1000,nmod=100,nsm=32)
     
     	double precision  tt(nbig),ttcat(nsm),w(nbig)
	integer idiecat(nsm),icat(nsm),nv(nmod),ct(nbig)
     	

	call izerv(ct,nbig)
	ndnd=0
	totnd=0
	do n=ndstart,ndend
	ct(n)=ce(jdex(n))
	if(td(jdex(n)).le.tmaxt.and.td(jdex(n)).gt.
     1	tmint.and.ct(n).eq.1) then
         w(n)=td(jdex(n))
	ndnd=ndnd+1
	else
	w(n)=0
	ct(n)=0
	end if
	w(n)=w(n)*ct(n)
	tt(n)=xplus(dmin1(tmaxt,td(jdex(n)))-tmint)
	end do
	
	crit0=ndnd*dlog(1.0*ndnd/tmnd)
	critvar=-1.0e10
	
	
	
	vv=rand(1)
	
	
	call izerv(nv,200)

	do m=1,mdim
50	mv=int(rand(1)*mdim)+1
	if(nv(mv).eq.1) goto 50
	if(vv.le.th) mv=0
	if(mv.ge.1) nv(mv)=1
	if(mv.ge.1.and.cat(mv).ge.2) goto 70
	ndl=0
	ndr=ndnd
	tml=0
	tmr=tmnd
	
	

	sl=0
	sr=0
	do n=ndstart,ndend
	if(mv.eq.0) then
	xt(n)=dmin1(tmaxt-tmint,
     1	xplus(td(jdex(n))-tmint))
	sr=sr+xt(n)
	else
	xt(n)=x(mv,jdex(n))
	end if
	v(n)=xt(n)
	
	end do

	do k=1,ns
	ncase(k)=k
	end do
	
	npop=ndend-ndstart+1
	
	call quicksort(v,ncase,ndstart,ndend,ns)
	
	sts=0
	do n=1,npop
	sts=sts+w(n)
	end do

	do n=ndstart,ndend-1
	
	if(mv.gt.0) then
	ndl=ndl+ct(ncase(n))
	ndr=ndr-ct(ncase(n))
	
	tml=tml+tt(ncase(n))
	tmr=tmr-tt(ncase(n))
	
	end if
	
	       
	if(mv.eq.0) then
	if(v(n).gt.0) then
	ndl=ndl+ct(ncase(n))
	ndr=ndr-ct(ncase(n))
	end if
	z=.5*(v(n)+v(n+1))
	sl=sl+v(n)
	sr=sr-v(n)
	nst=n-ndstart
	tml=sl+(npop-nst-1)*z
	tmr=sr-(npop-nst-1)*z
	
	end if
	
	
	if(ndl*ndr.ge.1) then
	if (dmin1(tml,tmr).ge..1) then
	if (v(n).lt.v(n+1)) then
	s=ndl*dlog(1.0*ndl/tml)+ndr*dlog(1.0*ndr/tmr)
	if(s.gt.critvar) then
	critvar=s
	nbest=n
	ubest=.5*(v(n)+v(n+1))
	ndlt=ndl
	ndrt=ndr
	
	tmlt=tml
	tmrt=tmr
	mbest=mv
	
	
	do k=ndstart,ndend
	ut(k)=xt(k)
	end do
	end if
	end if
	end if
	end if
	end do

70 	continue
	
	if(mv.ge.1.and.cat(mv).ge.2) then
	
	call izerv(idiecat,32)
	call zervr(ttcat,32)
	do n=ndstart,ndend

	l=dnint(x(mv,jdex(n)))
	if(l.le.0.or.l.gt.32) write(*,*) "category error"
	idiecat(l)=idiecat(l)+ct(n)
	ttcat(l)=ttcat(l)+tt(n)
	end do
	lcat=cat(mv)
	do n=1, (2**(lcat-1))-1
	call unpack(lcat,n,icat)
	ndl=0
	tml=0
	do l=1,lcat
	if(icat(l).eq.1) then
	ndl=ndl+idiecat(l)
	tml=tml+ttcat(l)
	end if
	end do
	ndr=ndnd-ndl
	tmr=tmnd-tml
	if(ndl*ndr.ge.1) then
	if(dmin1(tml,tmr).ge.1.0e-6) then
	crit=ndl*dlog(1.0*ndl/tml)+ndr*dlog(1.0*ndr/tmr)
	if(crit.gt.critvar) then
	critvar=crit
	ubest=n
	mbest=mv
	ndlt=ndl
	ndrt=ndr
	tmlt=tml
	tmrt=tmr
	end if
	end if
	end if
	end do !n
	end if
90 	continue
	if(mv.eq.0) call exit
	end do !mdim
	
	
c	&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

c	CHANGE THE INDEXING

	if(critvar.le.-1.0e9) then
	jstat=1
	return
	end if
	
	critinc=critvar-crit0
	jcure=0
	if(mbest.eq.0) ubest=ubest+tmint
	if(jcure.eq.1) then
	nl=ndstart-1
	do nsp=ndstart,ndend
	if(td(nsp).le.ubest.and.td(n).gt.tmint.
     1	or.td(n).gt.tmaxt) then
     	nl=nl+1
	ncase(nl)=jdex(nsp)
	end if
	end do
	
	ndendl=max(nl,ndstart)
	nr=ndendl
	do nsp=ndstart,ndend
	if(ut(nsp).gt.ubest) then
	nr=nr+1
	if(nr.gt.ns) goto 375
	ncase(nr)=jdex(nsp)
	end if
	end do
	end if 
375	continue
	
     
	
	if(mbest.ge.1) then
	if (cat(mbest).eq.1) then
	nl=ndstart-1
	do nsp=ndstart,ndend
	if(ut(nsp).le.ubest) then
	nl=nl+1
	ncase(nl)=jdex(nsp)
	end if
	end do
	
	ndendl=max(nl,ndstart)
	nr=ndendl
	do nsp=ndstart,ndend
	if(ut(nsp).gt.ubest) then
	nr=nr+1
	if(nr.gt.ns) goto 765
	ncase(nr)=jdex(nsp)
	end if
	end do
        end if
        end if
765     continue

	if(mbest.ge.1) then
	if (cat(mbest).ge.2) then
	lcat=cat(mbest)
	n=nint(ubest)
	call unpack(lcat,n,icat)
	nl=ndstart-1
	do nsp=ndstart,ndend
	l=nint(x(mbest,jdex(nsp)))
	if(icat(l).eq.1) then
	nl=nl+1
	ncase(nl)=jdex(nsp)
	end if
	end do
	ndendl=max(nl,ndstart)
	nr=ndendl
	do nsp=ndstart,ndend
	l=nint(x(mbest,jdex(nsp)))
	if(icat(l).eq.0) then
	nr=nr+1
	if(nr.gt.ns) goto 775
	ncase(nr)=jdex(nsp)
	end if
	end do
775	continue
	end if !cat
	end if !mv>0
	
	if(mbest.ge.1) then
	if(ndendl.eq.ndend) ndendl=ndend-1
	do n=ndstart,ndend
	jdex(n)=ncase(n)
	end do
	end if
	
	
	end
	
	
	
	subroutine where(x,xsplit,nalive,msplit,treemap,
     1	nodestatus,ns,mdim,nrnodes,ndbigtree,cat,noc)
     	implicit double precision (a-h,o-z)
	
	double precision x(mdim,ns),xsplit(nrnodes)
     
     
	integer nalive(nrnodes),msplit(nrnodes),noc(nrnodes,ns),
     1	treemap(2,nrnodes),nodestatus(nrnodes),cat(0:mdim),icat(32)
     	
	call zerm(noc,nrnodes,ns)
	
	do n=1,ns
	
	call izerv(nalive,nrnodes)
	nalive(1)=1
	
	do kt=1,ndbigtree
	if(nalive(kt).eq.1.and.nodestatus(kt).eq.1) then
	nalive(kt)=0
	mv=msplit(kt)
	ubest=xsplit(kt)

	if(mv.ge.1) then
	if(cat(mv).eq.1) then
	if (x(mv,n).le.ubest) then
	kt1=treemap(1,kt)
	else
	kt1=treemap(2,kt)
	end if
	end if
	
	if(cat(mv).ge.2) then
	lcat=cat(mv)
	np=nint(xsplit(kt))
	call unpack(lcat,np,icat)
	l=nint(x(mv,n))
	if(icat(l).eq.1) then
	kt1=treemap(1,kt)
	else
	kt1=treemap(2,kt)
	end if
	end if
	
	if(nodestatus(kt1).eq.-1) then
	nalive(kt1)=0
	noc(kt1,n)=1
     	else
	nalive(kt1)=1
	end if
	end if  !mv
	
	if(mv.eq.0) then
	kt1=treemap(1,kt)
	kt2=treemap(2,kt)
	
	if(nodestatus(kt1).eq.-1) then
	nalive(kt1)=0
	noc(kt1,n)=1
	else
	nalive(kt1)=1
	end if
	
	if(nodestatus(kt2).eq.-1) then
	nalive(kt2)=0
	noc(kt2,n)=1
	else
	nalive(kt2)=1
	end if
	end if
	
	end if
	end do !kt
	end do !ns
	end
	
	subroutine surv(eml,nodestatus,jin,noc,tmax,tmin,
     1	nd,tm,td,ndbigtree,ns,nrnodes,tim)
     	implicit double precision (a-h,o-z)
     
     	double precision eml(ns),tmin(nrnodes),tmax(nrnodes),td(ns),
     1	tm(nrnodes)
     
     	integer noc(nrnodes,ns),jin(ns),nd(nrnodes),
     1	nodestatus(nrnodes)
     
     	call zervr(eml,ns)
	do k=1,ndbigtree
	if(nodestatus(k).eq.-1) then
	do n=1,ns
	if(jin(n).eq.0.and.noc(k,n).eq.1) eml(n)=eml(n)+
     1	xplus(dmin1(tim,tmax(k))-tmin(k))*nd(k)/tm(k)
     	end do
	end if
	end do
	end
	
	subroutine errorpred(agtime,vt,ns,ntsm,azip,gin,td,
     &	ce,xt,ut,std)
     	implicit double precision(a-h,o-z)
     
	double precision vt(ntsm),agtime(ntsm,ns),azip(ns,ntsm),
     &	xt(ns),ut(ns),td(ns),gin(ns)
     	integer ce(ns)
	

	do jj=1,ntsm-1
	do n=1,ns
	if(jj.gt.ntsm/2) then
	azip(n,jj)=vt(ntsm)
	else
	azip(n,jj)=vt(1)
	end if
	end do
	end do
	
	
	do jj=1,ntsm-1
	rg=dble(jj)/ntsm
	arg=-dlog(rg)
	do n=1,ns
	do k=1,ntsm-1
	if(agtime(k,n).le.arg.and.agtime(k+1,n).gt.arg) then
	azip(n,jj)=vt(k)
	goto 78
	end if
	end do !k
78	continue
	end do !n
	end do !jj
	
	
	do n=1,ns
	gin(n)=0
	do jj=1,ntsm-1
	gin(n)=gin(n)+azip(n,jj)
	end do
	gin(n)=gin(n)/(ntsm-1)
	end do
	
	
	std=0
	ntp=0
	tv=0
	do n=1,ns
	if(ce(n).eq.1) then
	ntp=ntp+1
	if(ntp.ge.1) then
	xt(ntp)=gin(n)
	ut(ntp)=td(n)
	std=std+2*dabs(xt(ntp)-ut(ntp))/(xt(ntp)+ut(ntp))
	end if
	end if
	end do
	std=std/ntp
	end
	
	
	
	subroutine proxim(prox,nvec,noc,ns,nrnodes)
	implicit double precision (a-h,o-z)
	double precision prox(ns,ns)
        integer nvec(nrnodes),noc(nrnodes,ns)
     	call izerv(nvec,nrnodes)
	do n1=1,ns
        nj=0
	do k=1,nrnodes
	if(noc(k,n1).eq.1) then
	nj=nj+1
	nvec(nj)=k
	end if
	end do
	do n2=1,ns
	isame=0
	do jt=1,nj
	if(noc(nvec(jt),n2).eq.1) then
	isame=isame+1
	call exit
	end if
	end do
	prox(n1,n2)=prox(n1,n2)+isame
	end do
	end do

	end 
	
c       SUNROUTINE  SCALE

	subroutine scale(s,xsc,y,u,dl,ns,nn,red)
	implicit double precision (a-h,o-z)
	double precision s(ns,ns),y(ns),u(ns),dl(ns),
     1	xsc(ns,nn),red(ns)	
	
	do j=1,ns
	red(j)=0
	do i=1,ns
	red(j)=red(j)+s(i,j)
	end do
	red(j)=red(j)/ns
	end do
	
	sred=0
	do j=1,ns
	sred=sred+red(j)
	end do
	sred=sred/ns
	
	do i=1,ns
	do j=1,ns
	s(i,j)=(s(i,j)-red(j)-red(i)+sred)/2
	end do
	end do
	
      
      	do it=1,nn
	
	do n=1,ns
	y(n)=2*(rand(1)-.5)
	end do
	
	
	
	do jit=1,1000
	
	y2=0
	do n=1,ns
	y2=y2+y(n)**2
	end do
	y2=dsqrt(y2)
	do n=1,ns
	u(n)=y(n)/y2
	end do
	
	
	do n=1,ns
	y(n)=0
	do k=1,ns
	y(n)=y(n)+s(n,k)*u(k)
	end do
	
	end do
	
	
	ra=0
	do n=1,ns
	ra=ra+y(n)*u(n)
	end do
	ynorm=0
	do n=1,ns
	ynorm=ynorm+(y(n)-ra*u(n))**2
	
	end do
	
	
	if(ynorm.lt.1.0e-7*ra)then
	do n=1,ns
	xsc(n,it)=dsqrt(ra)*u(n)
	
	end do
	dl(it)=ra
	call exit
	end if
	end do
	
	do n=1,ns
	do k=1,ns
	s(n,k)=s(n,k)-ra*u(n)*u(k)
	end do
	end do
	write(*,*) it,real(ra)
	end do !nn
	
	end
	
		


   
	
	
	subroutine regress(agtime,vt,ntsm,mdim,ns,cov,b,y,x,sgn,
     1	xy,xt,yt,xtt,wtt,ytt,ztt,ncx,ivncx,sot,nox,ivnox,cat,
     1	ut,ireg,avrsq,avb,atb,trsq)
     
   
     	implicit double precision(a-h,o-z)
	double precision b(mdim),cov(mdim,mdim),y(ns),sgn(mdim),
     1	x(mdim,ns),agtime(ntsm,ns),xt(ns),yt(ns),ztt(ns),ut(ns),
     1	xy(mdim),vt(ntsm),xtt(ns),ytt(ns),sot(mdim,ns),wtt(ns),
     1	avrsq(ntsm),avb(ntsm,mdim),zx(10,1000),b0(20),atb(mdim)
     
     	integer ncx(ns),ivncx(ns),nox(mdim,ns),ivnox(mdim,ns),
     1	cat(0:mdim)
	
	
	do m=1,mdim
	do n=1,ns
	sot(m,n)=x(m,n)
	end do
	end do
	
	
	
	call varnorm(x,mdim,ns)
	
	
	
	call zermr(cov,mdim,mdim)
	do i=1,mdim
	do j=1,mdim
	do n=1,ns
	cov(i,j)=cov(i,j)+x(i,n)*x(j,n)
	end do
	end do
	end do
	
	do j=1,mdim
	cov(j,j)=cov(j,j)*(1+.001)
	end do
	
	call rev(cov,mdim)

	do k=1,ntsm-1
	
     	do n=1,ns
	y(n)=(agtime(k,n))
	end do
	
	call stand(y,ns)
	
	do m=1,mdim
	xy(m)=0
	do n=1,ns
	xy(m)=xy(m)+y(n)*x(m,n)
	end do
	end do
	
	do m=1,mdim
	b(m)=0
	do j=1,mdim
	b(m)=b(m)+cov(m,j)*xy(j)
	end do
	end do

	
	c=0
	do m=1,mdim
	c=c+xy(m)*b(m)
	end do
	rsq=c/ns
	
	
	do n=1,ns
	xt(n)=0
	do m=1,mdim
	xt(n)=xt(n)+b(m)*x(m,n)
	zx(m,n)=x(m,n)
	end do
	yt(n)=y(n)-xt(n)
	end do
	
	do jit=1,3
	do m=1,mdim
	do n=1,ns
	yt(n)=yt(n)+b(m)*zx(m,n)
	end do
	if(cat(m).ge.2)  then
	do n=1,ns
	xtt(n)=sot(m,n)
	end do
	ndm=cat(m)
	call catress(xtt,yt,ut,ndm,ns,rss)
	call stand(ut,ns)
	b(m)=0
	do n=1,ns
	b(m)=b(m)+yt(n)*ut(n)
	zx(m,n)=ut(n)
	end do
	b(m)=b(m)/ns
	do n=1,ns
	yt(n)=yt(n)-b(m)*zx(m,n)
	end do
	end if
	
	if(cat(m).eq.1) then
	b(m)=0
	do n=1,ns
	b(m)=b(m)+yt(n)*zx(m,n)
	end do
	b(m)=b(m)/ns
	do n=1,ns
	yt(n)=yt(n)-b(m)*zx(m,n)
	end do
	end if
	
	end do !m
	end do !jit
	rss=0
	do n=1,ns
	rss=rss+yt(n)*yt(n)
	end do
	
	
	rsq=1-(rss/ns)
	if(ireg.eq.1) goto 223
	
	
	
c	begin ireg=2 &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
	do m=1,mdim
	
	if(cat(m).eq.1) then
	do n=1,ns
	xtt(n)=x(m,n)
	ncx(n)=n
	end do
	call quicksort(xtt,ncx,1,ns,ns)
	call revncase(ncx,ns,ivncx)
	do n=1,ns
	nox(m,n)=ncx(n)
	ivnox(m,n)=ivncx(n)
	end do
	end if
	
	if(cat(m).ge.2) then
	do n=1,ns
	nox(m,n)=n
	ivnox(m,n)=n
	end do
	end if
	
	end do !m
	
	do m=1,mdim
	do n=1,ns
	zx(m,n)=b(m)*zx(m,n)
	end do 
	end do
	

	do n=1,ns
	xt(n)=0
	do m=1,mdim
	xt(n)=xt(n)+zx(m,n)
	end do
	end do
	do n=1,ns
	yt(n)=y(n)-xt(n)
	end do
	
	rss=0
	do n=1,ns
	rss=rss+yt(n)*yt(n)
	end do
	
	do n=1,ns
	ytt(n)=yt(n)
	end do
	
	do jit=1,5
	
	do j=1,mdim
	rss1=0
	do n=1,ns
	ytt(n)=ytt(nox(j,n))+zx(j,nox(j,n))
	ztt(n)=ytt(n)
	wtt(n)=-ytt(n)
	end do 
	
	if(cat(j).ge.2)  then
	do n=1,ns
	xtt(n)=sot(j,n)
	end do
	ndm=cat(j)
	call catress(xtt,ytt,ut,ndm,ns,rss0)
	if(rss0.ge.rss) then
	do n=1,ns
	ytt(n)=ytt(n)-zx(j,n)
	end do
	goto 999
	else
	do n=1,ns
	zx(j,n)=ut(n)
	end do
	sgn(j)=0
	do n=1,ns
	sgn(j)=sgn(j)+ytt(n)*ut(n)
	ytt(n)=ytt(n)-ut(n)
	end do
	goto 999
	end if
	end if
	
	
	call cmon(ns,ztt)
	rssp=0
	do n=1,ns
	rssp=rssp+(ytt(n)-ztt(n))**2
	end do
	
	call cmon(ns,wtt)
	rssm=0
	do n=1,ns
	wtt(n)=-wtt(n)
	rssm=rssm+(ytt(n)-wtt(n))**2
	end do
	if(dmin1(rssm,rssp).ge.rss) then
     	do n=1,ns
	ytt(n)=ytt(n)-zx(j,nox(j,n))
	end do
	rss=rss
	goto 999
	end if
	if(rssm.le.rssp) then
	rss=rssm
	sgn(j)=0
	do n=1,ns
	sgn(j)=sgn(j)+ytt(n)*wtt(n)
	ytt(n)=ytt(n)-wtt(n)
	zx(j,n)=wtt(ivnox(j,n))
	end do
	else
	rss=rssp
	sgn(j)=0
	do n=1,ns
	sgn(j)=sgn(j)+ytt(n)*ztt(n)
	ytt(n)=ytt(n)-ztt(n)
	zx(j,n)=ztt(ivnox(j,n))
	end do
	end if
	
999	continue
	
	do n=1,ns
	ztt(n)=ytt(ivnox(j,n))
	end do
	do n=1,ns
	ytt(n)=ztt(n)
	end do
	end do !j
	end do !jit
	
	do m=1,mdim
	do n=1,ns
	xt(n)=zx(m,n)
	end do
	call standv(xt,ns,sd)
	if(cat(m).eq.1) then
	if (b(m).gt.0) then
	bt=sd 
	else
	bt=-sd
	end if
	b(m)=bt
	else
	if (sgn(m).gt.0) then
	bt=sd
	else
	bt=-sd
	end if
	b(m)=bt
	end if
	end do
	rsq=1-(rss/ns)
223	continue

	
	
	avrsq(k)=rsq
	trsq=trsq+rsq
	do m=1,mdim
	avb(k,m)=b(m)
	atb(m)=atb(m)+b(m)
	end do
	if(k.eq.ntsm-1) call exit
	end do !k
	trsq=trsq/(ntsm-1)
	do m=1,mdim
	atb(m)=atb(m)/(ntsm-1)
	end do
	
	end
	
	subroutine corr(agtime,x,vt,ntsm,mdim,ns,
     1  ncx,xt,yt,acor,ty,ytt,zt,ut,cat,icor,sot,avcor)
     
     	
	implicit double precision(a-h,o-z)
	double precision x(mdim,ns),agtime(ntsm,ns),
     1	xt(ns),yt(ns),vt(ntsm),acor(ntsm,mdim),ty(ns),
     1	ytt(ns),ut(ns),zz(100),zt(ns),sot(mdim,ns),
     1	avcor(mdim)
     
     
     	integer ncx(ns),cat(0:mdim)
	
	do m=1,mdim
	do n=1,ns
	sot(m,n)=x(m,n)
	end do
	end do
	
	call varnorm(x,mdim,ns)
	
	if(icor.eq.1) then
	do m=1,mdim
	do n=1,ns
	zt(n)=x(m,n)
	end do
	do k=1,ntsm
	do n=1,ns
	yt(n)=agtime(k,n)
	end do
	call stand(yt,ns)
	if(cat(m).ge.2) then
	do n=1,ns
	zt(n)=sot(m,n)
	end do
	ndm=cat(m)
	call acatc(zt,yt,ns,cc,ndm)
	else
	cc=0
	do n=1,ns
	cc=cc+zt(n)*yt(n)
	end do
	cc=cc/ns
	end if
	acor(k,m)=cc
	avcor(m)=avcor(m)+cc
	end do !k
	
	end do !m
	do m=1,mdim
	avcor(m)=avcor(m)/(ntsm-1)
	end do
	end if
	
	if(icor.eq.2) then
	
	do m=1,mdim
	if(cat(m).eq.1) then
	do n=1,ns
	xt(n)=x(m,n)
	zt(n)=xt(n)
	ncx(n)=n
	end do
	call quicksort(xt,ncx,1,ns,ns)
	end if
	
	if(cat(m).ge.2) then
	do n=1,ns
	zt(n)=x(m,n)
	end do
	end if
	
	
	
	do k=1,ntsm
	do n=1,ns
	yt(n)=agtime(k,n)
	end do
	call stand(yt,ns)
	
	ndm=cat(m)
	if(cat(m).gt.1) call acatc(zt,yt,ns,cc,ndm)
	if(cat(m).eq.1) then
	do n=1,ns
	ytt(n)=yt(ncx(n))
	ty(n)=ytt(n)
	end do
	call cmon(ns,ty)
	call stand(ty,ns)
	cc1=0
	do n=1,ns
	cc1=cc1+ytt(n)*ty(n)
	end do
	do n=1,ns
	ty(n)=-ytt(n)
	end do
	call cmon(ns,ty)
	call stand(ty,ns)
	cc2=0
	do n=1,ns
	cc2=cc2-ytt(n)*ty(n)
	end do
	cc=(cc1-cc2)/ns
	end if
	acor(k,m)=cc
	avcor(m)=avcor(m)+cc
	end do !k
	end do !m
	do m=1,mdim
	avcor(m)=avcor(m)/(ntsm-1)	
	end do
	end if
	end
	
	subroutine cmon(jr,rlamm)
	implicit double precision (a-h,o-z)
	dimension rlamm(jr)
	
	do j=2,jr
	if(rlamm(j-1).gt.rlamm(j)) then
	k1=0
	av=rlamm(j)
	do k=1,j-1
	k1=k1+1
	av=(k*av+rlamm(j-k))/(k+1)
	if (j.ge.3) then
	if(k.le.j-2) then
	
	if(av.ge.rlamm(j-k-1)) call exit
	end if
	end if
	end do
	do i=j-k1,j
	rlamm(i)=av
	end do
	end if
	end do
	
	end
	
	subroutine revncase(ncase,ns,ivncase)
	integer ncase(ns),ivncase(ns)
	do k=1,ns
	do n=1,ns
	if(ncase(n).eq.k) ivncase(k)=n
	end do
	end do
	end

	
	subroutine reducetime(ce,td,otd,vt,ntsm,ns,ncase)
	implicit double precision (a-h,o-z)
	double precision td(ns),otd(ns),vt(ntsm)
	integer ce(ns),ncase(ns)
	
	tmx=0
	ntt=0
	call zervr(otd,ns)
	do n=1,ns
	if(ce(n).eq.1) then
	ntt=ntt+1
	otd(ntt)=td(n)
	end if
	end do
	ndie=ntt
	call quicksort(otd,ncase,1,ndie,ns)
	ntt=max(ntt,ntsm)
	do k=1,ntsm
	vt(k)=otd(idnint(dble(k*ntt)/ntsm))
	end do
	end
	
	subroutine sortdie(x,td,jce,xt,ncase,mdim,ns)
     	implicit double precision(a-h,o-z)
	
	double precision x(mdim,ns),td(ns),xt(ns)
	integer jce(ns),ncase(ns)
	
	do n=1,ns
	ncase(n)=n
	end do
	call quicksort(td,ncase,1,ns,ns)
	
	do m=1,mdim
	do n=1,ns
	xt(n)=x(m,n)
	end do
	do n=1,ns
	x(m,n)=xt(ncase(n))
	end do
	end do
	
	do n=1,ns
	xt(n)=jce(n)
	end do
	do n=1,ns
	jce(n)=nint(xt(ncase(n)))
	end do
	end
	
	subroutine acatc(xt,yt,ns,cc,ndm)
	implicit double precision(a-h,o-z)
	parameter(mc=10)
	double precision xt(ns),yt(ns),s(mc),p(mc)
     
     	call stand(yt,ns)
	
	do j=1,ndm
	s(j)=0
	p(j)=0
	do n=1,ns
	if(idnint(xt(n)).eq.j)then
	s(j)=s(j)+yt(n)
	p(j)=p(j)+1
	end if
	end do
	end do
	
	cc=0
	do j=1,ndm
	if(p(j).gt.0) then
	s(j)=s(j)/p(j)
	else
	s(j)=0
	end if
	p(j)=p(j)/ns
	cc=cc+s(j)*s(j)*p(j)
	end do
	
	cc=dsqrt(cc)
	
	end
	
	subroutine catress(xt,yt,ut,ndm,ns,rss)
	implicit double precision(a-h,o-z)
	
	double precision xt(ns),yt(ns),ut(ns)
	parameter(nsm=32)
	double precision s(nsm),nc(nsm)
	
	
	do j=1,ndm
	s(j)=0
	nc(j)=0
	do n=1,ns
	if(idnint(xt(n)).eq.j)then
	s(j)=s(j)+yt(n)
	nc(j)=nc(j)+1
	end if
	end do
	if(nc(j).ge.1) then
	s(j)=s(j)/nc(j)
	else
	s(j)=0
	end if
	end do
	
	rss=0
	do n=1,ns
	ut(n)=s(idnint(xt(n)))
	rss=rss+(yt(n)-ut(n))**2
	end do
	end
	
	
	

	
	subroutine rev(sw,m2)
	implicit double precision(a-h,o-z)
	double precision sw(m2,m2)
	integer im(0:200)
	iflag=0
	do 5 k=1,m2
	im(k)=0
5	continue
	do 10 ii= 1,m2
	rmax=0
	do 20 k=1,m2
	if(im(k).eq.0) then
		if(abs(sw(k,k)).gt.rmax) then
			ll=k
			rmax=abs(sw(k,k))
		end if
	end if
20	continue
	im(ll)=1
	td=sw(ll,ll)
	if (td.eq.0.0) then 
		write(*,*) "zero-inverse",ll
		iflag=1
	return
	end if
	do 30 j= 1,m2
	sw(ll,j)=sw(ll,j)/td
30	continue 
	do 40 i= 1,m2
	if (i.eq.ll) goto 40
	ct=sw(i,ll)
	do 50 j= 1,m2
	sw(i,j)=sw(i,j)-ct*sw(ll,j)
50	continue 
	sw(i,ll)=-ct/td
40	continue 
	sw(ll,ll)=1/td
10	continue

	end
	
	

	
	 subroutine unpack(l,npack,icat)
	
c	npack is a long integer.  The sub. returns icat, an integer of zeroes and
c	ones corresponding to the coefficients in the binary expansion of npack.
	
	integer icat(32)
	call izerv(icat,32)
	n=npack
	icat(1)=mod(n,2)
	do 10 k=2,l
	n=(n-icat(k-1))/2
	icat(k)=mod(n,2)
10	continue
	end
	
	subroutine stand(z,ns)
	implicit double precision(a-h,o-z)
	double precision z(ns)
	
	av=0
	sq=0
	do n=1,ns
	av=av+z(n)
	sq=sq+z(n)*z(n)
	end do
	av=av/ns
	sq=sq/ns
	sd=dsqrt(sq-av*av)
	if(sd.le.1.0e-6) then
	do n=1,ns
	z(n)=0
	end do
	else
	do n=1,ns
	z(n)=(z(n)-av)/sd
	end do
	end if
	end
	
	subroutine standv(z,ns,sd)
	implicit double precision(a-h,o-z)
	double precision z(ns)
	
	av=0
	sq=0
	do n=1,ns
	
	av=av+z(n)
	sq=sq+z(n)*z(n)
	end do
	av=av/ns
	sq=sq/ns
	sd=dsqrt(sq-av*av)
	if(sd.le.1.0e-6) then
	do n=1,ns
	z(n)=0
	end do
	else
	do n=1,ns
	z(n)=(z(n)-av)/sd
	end do
	end if
	end
	
	
      
      	
	
	subroutine varnorm(x,mdim,ns)
	implicit double precision(a-h,o-z)
	double precision x(mdim,ns)
	
	
		do m=1,mdim
		avt=0
		sqt=0
		do n=1,ns
		sqt=sqt+((n-1)*(x(m,n)-avt)**2)/n
		avt=((n-1)*avt+x(m,n))/n
		end do
		sdt=dsqrt(sqt/ns)
		do n=1,ns
		x(m,n)=(x(m,n)-avt)/sdt
		end do
		end do
		end
		
	
	
c	MISCELLANOUS SMALL SUBROUTINES

	double precision function xplus(x)
	double precision x
	if(x.gt.0) then
	xplus=x
	else
	xplus=0
	end if
	end
	
	subroutine izerv(ix,m1)
	integer ix(m1)
	do 10 n=1,m1
	ix(n)=0
10	continue
	end
	
	subroutine zervr(rx,m1)
	double precision rx(m1)
	do 10 n=1,m1
	rx(n)=0
10	continue
	end
	
	subroutine zerm(mx,m1,m2)
	integer mx(m1,m2)
	do 10 i=1,m1
	do 20 j=1,m2
	mx(i,j)=0
20	continue
10	continue
	end
	
	subroutine zermr(rx,m1,m2)
	double precision  rx(m1,m2)
	do 10 i=1,m1
	do 20 j=1,m2
	rx(i,j)=0
20	continue
10	continue
	end
	
	
	
        double precision function rand(j)
	double precision dseed,u
	save dseed
	data dseed /17395/
	call lrnd(dseed,u)
	rand=u
	end

	subroutine lrnd(dseed,u)
	double precision dseed, d31m1,u
	data d31m1 /2147483647/
	dseed=dmod(16087*dseed,d31m1)
	u=dseed/d31m1
	return
	end
	
	double precision function rnorm(j)
        implicit double precision (a-h,o-z)
	double precision u,v
        u=rand(1)
	v=rand(1)
	rnorm=dsqrt(-2*dlog(u))*dcos(6.28318531*v)
	end 
	
		subroutine quicksort (v,iperm,ii,jj,kk)
c
c     puts into iperm the permutation vector which sorts v into
c     increasing order.  only elementest from ii to jj are considered.
c     array iu(k) and array il(k) permit sorting up to 2**(k+1)-1 elements
c
c     this is a modification of acm algorithm #347 by r. c. singleton,
c     which is a modified hoare quicksort.
      implicit double precision(a-h,o-z)
      integer iperm(kk),iu(64),il(64)
      double precision v(kk)
      integer t,tt
      
      m=1
      i=ii
      j=jj
 10   if (i.ge.j) go to 80
 20   k=i
      ij=(j+i)/2
      t=iperm(ij)
      vt=v(ij)
      if (v(i).le.vt) go to 30
      iperm(ij)=iperm(i)
      iperm(i)=t
      t=iperm(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
 30   l=j
      if (v(j).ge.vt) go to 50
      iperm(ij)=iperm(j)
      iperm(j)=t
      t=iperm(ij)
      v(ij)=v(j)
      v(j)=vt
      vt=v(ij)
      if (v(i).le.vt) go to 50
      iperm(ij)=iperm(i)
      iperm(i)=t
      t=iperm(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
      go to 50
 40   iperm(l)=iperm(k)
      iperm(k)=tt
      v(l)=v(k)
      v(k)=vtt
 50   l=l-1
      if (v(l).gt.vt) go to 50
      tt=iperm(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=iperm(i+1)
      vt=v(i+1)
      if (v(i).le.vt) go to 100
      k=i
 110  iperm(k+1)=iperm(k)
      v(k+1)=v(k)
      k=k-1
      if (vt.lt.v(k)) go to 110
      iperm(k+1)=t
      v(k+1)=vt
      go to 100
      end
      
	
	
	
	
