| C |
| C Modified from the SMART package by J.H. Friedman, 10/10/84 |
| C Main change is to add spline smoothing modified from BRUTO, |
| C calling code written for smooth.spline in S. |
| C |
| C B.D. Ripley (ripley@stats.ox.ac.uk) 1994-7. |
| C |
| C |
| subroutine smart(m,mu,p,q,n, w,x,y,ww,smod,nsmod, |
| & sp,nsp,dp,ndp,edf) |
| |
| integer m,mu,p,q,n, nsmod, nsp,ndp |
| double precision x(p,n),y(q,n),w(n),ww(q),smod(nsmod), |
| & sp(nsp),edf(m),dp(ndp) |
| smod(1)=m |
| smod(2)=p |
| smod(3)=q |
| smod(4)=n |
| call smart1(m,mu,p,q,n, w,x,y,ww, smod(6),smod(q+6), |
| & smod(q+7),smod(q+7+p*m),smod(q+7+m*(p+q)), |
| & smod(q+7+m*(p+q+n)),smod(q+7+m*(p+q+2*n)), |
| & sp,sp(q*n+1),sp(n*(q+15)+1),sp(n*(q+15)+q+1), |
| & dp,smod(5),edf) |
| return |
| end |
| |
| subroutine smart1(m,mu,p,q,n, w,x,y,ww, yb,ys, |
| & a,b,f, |
| & t,asr, |
| & r,sc,bt,g, |
| & dp,flm,edf) |
| |
| integer m,mu,p,q,n |
| double precision w(n),x(p,n),y(*),ww(q), yb(q), ys |
| double precision a(p,m),b(q,m),f(n,m),t(n,m), asr(15),asr1 |
| double precision r(q,n),sc(n,15),bt(q),g(p,3) |
| double precision dp(*), flm,edf(m) |
| C ^^^ really (ndb) of smart(.) |
| integer i,j,l, lm |
| double precision sw,s |
| c Common Vars |
| double precision span,alpha,big |
| integer ifl,lf |
| common /pprpar/ ifl,lf,span,alpha,big |
| |
| double precision conv, cutmin,fdel,cjeps |
| integer maxit,mitone, mitcj |
| common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj |
| |
| sw=0d0 |
| do j=1,n |
| sw=sw+w(j) |
| end do |
| do j=1,n |
| do i=1,q |
| r(i,j)=y(q*(j-1)+i) |
| end do |
| end do |
| do i=1,q |
| s=0d0 |
| do j=1,n |
| s=s+w(j)*r(i,j) |
| end do |
| yb(i)=s/sw |
| end do |
| c yb is vector of means |
| do j=1,n |
| do i=1,q |
| r(i,j)=r(i,j)-yb(i) |
| end do |
| end do |
| ys=0.d0 |
| do i=1,q |
| s=0.d0 |
| do j=1,n |
| s=s+w(j)*r(i,j)**2 |
| end do |
| ys=ys+ww(i)*s/sw |
| end do |
| if(ys .gt. 0d0) goto 311 |
| c ys is the overall standard deviation -- quit if zero |
| return |
| |
| 311 continue |
| ys=sqrt(ys) |
| s=1.d0/ys |
| do j=1,n |
| do i=1,q |
| r(i,j)=r(i,j)*s |
| end do |
| end do |
| |
| c r is now standardized residuals |
| c subfit adds up to m terms one at time; lm is the number fitted. |
| call subfit(m,p,q,n,w,sw,x,r,ww,lm,a,b,f,t,asr(1),sc,bt,g,dp,edf) |
| if(lf.le.0) go to 9999 |
| call fulfit(lm,lf,p,q,n,w,sw,x,r,ww,a,b,f,t,asr,sc,bt,g,dp,edf) |
| C REPEAT |
| 371 continue |
| do l=1,lm |
| sc(l,1)=l+0.1d0 |
| s=0d0 |
| do i=1,q |
| s=s+ww(i)*abs(b(i,l)) |
| end do |
| sc(l,2)=-s |
| end do |
| call sort(sc(1,2),sc,1,lm) |
| do j=1,n |
| do i=1,q |
| r(i,j)=y(q*(j-1)+i) |
| end do |
| end do |
| |
| do i=1,q |
| do j=1,n |
| r(i,j)=r(i,j)-yb(i) |
| s=0.d0 |
| do l=1,lm |
| s=s+b(i,l)*f(j,l) |
| end do |
| r(i,j)=r(i,j)/ys-s |
| end do |
| end do |
| |
| if(lm.le.mu) goto 9999 |
| c back to integer: |
| l=int(sc(lm,1)) |
| asr1=0d0 |
| do j=1,n |
| do i=1,q |
| r(i,j)=r(i,j)+b(i,l)*f(j,l) |
| asr1=asr1+w(j)*ww(i)*r(i,j)**2 |
| end do |
| end do |
| |
| asr1=asr1/sw |
| asr(1)=asr1 |
| if(l .ge. lm) goto 591 |
| do i=1,p |
| a(i,l)=a(i,lm) |
| end do |
| do i=1,q |
| b(i,l)=b(i,lm) |
| end do |
| do j=1,n |
| f(j,l)=f(j,lm) |
| t(j,l)=t(j,lm) |
| end do |
| |
| 591 continue |
| lm=lm-1 |
| call fulfit(lm,lf,p,q,n,w,sw,x,r,ww,a,b,f,t,asr,sc,bt,g,dp,edf) |
| goto 371 |
| C END REPEAT |
| 9999 continue |
| flm=lm |
| return |
| end |
| |
| subroutine subfit(m,p,q,n,w,sw,x,r,ww,lm,a,b,f,t,asr,sc, |
| & bt,g,dp,edf) |
| c Args |
| integer m,p,q,n, lm |
| double precision w(n),sw, x(p,n),r(q,n),ww(q),a(p,m),b(q,m), |
| & f(n,m), t(n,m), asr(15), sc(n,15), bt(q), g(p,3), edf(m) |
| double precision dp(*) |
| c Var |
| integer i,j,l, iflsv |
| double precision asrold |
| c Common Vars |
| double precision span,alpha,big |
| integer ifl,lf |
| common /pprpar/ ifl,lf,span,alpha,big |
| |
| double precision conv, cutmin,fdel,cjeps |
| integer maxit,mitone, mitcj |
| common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj |
| |
| asr(1)=big |
| lm=0 |
| do 100 l=1,m |
| call rchkusr() |
| lm=lm+1 |
| asrold=asr(1) |
| call newb(lm,q,ww,b) |
| c does 'edf' mean 'edf(1)' or 'edf(l)'? |
| call onetrm(0,p,q,n,w,sw,x,r,ww,a(1,lm),b(1,lm), |
| & f(1,lm),t(1,lm),asr(1),sc,g,dp,edf(1)) |
| do 20 j=1,n |
| do 10 i=1,q |
| r(i,j)=r(i,j)-b(i,lm)*f(j,lm) |
| 10 continue |
| 20 continue |
| if(lm.eq.1) goto 100 |
| if(lf.gt.0) then |
| if(lm.eq.m) return |
| iflsv=ifl |
| ifl=0 |
| call fulfit(lm,1,p,q,n,w,sw,x,r,ww,a,b,f,t,asr,sc,bt, |
| & g,dp, edf) |
| ifl=iflsv |
| endif |
| if(asr(1).le.0d0.or.(asrold-asr(1))/asrold.lt.conv) return |
| 100 continue |
| return |
| end |
| |
| subroutine fulfit(lm,lbf,p,q,n,w,sw,x,r,ww,a,b,f,t, |
| & asr,sc,bt,g,dp,edf) |
| c Args |
| integer lm,lbf,p,q,n |
| double precision w(n),sw,x(p,n),r(q,n),ww(q),a(p,lm),b(q,lm), |
| & f(n,lm),t(n,lm),asr(1+lm), sc(n,15),bt(q),g(p,3), edf(lm) |
| double precision dp(*) |
| c Var |
| double precision asri, fsv, asrold |
| integer i,j,iter,lp,isv |
| c Common Vars |
| double precision span,alpha,big |
| integer ifl,lf |
| common /pprpar/ ifl,lf,span,alpha,big |
| |
| double precision conv, cutmin,fdel,cjeps |
| integer maxit,mitone, mitcj |
| common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj |
| |
| if(lbf.le.0) return |
| asri=asr(1) |
| fsv=cutmin |
| isv=mitone |
| if(lbf .lt. 3) then |
| cutmin=1d0 |
| mitone=lbf-1 |
| endif |
| iter=0 |
| C Outer loop: |
| 1000 continue |
| asrold=asri |
| iter=iter+1 |
| do 100 lp=1,lm |
| do 10 i=1,q |
| bt(i)=b(i,lp) |
| 10 continue |
| do 20 i=1,p |
| g(i,3)=a(i,lp) |
| 20 continue |
| do 35 j=1,n |
| do 30 i=1,q |
| r(i,j)=r(i,j)+bt(i)*f(j,lp) |
| 30 continue |
| 35 continue |
| |
| call onetrm(1,p,q,n,w,sw,x,r,ww,g(1,3),bt,sc(1,14),sc(1,15), |
| & asri,sc,g,dp,edf(lp)) |
| if(asri .lt. asrold) then |
| do 40 i=1,q |
| b(i,lp)=bt(i) |
| 40 continue |
| do 50 i=1,p |
| a(i,lp)=g(i,3) |
| 50 continue |
| do 60 j=1,n |
| f(j,lp)=sc(j,14) |
| t(j,lp)=sc(j,15) |
| 60 continue |
| else |
| asri=asrold |
| endif |
| do 85 j=1,n |
| do 80 i=1,q |
| r(i,j)=r(i,j)-b(i,lp)*f(j,lp) |
| 80 continue |
| 85 continue |
| 100 continue |
| if((iter .le. maxit) .and. ((asri .gt. 0d0) |
| & .and. ((asrold-asri)/asrold .ge. conv))) goto 1000 |
| cutmin=fsv |
| mitone=isv |
| if(ifl .gt. 0) then |
| asr(1+lm) = asri |
| asr(1) = asri |
| endif |
| return |
| end |
| |
| subroutine onetrm(jfl,p,q,n,w,sw,x,y,ww,a,b,f,t,asr, |
| & sc,g,dp,edf) |
| c Args |
| integer jfl,p,q,n |
| double precision w(n),sw, x(p,n),y(q,n),ww(q),a(p),b(q),f(n),t(n), |
| & asr, sc(n,13),g(p,2), edf |
| double precision dp(*) |
| c Var |
| double precision asrold,s |
| integer i,j,iter |
| c Common Vars |
| double precision span,alpha,big |
| integer ifl,lf |
| common /pprpar/ ifl,lf,span,alpha,big |
| |
| double precision conv, cutmin,fdel,cjeps |
| integer maxit,mitone, mitcj |
| common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj |
| |
| iter=0 |
| asr=big |
| C REPEAT |
| 1000 continue |
| iter=iter+1 |
| asrold=asr |
| do 11 j=1,n |
| s=0d0 |
| do 21 i=1,q |
| s=s+ww(i)*b(i)*y(i,j) |
| 21 continue |
| sc(j,13)=s |
| 11 continue |
| call oneone(max0(jfl,iter-1),p,n,w,sw,sc(1,13),x,a,f,t, |
| & asr,sc,g,dp,edf) |
| do 31 i=1,q |
| s=0d0 |
| do 41 j=1,n |
| s=s+w(j)*y(i,j)*f(j) |
| 41 continue |
| b(i)=s/sw |
| 31 continue |
| asr=0d0 |
| do 51 i=1,q |
| s=0d0 |
| do 61 j=1,n |
| s=s+w(j)*(y(i,j)-b(i)*f(j))**2 |
| 61 continue |
| asr=asr+ww(i)*s/sw |
| 51 continue |
| if((q .ne. 1) .and. (iter .le. maxit) .and. (asr .gt. 0d0) |
| & .and. (asrold-asr)/asrold .ge. conv) goto 1000 |
| return |
| end |
| |
| subroutine oneone(ist,p,n, w,sw,y,x,a,f,t,asr,sc,g,dp,edf) |
| c Args |
| integer ist,p,n |
| double precision w(n),sw,y(n),x(p,n),a(p),f(n),t(n),asr, |
| & sc(n,12), g(p,2), edf, dp(*) |
| c Var |
| integer i,j,k,iter |
| double precision sml, s,v,cut,asrold |
| c Common Vars |
| double precision span,alpha,big |
| integer ifl,lf |
| common /pprpar/ ifl,lf,span,alpha,big |
| |
| double precision conv, cutmin,fdel,cjeps |
| integer maxit,mitone, mitcj |
| common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj |
| |
| sml=1d0/big |
| if(ist .le. 0) then |
| if(p .le. 1) a(1)=1d0 |
| do 10 j=1,n |
| sc(j,2)=1d0 |
| 10 continue |
| call pprdir(p,n,w,sw,y,x,sc(1,2),a,dp) |
| endif |
| s=0d0 |
| do 20 i=1,p |
| g(i,1)=0d0 |
| s=s+a(i)**2 |
| 20 continue |
| s=1d0/sqrt(s) |
| do 30 i=1,p |
| a(i)=a(i)*s |
| 30 continue |
| iter=0 |
| asr=big |
| cut=1d0 |
| C REPEAT ----------------------------- |
| 100 continue |
| iter=iter+1 |
| asrold=asr |
| C REPEAT [inner loop] ----- |
| 60 continue |
| s=0d0 |
| do 70 i=1,p |
| g(i,2)=a(i)+g(i,1) |
| s=s+g(i,2)**2 |
| 70 continue |
| s=1.d0/sqrt(s) |
| do 80 i=1,p |
| g(i,2)=g(i,2)*s |
| 80 continue |
| do 90 j=1,n |
| sc(j,1)=j+0.1d0 |
| s=0.d0 |
| do 91 i=1,p |
| s=s+g(i,2)*x(i,j) |
| 91 continue |
| sc(j,11)=s |
| 90 continue |
| call sort(sc(1,11),sc,1,n) |
| do 110 j=1,n |
| k=int(sc(j,1)) |
| sc(j,2)=y(k) |
| sc(j,3)=max(w(k),sml) |
| 110 continue |
| call supsmu(n,sc(1,11),sc(1,2),sc(1,3),1,span,alpha, |
| & sc(1,12),sc(1,4), edf) |
| s=0d0 |
| do 120 j=1,n |
| s=s+sc(j,3)*(sc(j,2)-sc(j,12))**2 |
| 120 continue |
| s=s/sw |
| if(s .lt. asr) goto 140 |
| cut=cut*0.5d0 |
| if(cut.lt.cutmin) goto 199 |
| do 150 i=1,p |
| g(i,1)=g(i,1)*cut |
| 150 continue |
| go to 60 |
| C -------- |
| 140 continue |
| asr=s |
| cut=1d0 |
| do 160 i=1,p |
| a(i)=g(i,2) |
| 160 continue |
| do 170 j=1,n |
| k=int(sc(j,1)) |
| t(k)=sc(j,11) |
| f(k)=sc(j,12) |
| 170 continue |
| if(asr.le.0d0.or.(asrold-asr)/asrold.lt.conv) goto 199 |
| if(iter.gt.mitone.or.p.le.1) goto 199 |
| call pprder(n,sc(1,11),sc(1,12),sc(1,3),fdel,sc(1,4),sc(1,5)) |
| do 180 j=1,n |
| k=int(sc(j,1)) |
| sc(j,5)=y(j)-f(j) |
| sc(k,6)=sc(j,4) |
| 180 continue |
| call pprdir(p,n,w,sw,sc(1,5),x,sc(1,6),g,dp) |
| |
| goto 100 |
| c-------------- |
| 199 continue |
| c-------------- |
| s=0d0 |
| v=s |
| do 210 j=1,n |
| s=s+w(j)*f(j) |
| 210 continue |
| s=s/sw |
| do 220 j=1,n |
| f(j)=f(j)-s |
| v=v+w(j)*f(j)**2 |
| 220 continue |
| if(v .gt. 0d0) then |
| v=1d0/sqrt(v/sw) |
| do 230 j=1,n |
| f(j)=f(j)*v |
| 230 continue |
| endif |
| return |
| end |
| |
| |
| subroutine pprdir(p,n,w,sw,r,x,d,e,g) |
| |
| integer p,n |
| double precision w(n),sw,r(n),x(p,n),d(n),e(p), g(*) |
| |
| double precision s |
| integer i,j,k,l,m1,m2 |
| |
| double precision conv, cutmin,fdel,cjeps |
| integer maxit,mitone, mitcj |
| common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj |
| |
| do 10 i=1,p |
| s=0d0 |
| do 15 j=1,n |
| s=s+w(j)*d(j)*x(i,j) |
| 15 continue |
| e(i)=s/sw |
| 10 continue |
| k=0 |
| m1=p*(p+1)/2 |
| m2=m1+p |
| do 20 j=1,p |
| s=0d0 |
| do 22 l=1,n |
| s=s+w(l)*r(l)*(d(l)*x(j,l)-e(j)) |
| 22 continue |
| g(m1+j)=s/sw |
| do 25 i=1,j |
| s=0d0 |
| do 27 l=1,n |
| s=s+w(l)*(d(l)*x(i,l)-e(i))*(d(l)*x(j,l)-e(j)) |
| 27 continue |
| k=k+1 |
| g(k)=s/sw |
| 25 continue |
| 20 continue |
| call ppconj(p,g,g(m1+1),g(m2+1),cjeps,mitcj,g(m2+p+1)) |
| do 30 i=1,p |
| e(i)=g(m2+i) |
| 30 continue |
| return |
| end |
| |
| subroutine ppconj(p,g,c,x,eps,maxit,sc) |
| integer p,maxit |
| double precision g(*),c(p),x(p),eps,sc(p,4) |
| |
| integer i,j,im1,iter,nit |
| double precision beta,h,s,alpha,t |
| |
| do 1 i=1,p |
| x(i)=0d0 |
| sc(i,2)=0d0 |
| 1 continue |
| nit=0 |
| C REPEAT |
| 11321 continue |
| nit=nit+1 |
| h=0d0 |
| beta=0d0 |
| do 11331 i=1,p |
| sc(i,4)=x(i) |
| s=g(i*(i-1)/2+i)*x(i) |
| im1=i-1 |
| j=1 |
| goto 11343 |
| 11341 j=j+1 |
| 11343 if(j.gt.im1) goto 11342 |
| s=s+g(i*(i-1)/2+j)*x(j) |
| goto 11341 |
| 11342 continue |
| j=i+1 |
| goto 11353 |
| 11351 j=j+1 |
| 11353 if(j.gt.p) goto 11352 |
| s=s+g(j*(j-1)/2+i)*x(j) |
| goto 11351 |
| 11352 continue |
| sc(i,1)=s-c(i) |
| h=h+sc(i,1)**2 |
| 11331 continue |
| if(h.le.0d0) goto 11322 |
| do 11361 iter=1,p |
| do 11371 i=1,p |
| sc(i,2)=beta*sc(i,2)-sc(i,1) |
| 11371 continue |
| t=0d0 |
| do 11381 i=1,p |
| s=g(i*(i-1)/2+i)*sc(i,2) |
| im1=i-1 |
| j=1 |
| goto 11393 |
| 11391 j=j+1 |
| 11393 if(j.gt.im1) goto 11392 |
| s=s+g(i*(i-1)/2+j)*sc(j,2) |
| goto 11391 |
| 11392 continue |
| j=i+1 |
| goto 11403 |
| 11401 j=j+1 |
| 11403 if(j.gt.p) goto 11402 |
| s=s+g(j*(j-1)/2+i)*sc(j,2) |
| goto 11401 |
| 11402 continue |
| sc(i,3)=s |
| t=t+s*sc(i,2) |
| 11381 continue |
| alpha=h/t |
| s=0d0 |
| do 11411 i=1,p |
| x(i)=x(i)+alpha*sc(i,2) |
| sc(i,1)=sc(i,1)+alpha*sc(i,3) |
| s=s+sc(i,1)**2 |
| 11411 continue |
| if(s.le.0d0) goto 11362 |
| beta=s/h |
| h=s |
| 11361 continue |
| 11362 continue |
| s=0d0 |
| do 11421 i=1,p |
| s=dmax1(s,dabs(x(i)-sc(i,4))) |
| 11421 continue |
| if((s .ge. eps) .and. (nit .lt. maxit)) goto 11321 |
| 11322 continue |
| return |
| end |
| |
| subroutine pprder (n,x,s,w,fdel,d,sc) |
| integer n |
| double precision x(n),s(n),w(n), fdel, d(n),sc(n,3) |
| |
| integer i,j,bl,el,bc,ec,br,er |
| double precision scale, del |
| c |
| c unnecessary initialization of bl el ec to keep g77 -Wall happy |
| c |
| bl = 0 |
| el = 0 |
| ec = 0 |
| c |
| if(x(n) .gt. x(1)) goto 11441 |
| do 11451 j=1,n |
| d(j)=0d0 |
| 11451 continue |
| return |
| 11441 continue |
| i=n/4 |
| j=3*i |
| scale=x(j)-x(i) |
| 11461 if(scale.gt.0d0) goto 11462 |
| if(j.lt.n) j=j+1 |
| if(i.gt.1) i=i-1 |
| scale=x(j)-x(i) |
| goto 11461 |
| 11462 continue |
| del=fdel*scale*2d0 |
| do 11471 j=1,n |
| sc(j,1)=x(j) |
| sc(j,2)=s(j) |
| sc(j,3)=w(j) |
| 11471 continue |
| call pool (n,sc,sc(1,2),sc(1,3),del) |
| bc=0 |
| br=bc |
| er=br |
| 11481 continue |
| br=er+1 |
| er=br |
| 11491 if(er .ge. n) goto 11492 |
| if(sc(br,1) .ne. sc(er+1,1)) goto 11511 |
| er=er+1 |
| goto 11521 |
| 11511 continue |
| goto 11492 |
| 11521 continue |
| goto 11491 |
| 11492 continue |
| if(br .ne. 1) goto 11541 |
| bl=br |
| el=er |
| goto 11481 |
| 11541 continue |
| if(bc .ne. 0) goto 11561 |
| bc=br |
| ec=er |
| do 11571 j=bl,el |
| d(j)=(sc(bc,2)-sc(bl,2))/(sc(bc,1)-sc(bl,1)) |
| 11571 continue |
| goto 11481 |
| 11561 continue |
| c sanity check needed for PR#13517 |
| if(br.gt.n) call rexit('br is too large') |
| do 11581 j=bc,ec |
| d(j)=(sc(br,2)-sc(bl,2))/(sc(br,1)-sc(bl,1)) |
| 11581 continue |
| if(er .ne. n) goto 11601 |
| do 11611 j=br,er |
| d(j)=(sc(br,2)-sc(bc,2))/(sc(br,1)-sc(bc,1)) |
| 11611 continue |
| goto 11482 |
| 11601 continue |
| bl=bc |
| el=ec |
| bc=br |
| ec=er |
| goto 11481 |
| 11482 continue |
| return |
| end |
| |
| subroutine pool (n,x,y,w,del) |
| integer n |
| double precision x(n),y(n),w(n),del |
| |
| integer i,bb,eb,br,er,bl,el |
| double precision px, py, pw |
| |
| bb=0 |
| eb=bb |
| 11621 if(eb.ge.n) goto 11622 |
| bb=eb+1 |
| eb=bb |
| 11631 if(eb .ge. n) goto 11632 |
| if(x(bb) .ne. x(eb+1)) goto 11651 |
| eb=eb+1 |
| goto 11661 |
| 11651 continue |
| goto 11632 |
| 11661 continue |
| goto 11631 |
| 11632 continue |
| if(eb .ge. n) goto 11681 |
| if(x(eb+1)-x(eb) .ge. del) goto 11701 |
| br=eb+1 |
| er=br |
| 11711 if(er .ge. n) goto 11712 |
| if(x(er+1) .ne. x(br)) goto 11731 |
| er=er+1 |
| goto 11741 |
| 11731 continue |
| goto 11712 |
| 11741 continue |
| goto 11711 |
| 11712 continue |
| C avoid bounds error: this was .and. but order is not guaranteed |
| if(er.lt.n) then |
| if(x(er+1)-x(er).lt.x(eb+1)-x(eb)) goto 11621 |
| endif |
| eb=er |
| pw=w(bb)+w(eb) |
| px=(x(bb)*w(bb)+x(eb)*w(eb))/pw |
| py=(y(bb)*w(bb)+y(eb)*w(eb))/pw |
| do 11751 i=bb,eb |
| x(i)=px |
| y(i)=py |
| w(i)=pw |
| 11751 continue |
| 11701 continue |
| 11681 continue |
| 11761 continue |
| if(bb.le.1) goto 11762 |
| if(x(bb)-x(bb-1).ge.del) goto 11762 |
| bl=bb-1 |
| el=bl |
| 11771 if(bl .le. 1) goto 11772 |
| if(x(bl-1) .ne. x(el)) goto 11791 |
| bl=bl-1 |
| goto 11801 |
| 11791 continue |
| goto 11772 |
| 11801 continue |
| goto 11771 |
| 11772 continue |
| bb=bl |
| pw=w(bb)+w(eb) |
| px=(x(bb)*w(bb)+x(eb)*w(eb))/pw |
| py=(y(bb)*w(bb)+y(eb)*w(eb))/pw |
| do 11811 i=bb,eb |
| x(i)=px |
| y(i)=py |
| w(i)=pw |
| 11811 continue |
| goto 11761 |
| 11762 continue |
| goto 11621 |
| 11622 continue |
| return |
| end |
| |
| subroutine newb(lm,q,ww,b) |
| integer lm, q |
| double precision ww(q),b(q,lm) |
| |
| integer i,lm1,l,l1 |
| double precision s,t,sml |
| c Common |
| double precision span,alpha,big |
| integer ifl,lf |
| common /pprpar/ ifl,lf,span,alpha,big |
| |
| |
| sml=1d0/big |
| if(q .ne. 1) goto 11831 |
| b(1,lm)=1d0 |
| return |
| 11831 continue |
| if(lm .ne. 1) goto 11851 |
| do 11861 i=1,q |
| b(i,lm)=i |
| 11861 continue |
| return |
| 11851 continue |
| lm1=lm-1 |
| do 11871 i=1,q |
| b(i,lm)=0d0 |
| 11871 continue |
| t=0d0 |
| do 11881 i=1,q |
| s=0d0 |
| do 11891 l=1,lm1 |
| s=s+abs(b(i,l)) |
| 11891 continue |
| b(i,lm)=s |
| t=t+s |
| 11881 continue |
| do 11901 i=1,q |
| b(i,lm)=ww(i)*(t-b(i,lm)) |
| 11901 continue |
| l1=1 |
| if(lm.gt.q) l1=lm-q+1 |
| do 11911 l=l1,lm1 |
| s=0d0 |
| t=s |
| do 11921 i=1,q |
| s=s+ww(i)*b(i,lm)*b(i,l) |
| t=t+ww(i)*b(i,l)**2 |
| 11921 continue |
| s=s/sqrt(t) |
| do 11931 i=1,q |
| b(i,lm)=b(i,lm)-s*b(i,l) |
| 11931 continue |
| 11911 continue |
| do 11941 i=2,q |
| if(abs(b(i-1,lm)-b(i,lm)).gt.sml) return |
| 11941 continue |
| do 11951 i=1,q |
| b(i,lm)=i |
| 11951 continue |
| return |
| end |
| |
| block data bkppr |
| |
| c Common Vars |
| double precision span,alpha,big |
| integer ifl,lf |
| common /pprpar/ ifl,lf,span,alpha,big |
| |
| double precision conv, cutmin,fdel,cjeps |
| integer maxit,mitone, mitcj |
| common /pprz01/ conv,maxit,mitone,cutmin,fdel,cjeps,mitcj |
| |
| double precision df, gcvpen |
| integer ismethod |
| logical trace |
| common /spsmooth/ df, gcvpen, ismethod, trace |
| |
| data df, gcvpen, ismethod, trace /4d0, 1d0, 0, .false./ |
| |
| data ifl,maxit, conv, mitone, cutmin, fdel, |
| & span,alpha, big, cjeps, mitcj, lf |
| & /6, 20, .005, 20, 0.1, 0.02, |
| & 0.0, 0.0,1.0e20,0.001, 1, 2/ |
| end |
| |
| subroutine setppr(span1, alpha1, optlevel, ism, df1, gcvpen1) |
| c Put 'parameters' into Common blocks |
| integer optlevel,ism |
| double precision span1,alpha1, df1, gcvpen1 |
| |
| double precision span,alpha,big |
| integer ifl,lf |
| common /pprpar/ ifl,lf,span,alpha,big |
| |
| double precision df, gcvpen |
| integer ismethod |
| logical trace |
| common /spsmooth/ df, gcvpen, ismethod, trace |
| |
| span = span1 |
| lf = optlevel |
| alpha = alpha1 |
| if(ism .ge. 0) then |
| ismethod = ism |
| trace = .false. |
| else |
| ismethod = -(ism+1) |
| trace = .true. |
| end if |
| df = df1 |
| gcvpen = gcvpen1 |
| return |
| end |
| |
| subroutine fsort(mu,n,f,t,sp) |
| c |
| integer mu, n |
| double precision f(n,mu),t(n,mu),sp(n,2) |
| c |
| integer l,j,k |
| |
| do 100 l=1,mu |
| do 10 j=1,n |
| sp(j,1)=j+0.1d0 |
| sp(j,2)=f(j,l) |
| 10 continue |
| call sort(t(1,l),sp,1,n) |
| do 20 j=1,n |
| k=int(sp(j,1)) |
| f(j,l)=sp(k,2) |
| 20 continue |
| 100 continue |
| return |
| end |
| |
| subroutine pppred(np,x,smod,y,sc) |
| |
| integer np |
| double precision x(np,*),y(np,*),smod(*), sc(*) |
| |
| integer p,q, place,low,high, i,j,l,m,n, |
| + inp,ja,jb,jf,jt,jfl,jfh,jtl,jth, mu |
| double precision ys, s, t |
| |
| m= int(smod(1)+0.1d0) |
| p= int(smod(2)+0.1d0) |
| q= int(smod(3)+0.1d0) |
| n= int(smod(4)+0.1d0) |
| mu=int(smod(5)+0.1d0) |
| ys=smod(q+6) |
| ja=q+6 |
| jb=ja+p*m |
| jf=jb+m*q |
| jt=jf+n*m |
| call fsort(mu,n,smod(jf+1),smod(jt+1),sc) |
| do 100 inp = 1, np |
| ja=q+6 |
| jb=ja+p*m |
| jf=jb+m*q |
| jt=jf+n*m |
| do 81 i=1,q |
| y(inp,i)=0d0 |
| 81 continue |
| do 91 l=1,mu |
| s=0d0 |
| do 12201 j=1,p |
| s=s+smod(ja+j)*x(inp,j) |
| 12201 continue |
| if(s .gt. smod(jt+1)) goto 12221 |
| place=1 |
| go to 12230 |
| 12221 continue |
| if(s .lt. smod(jt+n)) goto 12251 |
| place=n |
| go to 12230 |
| |
| 12251 continue |
| low=0 |
| high=n+1 |
| C WHILE |
| 12261 if(low+1.ge.high) goto 12262 |
| place=(low+high)/2 |
| t=smod(jt+place) |
| if(s.eq.t) goto 12230 |
| if(s .lt. t) then |
| high=place |
| else |
| low=place |
| endif |
| goto 12261 |
| C END |
| 12262 continue |
| jfl=jf+low |
| jfh=jf+high |
| jtl=jt+low |
| jth=jt+high |
| t=smod(jfl)+(smod(jfh)-smod(jfl))*(s-smod(jtl)) / |
| & (smod(jth)-smod(jtl)) |
| go to 12300 |
| 12230 continue |
| t=smod(jf+place) |
| 12300 continue |
| do 12311 i=1,q |
| y(inp,i)=y(inp,i)+smod(jb+i)*t |
| 12311 continue |
| ja=ja+p |
| jb=jb+q |
| jf=jf+n |
| jt=jt+n |
| 91 continue |
| do 12321 i=1,q |
| y(inp,i)=ys*y(inp,i)+smod(i+5) |
| 12321 continue |
| 100 continue |
| return |
| end |
| |
| c Called from R's supsmu() |
| subroutine setsmu (tr) |
| integer tr |
| |
| double precision df, gcvpen |
| integer ismethod |
| logical trace |
| common /spsmooth/ df, gcvpen, ismethod, trace |
| |
| ismethod = 0 |
| trace = tr .ne. 0 |
| return |
| end |
| |
| subroutine supsmu (n,x,y,w,iper,span,alpha,smo,sc,edf) |
| c |
| c------------------------------------------------------------------ |
| c |
| c super smoother (Friedman, 1984). |
| c |
| c version 10/10/84 |
| c |
| c coded and copywrite (c) 1984 by: |
| c |
| c Jerome H. Friedman |
| c department of statistics |
| c and |
| c stanford linear accelerator center |
| c stanford university |
| c |
| c all rights reserved. |
| c |
| c |
| c input: |
| c n : number of observations (x,y - pairs). |
| c x(n) : ordered abscissa values. |
| c y(n) : corresponding ordinate (response) values. |
| c w(n) : weight for each (x,y) observation. |
| c iper : periodic variable flag. |
| c iper=1 => x is ordered interval variable. |
| c iper=2 => x is a periodic variable with values |
| c in the range (0.0,1.0) and period 1.0. |
| c span : smoother span (fraction of observations in window). |
| c span=0.0 <=> "cv" : automatic (variable) span selection. |
| c alpha : controls high frequency (small span) penality |
| c used with automatic span selection (bass tone control). |
| c (alpha.le.0.0 or alpha.gt.10.0 => no effect.) |
| c output: |
| c smo(n) : smoothed ordinate (response) values. |
| c scratch: |
| c sc(n,7) : internal working storage. |
| c |
| c note: |
| c for small samples (n < 40) or if there are substantial serial |
| c correlations between observations close in x - value, then |
| c a prespecified fixed span smoother (span > 0) should be |
| c used. reasonable span values are 0.2 to 0.4. |
| c |
| c------------------------------------------------------------------ |
| |
| c Args |
| integer n, iper |
| double precision x(n),y(n),w(n), smo(n),sc(n,7) |
| double precision span, alpha, edf |
| c Var |
| double precision sy,sw, a,h(n),f, scale,vsmlsq,resmin |
| integer i,j, jper |
| |
| double precision spans(3), big,sml,eps |
| common /spans/ spans /consts/ big,sml,eps |
| |
| double precision df, gcvpen |
| integer ismethod |
| logical trace |
| common /spsmooth/ df, gcvpen, ismethod, trace |
| c Called from R's supsmu(), ismethod = 0, always (but not when called from ppr) |
| |
| if (x(n).gt.x(1)) go to 30 |
| c x(n) <= x(1) : boundary case: smo[.] := weighted mean( y ) |
| sy=0d0 |
| sw=sy |
| do 10 j=1,n |
| sy=sy+w(j)*y(j) |
| sw=sw+w(j) |
| 10 continue |
| a=0d0 |
| if (sw.gt.0d0) a=sy/sw |
| do 20 j=1,n |
| smo(j)=a |
| 20 continue |
| return |
| |
| C Normal Case |
| 30 continue |
| if (ismethod .ne. 0) then ! possible only when called from ppr() |
| call spline(n, x, y, w, smo, edf, sc) |
| else |
| i=n/4 |
| j=3*i |
| scale=x(j)-x(i) ! = IQR(x) |
| 40 if (scale.gt.0d0) go to 50 |
| if (j.lt.n) j=j+1 |
| if (i.gt.1) i=i-1 |
| scale=x(j)-x(i) |
| go to 40 |
| 50 vsmlsq=(eps*scale)**2 |
| jper=iper |
| if (iper.eq.2.and.(x(1).lt.0d0.or.x(n).gt.1d0)) jper=1 |
| if (jper.lt.1.or.jper.gt.2) jper=1 |
| if (span .gt. 0d0) then |
| call smooth (n,x,y,w,span,jper,vsmlsq,smo,sc) |
| return |
| end if |
| C else "cv" (crossvalidation) from three spans[] |
| do 70 i=1,3 |
| call smooth (n,x,y,w,spans(i),jper,vsmlsq, |
| & sc(1,2*i-1),sc(1,7)) |
| call smooth (n,x,sc(1,7),w,spans(2),-jper,vsmlsq, |
| & sc(1,2*i),h) |
| 70 continue |
| do 90 j=1,n |
| resmin=big |
| do 80 i=1,3 |
| if (sc(j,2*i).ge.resmin) go to 80 |
| resmin=sc(j,2*i) |
| sc(j,7)=spans(i) |
| 80 continue |
| if (alpha.gt.0d0 .and. alpha.le.10d0 .and. |
| & resmin.lt.sc(j,6) .and. resmin.gt.0d0) |
| & sc(j,7)= sc(j,7)+(spans(3)-sc(j,7)) * |
| & max(sml,resmin/sc(j,6))**(10d0-alpha) |
| 90 continue |
| |
| call smooth (n,x,sc(1,7),w,spans(2),-jper,vsmlsq,sc(1,2),h) |
| do 110 j=1,n |
| if (sc(j,2).le.spans(1)) sc(j,2)=spans(1) |
| if (sc(j,2).ge.spans(3)) sc(j,2)=spans(3) |
| f=sc(j,2)-spans(2) |
| if (f.ge.0d0) go to 100 |
| f=-f/(spans(2)-spans(1)) |
| sc(j,4)=(1d0-f)*sc(j,3)+f*sc(j,1) |
| go to 110 |
| 100 f=f/(spans(3)-spans(2)) |
| sc(j,4)=(1d0-f)*sc(j,3)+f*sc(j,5) |
| 110 continue |
| call smooth (n,x,sc(1,4),w,spans(1),-jper,vsmlsq,smo,h) |
| edf = 0 |
| endif |
| return |
| end |
| |
| subroutine smooth (n,x,y,w,span,iper,vsmlsq,smo,acvr) |
| c Args |
| integer n, iper |
| double precision x(n),y(n),w(n), span,vsmlsq, smo(n),acvr(n) |
| c Var |
| integer i,j, in,out, jper,ibw,it, j0 |
| double precision xm,ym,var,cvar, fbw,fbo,xti,xto,tmp, a,h,sy,wt |
| |
| c will use 'trace': |
| double precision df, gcvpen |
| integer ismethod |
| logical trace |
| common /spsmooth/ df, gcvpen, ismethod, trace |
| |
| xm=0d0 |
| ym=xm |
| var=ym |
| cvar=var |
| fbw=cvar |
| jper=iabs(iper) |
| ibw=int(0.5d0*span*n+0.5d0) |
| if (ibw.lt.2) ibw=2 |
| it=2*ibw+1 |
| if (it .gt. n) it = n |
| do i=1,it |
| j=i |
| if (jper.eq.2) j=i-ibw-1 |
| if (j.ge.1) then |
| xti=x(j) |
| else ! if (j.lt.1) then |
| j=n+j |
| xti=x(j)-1d0 |
| end if |
| wt=w(j) |
| fbo=fbw |
| fbw=fbw+wt |
| if (fbw.gt.0d0) xm=(fbo*xm+wt*xti)/fbw |
| if (fbw.gt.0d0) ym=(fbo*ym+wt*y(j))/fbw |
| tmp=0d0 |
| if (fbo.gt.0d0) tmp=fbw*wt*(xti-xm)/fbo |
| var =var +tmp*(xti-xm) |
| cvar=cvar+tmp*(y(j)-ym) |
| end do |
| |
| do 80 j=1,n |
| out=j-ibw-1 |
| in=j+ibw |
| if ((jper.ne.2) .and. (out.lt.1.or.in.gt.n)) go to 60 |
| if (out .lt. 1) then |
| out=n+out |
| xto=x(out)-1d0 |
| xti=x(in) |
| else if (in .gt. n) then |
| in=in-n |
| xti=x(in)+1d0 |
| xto=x(out) |
| else |
| xto=x(out) |
| xti=x(in) |
| end if |
| wt=w(out) |
| fbo=fbw |
| fbw=fbw-wt |
| tmp=0d0 |
| if (fbw.gt.0d0) tmp=fbo*wt*(xto-xm)/fbw |
| var = var-tmp*(xto-xm) |
| cvar=cvar-tmp*(y(out)-ym) |
| if (fbw.gt.0d0) xm=(fbo*xm-wt*xto)/fbw |
| if (fbw.gt.0d0) ym=(fbo*ym-wt*y(out))/fbw |
| wt=w(in) |
| fbo=fbw |
| fbw=fbw+wt |
| if (fbw.gt.0d0) xm=(fbo*xm+wt*xti)/fbw |
| if (fbw.gt.0d0) ym=(fbo*ym+wt*y(in))/fbw |
| tmp=0d0 |
| if (fbo.gt.0d0) tmp=fbw*wt*(xti-xm)/fbo |
| var = var+tmp*(xti-xm) |
| cvar=cvar+tmp*(y(in)-ym) |
| 60 a=0d0 |
| if (var.gt.vsmlsq) a=cvar/var |
| smo(j)=a*(x(j)-xm)+ym |
| if (iper.gt.0) then |
| h=0d0 |
| if (fbw.gt.0d0) h=1d0/fbw |
| if (var.gt.vsmlsq) h=h+(x(j)-xm)**2/var |
| acvr(j)=0d0 |
| a=1d0-w(j)*h |
| if (a.gt.0d0) then |
| acvr(j)=abs(y(j)-smo(j))/a |
| else if (j.gt.1) then |
| acvr(j)=acvr(j-1) |
| end if |
| end if |
| 80 continue |
| |
| if(trace) call smoothprt(span, iper, var, cvar) ! -> ./ksmooth.c |
| |
| c-- Recompute fitted values smo(j) as weighted mean for non-unique x(.) values: |
| j=1 |
| 90 j0=j |
| sy=smo(j)*w(j) |
| fbw=w(j) |
| if (j.ge.n) go to 110 |
| 100 if (x(j+1).le.x(j)) then |
| j=j+1 |
| sy=sy+w(j)*smo(j) |
| fbw=fbw+w(j) |
| if (j.lt.n) go to 100 |
| end if |
| 110 if (j.gt.j0) then |
| a=0d0 |
| if (fbw.gt.0d0) a=sy/fbw |
| do i=j0,j |
| smo(i)=a |
| end do |
| end if |
| j=j+1 |
| if (j.le.n) go to 90 |
| return |
| end |
| |
| |
| block data bksupsmu |
| double precision spans(3), big,sml,eps |
| common /spans/ spans /consts/ big,sml,eps |
| |
| data spans, big,sml,eps /0.05,0.2,0.5, 1.0e20,1.0e-7,1.0e-3/ |
| end |
| c--------------------------------------------------------------- |
| c |
| c this sets the compile time (default) values for various |
| c internal parameters : |
| c |
| c spans : span values for the three running linear smoothers. |
| c spans(1) : tweeter span. |
| c spans(2) : midrange span. |
| c spans(3) : woofer span. |
| c (these span values should be changed only with care.) |
| c big : a large representable floating point number. |
| c sml : a small number. should be set so that (sml)**(10.0) does |
| c not cause floating point underflow. |
| c eps : used to numerically stabilize slope calculations for |
| c running linear fits. |
| c |
| c these parameter values can be changed by declaring the |
| c relevant labeled common in the main program and resetting |
| c them with executable statements. |
| c |
| |
| c Only for ppr(*, ismethod != 0): Compute "smoothing" spline |
| c (rather, a penalized regression spline with at most 15 (inner) knots): |
| c----------------------------------------------------------------- |
| c |
| subroutine spline (n, x, y, w, smo, edf, sc) |
| c |
| c------------------------------------------------------------------ |
| c |
| c input: |
| c n : number of observations. |
| c x(n) : ordered abscissa values. |
| c y(n) : corresponding ordinate (response) values. |
| c w(n) : weight for each (x,y) observation. |
| c work space: |
| c sc(n,7) : used for dx(n), dy(n), dw(n), dsmo(n), lev(n) |
| c output: |
| c smo(n) : smoothed ordinate (response) values. |
| c edf : equivalent degrees of freedom |
| c |
| c------------------------------------------------------------------ |
| c Args |
| integer n |
| double precision x(n), y(n), w(n), smo(n), edf, sc(n,7) |
| |
| call splineAA(n, x, y, w, smo, edf, |
| + sc(n,1), ! dx |
| + sc(n,2), ! dy |
| + sc(n,3), ! dw |
| + sc(n,4), ! dsmo |
| + sc(n,5)) ! lev |
| |
| return |
| end |
| |
| |
| subroutine splineAA(n, x, y, w, smo, edf, dx, dy, dw, dsmo, lev) |
| c |
| c Workhorse of spline() above |
| c------------------------------------------------------------------ |
| c |
| c Additional input variables (no extra output, work): |
| c dx : |
| c dy : |
| c dw : |
| c dsmo: |
| c lev : "leverages", i.e., diagonal entries S_{i,i} of the smoother matrix |
| |
| c |
| c------------------------------------------------------------------ |
| c Args |
| integer n |
| double precision x(n), y(n), w(n), smo(n), edf, |
| + dx(n), dy(n), dw(n), dsmo(n), lev(n) |
| c Var |
| double precision knot(29), coef(25), work((17+25)*25) |
| double precision param(5), df1, lambda, crit, p, s |
| integer iparms(4), i, nk, ip, ier |
| |
| double precision df, gcvpen |
| integer ismethod |
| logical trace |
| common /spsmooth/ df, gcvpen, ismethod, trace |
| |
| c__no-more__ if (n .gt. 2500) call bdrsplerr() |
| do i = 1,n |
| dx(i) = (x(i)-x(1))/(x(n)-x(1)) |
| dy(i) = y(i) |
| dw(i) = w(i) |
| end do |
| nk = min(n,15) |
| knot(1) = dx(1) |
| knot(2) = dx(1) |
| knot(3) = dx(1) |
| knot(4) = dx(1) |
| knot(nk+1) = dx(n) |
| knot(nk+2) = dx(n) |
| knot(nk+3) = dx(n) |
| knot(nk+4) = dx(n) |
| do i = 5, nk |
| p = (n-1)*real(i-4)/real(nk-3) |
| ip = int(p) |
| p = p-ip |
| knot(i) = (1-p)*dx(ip+1) + p*dx(ip+2) |
| end do |
| c call dblepr('knots', 5, knot, nk+4) |
| C iparms(1:2) := (icrit, ispar) for ./sbart.c |
| if (ismethod .eq. 1) then |
| iparms(1) = 3 |
| df1 = df |
| else |
| iparms(1) = 1 |
| df1 = 0d0 |
| endif |
| c |
| iparms(2) = 0 ! ispar := 0 <==> estimate `spar' |
| iparms(3) = 500 ! maxit = 500 |
| iparms(4) = 0 ! spar (!= lambda) |
| c |
| param(1) = 0d0 ! = lspar : min{spar} |
| param(2) = 1.5d0 ! = uspar : max{spar} |
| c tol for 'spar' estimation: |
| param(3) = 1d-2 |
| c 'eps' (~= 2^-12 = sqrt(2^-24) ?= sqrt(machine eps)) in ./sbart.c : |
| param(4) = .000244 |
| |
| ier = 1 |
| call rbart(gcvpen,df1,dx,dy,dw,0.0d0,n,knot,nk,coef,dsmo,lev, |
| & crit,iparms,lambda,param, work,4,1,ier) |
| if(ier .gt. 0) call intpr('spline(.) TROUBLE:', 18, ier, 1) |
| do i = 1,n |
| smo(i) = dsmo(i) |
| end do |
| c call dblepr('smoothed',8, dsmo, n) |
| s = 0 |
| do i = 1, n |
| s = s + lev(i) |
| end do |
| edf = s |
| if(trace) call splineprt(df,gcvpen,ismethod, lambda, edf) |
| return |
| end |
| |
| |
| *********************************************************************** |
| |
| C=== This was 'sort()' in gamfit's mysort.f [or sortdi() in sortdi.f ] : |
| C |
| C=== FIXME: Translate to C and add to ../../../main/sort.c <<<<< |
| C |
| C a[] is double precision because the caller reuses a double (sometimes v[] itself!) |
| 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 |
| integer ii,jj |
| double precision v(*), a(jj) |
| c |
| integer iu(20),il(20) |
| integer t,tt, m,i,j,ij,k,l |
| double precision vt, vtt |
| |
| m=1 |
| i=ii |
| j=jj |
| 10 if (i.ge.j) go to 80 |
| 20 k=i |
| ij=(j+i)/2 |
| t=int(a(ij)) |
| vt=v(ij) |
| if (v(i).le.vt) go to 30 |
| a(ij)=a(i) |
| a(i)=t |
| t=int(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=int(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=int(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=int(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=int(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 |
