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(11, file='vet-lung',status='old')

	do n=1,ns

	read(11,*) 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


	

	

	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) goto 977

	end do !mdim
	
977	continue
	

	

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

	goto 977

	end if

	end do
	
977	continue

	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

	goto 977

	end if

	end do

977	continue
	

	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) goto 977

	end do !k
	
977	continue

	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)) goto 977

	end if

	end if

	end do
	
977	continue

	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