block data
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/terprm/maxtry,maxpro,ppconv,r2max
common/pltprm/nxplt,nyplt,plotrm
common/pursut/mode
common/inmode/inmode
common/irest/irest
common/consts/rbig,reps,ratio
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
data iprint/2/
data nprint,maxit/1,2/
data itap1,itap2,itap3,itap4/3*6,0/
data maxtry,maxpro,ppconv/5,10,0.15/
data nxplt,nyplt,plotrm/60,40,0./
data mode,inmode/2,2/
data irest/0/
data rbig,reps,ratio/1.e30,1.e-7,1.e-5/
data korder,maxkno,banfac/3,8,1.5/
data sambw,varbw,smofac(1),smofac(2),smofac(3)/0.5,1.0, 0.5,1.0,2.
*0/
end
subroutine plosur(direct)
dimension direct(1)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
common/iplopr/iplopr
iplopr=0
s=0.
s=s+direct(j)*predic((i-1)*npred+j)
stor4(i)=s
stor4(2*nobs+i)=smodel(predic((i-1)*npred+1))
call smooth(stor4(1),stor4(2*nobs+1),nobs,bandw,stor4(1),stor4(nob
*s+1), stor4(4*nobs+1),itag)
ntrim=plotrm*nobs
imin=ntrim+1
imax=nobs-ntrim
xmin=stor4(imin)
xmax=stor4(imax)
xmax=xmax+1.e-4*(xmax-xmin)
rmin=rbig
rmax=-rmin
s=stor4(2*nobs+itag(i))
if(s .ge. rmin)goto 10051
rmin=s
if(s .le. rmax)goto 10071
rmax=s
rmax=rmax+1.e-4*(rmax-rmin)
call psetup(min0(nxplt,100),xmin,xmax, min0(nyplt,50),rmin,rmax)
call plot(stor4(i),stor4(2*nobs+itag(i)))
call pprint
return
end
subroutine plospl(ipro)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
common/iplopr/iplopr
iplopr=ipro
write(itap3,10090)ipro,(axis((ipro-1)*npred+j),j=1,npred)
proj=0.
proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j)
stor4(i)=proj
itag(i)=i
call sort(stor4(1),itag,1,nobs)
ntrim=plotrm*nobs
imin=ntrim+1
imax=nobs-ntrim
xmin=stor4(imin)
xmax=stor4(imax)
xmax=xmax+1.e-4*(xmax-xmin)
r=bvalue(rkno((ipro-1)*maxkno+1),coeff((ipro-1)*maxkno+1),nkno(ipr
*o)-korder, korder,stor4(i),0)
stor4(nobs+i)=r
rmin=rbig
rmax=-rmin
if(stor4(nobs+i) .ge. rmin)goto 10151
rmin=stor4(nobs+i)
if(stor4(nobs+i) .le. rmax)goto 10171
rmax=stor4(nobs+i)
rmin=rmin-0.1*(rmax-rmin)
rmax=rmax+1.e-4*(rmax-rmin)
call psetup(min0(nxplt,100),xmin,xmax,min0(nyplt,50),rmin,rmax)
goto 10183
call plot(rkno((ipro-1)*maxkno+i),rmin)
goto 10181
call pprint
iplopr=0
write(itap1,10190)ipro,(axis((ipro-1)*npred+j),j=1,npred)
rmin=rbig
rmax=-rmin
r=stor4(3*nobs+itag(i))
if(r .ge. rmin)goto 10221
rmin=r
if(r .le. rmax)goto 10241
rmax=r
stor4(2*nobs+i)=r
rmin=rmin-0.1*(rmax-rmin)
rmax=rmax+1.e-4*(rmax-rmin)
call smooth(stor4(1),stor4(2*nobs+1),nobs,bandw,stor4(1),stor4(nob
*s+1),stor4(4*nobs+1),itag)
call psetup(min0(nxplt,100),xmin,xmax,min0(nyplt,50),rmin,rmax)
call plot(stor4(i),stor4(2*nobs+i))
goto 10263
call plot(rkno((ipro-1)*maxkno+i),rmin)
goto 10261
call pprint
return
*0(7(2x,f7.4)/))
*,f7.4)/))
end
subroutine save
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
if(itap4 .ne. 0)goto 10281
write(itap1,10290)itap4
return
no=nobs
np=npred
mk=maxkno
mp=maxpro
rewind itap4
write(itap4) npred,nobs,iprint,nprint,maxit,itap1,itap2,itap3,itap
*4, ntry,npro,maxtry,maxpro,ppconv,korder,maxkno,banfac,bandw, samb
*w,varbw,(smofac(i),i=1,3),nxplt,nyplt,plotrm
write(itap4) ((predic((j-1)*npred+i),i=1,np),j=1,no),(respon(i),i=
*1,no),((xlim(i,j), i=1,2),j=1,mp),((axis((j-1)*npred+i),i=1,np),j=
*1,mp),((rkno((j-1)*maxkno+i),i=1,mk), j=1,mp),((coeff((j-1)*maxkno
*+i),i=1,mk),j=1,mp),(stor4(3*nobs+i),i=1,no), (nkno(i),i=1,mp),asr
*old,mode,inmode,irest,rbig,reps,ratio
write(itap1,10290)itap4
return
end
subroutine restor
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
if(itap4 .eq. 0)goto 10311
rewind itap4
read(itap4) npred,nobs,iprint,nprint,maxit,itap1,itap2,itap3,itap4
*, ntry,npro,maxtry,maxpro,ppconv,korder,maxkno,banfac,bandw, sambw
*,varbw,(smofac(i),i=1,3),nxplt,nyplt,plotrm
no=nobs
np=npred
mk=maxkno
mp=maxpro
read(itap4) ((predic((j-1)*npred+i),i=1,np),j=1,no),(respon(i),i=1
*,no),((xlim(i,j), i=1,2),j=1,mp),((axis((j-1)*npred+i),i=1,np),j=1
*,mp),((rkno((j-1)*maxkno+i),i=1,mk), j=1,mp),((coeff((j-1)*maxkno+
*i),i=1,mk),j=1,mp),(stor4(3*nobs+i),i=1,no), (nkno(i),i=1,mp),asro
*ld,mode,inmode,irest,rbig,reps,ratio
irest=1
write(itap1,10290)itap4
return
end
subroutine masa
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
external smuuth
write(itap1,10090)
write(itap1,10190)nobs,npred
write(itap1,10320)mode,maxtry,maxpro,ppconv,maxit
write(itap1,10330)korder,maxkno,banfac
write(itap1,10340)iprint,nprint,plotrm
if(irest .ne. 0)goto 10361
ntry=1
npro=1
stor4(3*nobs+i)=respon(i)
inmode=mode
if(mode .ne. 3)goto 10391
inmode=1
sum=0.
sum=sum+respon(i)
sum=sum/float(nobs)
sstot=0.
sstot=sstot+(respon(i)-sum)**2
asrold=sstot/float(nobs)
write(itap1,10420)asrold
goto 10431
npro=npro+1
ntry=1
irest=0
bandw=banfac/(maxkno-korder)
write(itap1,10450)npro,ntry
if(inmode .ne. 1)goto 10471
call select(asr)
goto 10481
call gensta
call trimin(axis((npro-1)*npred+1),npred,asr,smuuth,maxit,nprint,
* stor4(1))
write(itap1,10490)asr,(axis((npro-1)*npred+i),i=1,npred)
call range(npro)
if(iprint .lt. 2)goto 10511
call plores
call knopla(i)
if(iprint .ne. 3)goto 10541
nki=itnkno(i)
write(itap1,10560)i,nki,(stor2(maxkno*maxpro+(i-1)*maxkno+j),j=1,n
*ki)
call fit(asr)
write(itap1,10570)asr
call accept(iacc,asr)
if(iacc .ne. 0)goto 10591
write(itap1,10600)
ntry=ntry+1
goto 10611
write(itap1,10620)
stor4(3*nobs+i)=stor4(4*nobs+i)
coeff((j-1)*maxkno+i)=stor2(2*maxpro*maxkno+(j-1)*maxkno+i)
rkno((j-1)*maxkno+i)=stor2(maxkno*maxpro+(j-1)*maxkno+i)
nkno(i)=itnkno(i)
npro=npro+1
ntry=1
asrold=asr
call decide(iterm)
if(iterm .ne. 1)goto 10681
goto 10442
goto 10441
npro=npro-1
if(iprint .lt. 2)goto 10701
call plospl(i)
return
*h (4/19/80)/ 24h0parameters for this run)
*ro ,i5/ 13h ppconv ,g12.6/ 13h maxit ,i5)
*fac ,g12.6)
*rm ,g12.6)
*6/ 20h solution direction:/20(7(2x,f7.4)/))
*sitions:/20(7(2x,g12.6)/))
end
function smodel(x)
dimension x(1)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
smod=0.
proj=0.
proj=proj+x(j)*axis((i-1)*npred+j)
if(proj .ge. xlim(1,i))goto 10751
proj=xlim(1,i)
if(proj .le. xlim(2,i))goto 10771
proj=xlim(2,i)
smod=smod+bvalue(rkno((i-1)*maxkno+1),coeff((i-1)*maxkno+1), nkno
*(i)-korder,korder,proj,0)
smodel=smod
return
end
subroutine gensta
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
if(ntry .ne. 1)goto 10791
fmin=rbig
axis((npro-1)*npred+i)=0.0
axis((npro-1)*npred+ipred)=1.0
asr=smuuth(axis((npro-1)*npred+1),npred)
if(asr .ge. fmin)goto 10831
ix=ipred
fmin=asr
if(npro .ne. 1)goto 10851
write(itap1,10090)ipred,asr
if(iprint .lt. 2)goto 10871
call range(npro)
call plores
axis((npro-1)*npred+i)=0.0
axis((npro-1)*npred+ix)=1.0
write(itap1,10190)ix,fmin
goto 10891
call genpts(1,npred,axis((npro-1)*npred+1))
sum=0.
sum=sum+axis((npro-1)*npred+i)**2
sum=sqrt(sum)
axis((npro-1)*npred+i)=axis((npro-1)*npred+i)/sum
write(itap1,10330)
return
*g12.6)
*ion;/ 6h asr= ,g12.6)
end
subroutine accept(iacc,asr)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
iacc=1
if((asrold-asr) .ge. (ppconv*asrold))goto 10931
iacc=0
return
end
subroutine decide(iterm)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
iterm=0
if(mode .ne. 1)goto 10951
if((npro .le. maxpro) .and. ((npro .le. npred) .and. (ntry .le. 1)
*))goto 10971
iterm=1
return
if(mode .ne. 2)goto 10991
if((npro .le. maxpro) .and. (ntry .le. maxtry))goto 11011
iterm=1
return
if(mode .ne. 3)goto 11031
if(inmode .ne. 1 .or. (npro .le. npred) .and. (ntry .le. 1))goto 1
*1051
inmode=2
ntry=1
if((npro .le. maxpro) .and. (ntry .le. maxtry))goto 11071
iterm=1
return
end
subroutine range(ipro)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
xmin=rbig
xmax=-xmin
proj=0.
proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j)
if(proj .ge. xmin)goto 11111
xmin=proj
if(proj .le. xmax)goto 11131
xmax=proj
xlim(1,ipro)=xmin
xlim(2,ipro)=xmax
return
end
subroutine select(asr)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
asrmin=rbig
imin=0
if(npro .le. 1)goto 11161
iused=0
goto 11173
if(axis((j-1)*npred+i) .le. 0.1)goto 11191
iused=1
goto 11172
goto 11171
if(iused .ne. 1)goto 11211
goto 11141
axis((npro-1)*npred+j)=0.
axis((npro-1)*npred+i)=1.
asr=smuuth(axis((npro-1)*npred+1),npred)
if(asr .ge. asrmin)goto 11241
asrmin=asr
imin=i
if(npro .ne. 1)goto 11261
write(itap1,10090)i,asr
if(iprint .lt. 2)goto 11281
call range(npro)
call plores
axis((npro-1)*npred+i)=0.
axis((npro-1)*npred+imin)=1.
asr=smuuth(axis((npro-1)*npred+1),npred)
*g12.6)
return
end
subroutine fit(asr)
dimension biatx(20)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
nbspli=0
nbspli=nbspli+itnkno(i)
nbspli=nbspli-(korder+1)*npro+1
xtx((j-1)*maxpro*maxkno+i)=0.
stor4(4*nobs+i)=0.
goto 11333
istart=0
stor2(j)=0.
goto 11353
proj=0.
proj=proj+axis((j-1)*npred+k)*predic((i-1)*npred+k)
call interv(stor2(maxkno*maxpro+(j-1)*maxkno+1),itnkno(j),proj,lef
*t,mflag)
call bsplvb(stor2(maxkno*maxpro+(j-1)*maxkno+1),korder,1,proj,left
*,biatx)
if(j .ne. 1)goto 11381
is=istart+left-korder
stor2(is+k)=biatx(k)
istart=istart+itnkno(1)-korder
goto 11401
is=istart+left-(korder+1)
if(left .ne. korder)goto 11421
stor2(is+k)=biatx(k)
goto 11441
stor2(is+k)=biatx(k)
istart=istart+itnkno(j)-(korder+1)
goto 11351
stor4(4*nobs+j)=stor4(4*nobs+j)+respon(i)*stor2(j)
xtx((k-1)*maxpro*maxkno+j)=xtx((k-1)*maxpro*maxkno+j)+stor2(j)*sto
*r2(k)
goto 11331
call tred2(maxkno*maxpro,nbspli,xtx,stor4(nobs+1),stor4(1),xtx)
call tql2(maxkno*maxpro,nbspli,stor4(nobs+1),stor4(1),xtx,ierr)
if(ierr .eq. 0)goto 11491
write(itap1,10090)npro,ntry,ierr
stop
if(iprint .ne. 3)goto 11511
write(itap1,10190)(stor4(nobs+i),i=1,nbspli)
small=ratio*stor4(nobs+nbspli)
ibig=1
ibig=ibig+1
goto 11521
stor4(i)=0.
stor4(i)=stor4(i)+xtx((i-1)*maxpro*maxkno+j)*stor4(4*nobs+j)
stor4(i)=stor4(i)/stor4(nobs+i)
stor4(4*nobs+i)=0.
stor4(4*nobs+i)=stor4(4*nobs+i)+xtx((j-1)*maxpro*maxkno+i)*stor4(j
*)
is=0
it4=itnkno(i)-korder
if(i .ne. 1)goto 11591
is=is+1
stor2(2*maxpro*maxkno+(i-1)*maxkno+j)=stor4(4*nobs+is)
goto 11611
stor2(2*maxpro*maxkno+(i-1)*maxkno+1)=0.
is=is+1
stor2(2*maxpro*maxkno+(i-1)*maxkno+j)=stor4(4*nobs+is)
rss=0.
stor4(4*nobs+i)=0.
proj=0.
proj=proj+axis((j-1)*npred+k)*predic((i-1)*npred+k)
stor4(4*nobs+i)=stor4(4*nobs+i)+bvalue(stor2(maxkno*maxpro+(j-1)*m
*axkno+1),stor2(2*maxpro*maxkno+(j-1)*maxkno+1), itnkno(j)-korder,
*korder,proj,0)
stor4(4*nobs+i)=respon(i)-stor4(4*nobs+i)
rss=rss+stor4(4*nobs+i)**2
asr=rss/float(nobs-nbspli)
return
*le to compute eigenvalue nr ,i3, 7h of xtx)
end
subroutine knopla(ipro)
real*8 sxy,sx,sy,sxx,syy,a,b,davsr
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
proj=0.
itag(i)=i
proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j)
stor4(i)=proj
call sort(stor4(1),itag,1,nobs)
i=0
i=i+1
m0=i
i=i+1
if(i.ge.nobs)goto 11692
goto 11691
if(i.eq.m0)goto 11681
call ranums(stor4(2*nobs+1),i-m0)
goto 11703
k=j+stor4(2*nobs+j-m0+1)*(i-j+1)
xt=stor4(j)
it=itag(j)
stor4(j)=stor4(k)
itag(j)=itag(k)
stor4(k)=xt
itag(k)=it
goto 11701
goto 11681
avvar=0.
ibw2=sambw*bandw*nobs/2.0
sxy=0.0
sx=sxy
sy=sx
sxx=sy
syy=sxx
x=stor4(j)
y=stor4(4*nobs+itag(j))
sxy=sxy+x*y
sx=sx+x
sy=sy+y
sxx=sxx+x*x
syy=syy+y*y
ibweff=ibw2
if(i .gt. nobs-ibw2)goto 11741
x=stor4(i+ibw2)
y=stor4(4*nobs+itag(i+ibw2))
ibweff=ibweff+1
sxy=sxy+x*y
sx=sx+x
sy=sy+y
sxx=sxx+x*x
syy=syy+y*y
if(i .le. ibw2+1)goto 11761
x=stor4(i-ibw2-1)
y=stor4(4*nobs+itag(i-ibw2-1))
ibweff=ibweff-1
sxy=sxy-x*y
sx=sx-x
sy=sy-y
sxx=sxx-x*x
syy=syy-y*y
denom=sxx-sx*sx/ibweff
a=0.0
if(denom.gt.0.0) a=(sxy-sx*sy/ibweff)/denom
b=(sy-a*sx)/ibweff
var=(syy+a*a*sxx+2.0*(a*(b*sx-sxy)-b*sy))/ibweff +b*b
stor4(2*nobs+i)=var
ibw2=varbw*bandw*nobs/2.0
sx=0.0
sx=sx+stor4(2*nobs+j)
ibweff=ibw2
if(i .gt. nobs-ibw2)goto 11801
sx=sx+stor4(2*nobs+i+ibw2)
ibweff=ibweff+1
if(i .le. ibw2+1)goto 11821
sx=sx-stor4(2*nobs+i-ibw2-1)
ibweff=ibweff-1
stor4(nobs+i)=sx/ibweff
avvar=0.
avvar=avvar+stor4(nobs+i)
avvar=avvar/float(nobs)
rhigh=avvar*smofac(3)
rlow=avvar*smofac(1)
avvar=0.
if(stor4(nobs+i) .le. rhigh)goto 11861
stor4(nobs+i)=rhigh
goto 11851
stor4(nobs+i)=rlow
stor4(nobs+i)=1./stor4(nobs+i)
avvar=avvar+stor4(nobs+i)
stor4(nobs+i)=stor4(nobs+i)/avvar
nkint=maxkno-2*korder
stor2(maxkno*maxpro+(ipro-1)*maxkno+i)=stor4(1)
nk=korder
nsame=korder
sum=0.
ind=0
ind=ind+1
sum=sum+stor4(nobs+ind)
if(sum.ge.float(i)/float(nkint+1))goto 11912
goto 11911
if(stor4(ind) .ne. stor2(maxkno*maxpro+(ipro-1)*maxkno+nk))goto 11
*931
nsame=nsame+1
goto 11941
nsame=1
if(nsame .le. korder)goto 11961
goto 11901
nk=nk+1
stor2(maxkno*maxpro+(ipro-1)*maxkno+nk)=stor4(ind)
goto 11973
if(stor2(maxkno*maxpro+(ipro-1)*maxkno+i) .eq. stor4(nobs))goto 11
*991
goto 11972
nk=nk-1
goto 11971
if(nk .ge. korder)goto 12011
write(itap1,10090)(axis((ipro-1)*npred+i),i=1,npred)
stop
d=(stor4(nobs)-stor4(1))*10.e-6
goto 12023
stor2(maxkno*maxpro+(ipro-1)*maxkno+i)=stor4(nobs)+d
goto 12021
itnkno(ipro)=nk+korder
return
*xis:/20(7(2x,f7.4)/))
end
subroutine trimin(x,nvar,fmin,fun,maxit,iprint,scrat)
dimension x(nvar),scrat(nvar,3)
common/ntapes/itap1,itap2,itap3,itap4
external fun
if(iprint .le. 0)goto 12041
write(itap2,10090)(x(i),i=1,nvar)
nit=1
nfun=0
xnorm=0.
goto 12063
xnorm=xnorm+x(i)**2
goto 12061
xnorm=sqrt(xnorm)
scrat(1,1)=sqrt((1.-x(1)/xnorm)/2.)
if(scrat(1,1) .ne. 0.)goto 12081
goto 12093
scrat(i,1)=0.
goto 12091
goto 12101
goto 12113
scrat(i,1)=-x(i)/(2.*xnorm*scrat(1,1))
goto 12111
goto 12123
goto 12133
scrat(j,2)=-2.*scrat(i,1)*scrat(j,1)
goto 12131
scrat(i,2)=scrat(i,2)+1
call cirmin(x,scrat(1,2),nvar,fun,nfun,scrat(1,3))
goto 12121
if(mod(nit,iprint) .ne. 0 .or. iprint .le. 0)goto 12151
fmin=fun(x,nvar)
write(itap2,10190)nit,nfun,fmin,(x(i),i=1,nvar)
nit=nit+1
if(nit .le. maxit)goto 12171
goto 12052
goto 12051
nit=nit-1
fmin=fun(x,nvar)
nfun=nfun+1
if(mod(nit,iprint) .eq. 0 .or. iprint .le. 0)goto 12191
write(itap2,10190)nit,nfun,fmin,(x(i),i=1,nvar)
if(iprint .le. 0)goto 12211
write(itap2,10320)
return
* 15h current axis =/10(5(2x,g12.6)/))
end
subroutine cirmin(x,rotvec,nvar,fun,nfun,scrat)
dimension x(nvar),rotvec(nvar),scrat(nvar),delta(6)
data delta/0.05,0.1,0.2,0.2,0.2,0.3/
data pi/3.14159/
idelta=0
u=0.
v=0.
w=0.
sign=1.
fu=fun(x,nvar)
f0=fu
nfun=nfun+1
idelta=idelta+1
v=u+delta(idelta)
acos=cos(pi*v)
asin=sin(pi*v)
goto 12223
scrat(i)=acos*x(i)+asin*rotvec(i)
goto 12221
fv=fun(scrat,nvar)
nfun=nfun+1
if(fv .lt. fu)goto 12241
sign=-1
idelta=0
uu=u
u=-v
v=uu
fuu=fu
fu=fv
fv=fuu
idelta=idelta+1
w=v+delta(idelta)
if(w .ge. 1.)goto 12271
sw=sign*w
acos=cos(sw*pi)
asin=sin(sw*pi)
goto 12283
scrat(i)=acos*x(i)+asin*rotvec(i)
goto 12281
fw=fun(scrat,nvar)
nfun=nfun+1
if(fw .le. fv)goto 12301
goto 12310
u=v
fu=fv
v=w
fv=fw
goto 12321
w=1.
fw=f0
goto 12310
goto 12251
a2=(fw-fu)/((w-u)*(w-v))-(fv-fu)/((v-u)*(w-v))
a1=(fv-fu)/(v-u)-a2*(v+u)
xmin=-a1/(2*a2)*sign
acos=cos(pi*xmin)
asin=sin(pi*xmin)
goto 12333
scrat(i)=acos*x(i)+asin*rotvec(i)
goto 12331
fsol=fun(scrat,nvar)
nfun=nfun+1
if(fsol .gt. fv)goto 12351
goto 12363
x(i)=scrat(i)
goto 12361
return
acos=cos(sign*pi*v)
asin=sin(sign*pi*v)
goto 12373
x(i)=acos*x(i)+asin*rotvec(i)
goto 12371
return
end
subroutine plores
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
common/iplopr/iplopr
iplopr=0
write(itap3,10190)(axis((npro-1)*npred+j),j=1,npred)
ntrim=plotrm*nobs
imin=ntrim+1
imax=nobs-ntrim
s=0.
s=s+axis((npro-1)*npred+j)*predic((i-1)*npred+j)
stor4(i)=s
itag(i)=i
call sort(stor4(1),itag,1,nobs)
xmin=stor4(imin)
xmax=stor4(imax)
xmax=xmax+1.e-4*(xmax-xmin)
rmin=rbig
rmax=-rmin
r=stor4(3*nobs+itag(i))
if(r .ge. rmin)goto 12421
rmin=r
if(r .le. rmax)goto 12441
rmax=r
stor4(nobs+i)=r-stor4(4*nobs+itag(i))
rmax=rmax+1.e-4*(rmax-rmin)
call psetup(min0(nxplt,100),xmin,xmax, min0(nyplt,50),rmin,rmax)
call plot(stor4(i),stor4(3*nobs+itag(i)))
call pprint
write(itap3,10090)
*38h moving average smooth along direction/ 20(7(2x,f7.4)/))
return
end
function ufun1(x,y)
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
common/iplopr/iplopr
if(iplopr .ne. 0)goto 12471
if(x .gt. stor4(1))goto 12491
ufun1=y-stor4(nobs+1)
return
if(x .lt. stor4(nobs))goto 12511
ufun1=y-stor4(nobs+nobs)
return
min=1
max=nobs
nt=(max+min)/2
if(x .ge. stor4(nt))goto 12541
max=nt
goto 12551
min=nt
goto 12521
if(stor4(max) .ne. stor4(min))goto 12571
ufun1=y-stor4(nobs+min)
return
ufun1=y-stor4(nobs+min)-(stor4(nobs+max)-stor4(nobs+min)) *(x-sto
*r4(min))/(stor4(max)-stor4(min))
goto 12581
xx=x
if(xx .ge. xlim(1,iplopr))goto 12601
xx=xlim(1,iplopr)
if(xx .le. xlim(2,iplopr))goto 12621
xx=xlim(2,iplopr)
r=bvalue(rkno((iplopr-1)*maxkno+1),coeff((iplopr-1)*maxkno+1), nk
*no(iplopr)-korder,korder,xx,0)
ufun1=r-y
return
end
function ufun2(x,y)
ufun2=1.
return
end
function smuuth(ax,np)
dimension ax(np)
dimension ibw(3),ibwt(3)
real*8 sxy(3),sx(3),sy(3),sxx(3),syy,a,b,davsr
common/inprm/npred,nobs
common/iprint/iprint
common/minpar/nprint,maxit
common/ntapes/itap1,itap2,itap3,itap4
common/ntry/ntry
common/npro/npro
common/terprm/maxtry,maxpro,ppconv
common/smothr/korder,maxkno,banfac
common/smopar/bandw,sambw,varbw,smofac(3)
common/pltprm/nxplt,nyplt,plotrm
common/predic/predic(1)
common/respon/respon(1)
common/stor1/xlim(2,1)
common/axis/axis(1)
common/stor2/stor2(1)
common/stor3/xtx(1)
common/stor4/stor4(1)
common/stor5/itnkno(1)
common/stor6/itag(1)
common/asrold/asrold
common/pursut/mode
common/irest/irest
common/consts/rbig,reps,ratio
common/coeff/coeff(1)
common/rkno/rkno(1)
common/nkno/nkno(1)
common/inmode/inmode
ic=0
if(ax(i).ne.0.0) ic=ic+1
if(ic .ne. 1)goto 12651
goto 12663
if(ax(id).ne.0.0)goto 12662
goto 12661
goto 12673
if(ic .ne. 1)goto 12691
stor4(nobs+i)=predic((i-1)*npred+id)
goto 12671
sum=0.
goto 12703
sum=sum+ax(j)*predic((i-1)*npred+j)
goto 12701
stor4(nobs+i)=sum
goto 12671
goto 12713
itag(i)=i
goto 12711
call sort(stor4(nobs+1),itag,1,nobs)
i=0
i=i+1
m0=i
i=i+1
if(i.ge.nobs)goto 12732
goto 12731
if(i.eq.m0)goto 12721
call ranums(stor4(2*nobs+1),i-m0)
goto 12743
k=j+stor4(2*nobs+j-m0+1)*(i-j+1)
xt=stor4(nobs+j)
it=itag(j)
stor4(nobs+j)=stor4(nobs+k)
itag(j)=itag(k)
stor4(nobs+k)=xt
itag(k)=it
goto 12741
goto 12721
avvar=0.
ibw2=sambw*bandw*nobs/2.0
sxy(1)=0.0
sx(1)=sxy(1)
sy(1)=sx(1)
sxx(1)=sy(1)
syy=sxx(1)
x=stor4(nobs+j)
y=stor4(3*nobs+itag(j))
sxy(1)=sxy(1)+x*y
sx(1)=sx(1)+x
sy(1)=sy(1)+y
sxx(1)=sxx(1)+x*x
syy=syy+y*y
ibweff=ibw2
if(i .gt. nobs-ibw2)goto 12781
x=stor4(nobs+i+ibw2)
y=stor4(3*nobs+itag(i+ibw2))
ibweff=ibweff+1
sxy(1)=sxy(1)+x*y
sx(1)=sx(1)+x
sy(1)=sy(1)+y
sxx(1)=sxx(1)+x*x
syy=syy+y*y
if(i .le. ibw2+1)goto 12801
x=stor4(nobs+i-ibw2-1)
y=stor4(3*nobs+itag(i-ibw2-1))
ibweff=ibweff-1
sxy(1)=sxy(1)-x*y
sx(1)=sx(1)-x
sy(1)=sy(1)-y
sxx(1)=sxx(1)-x*x
syy=syy-y*y
denom=sxx(1)-sx(1)*sx(1)/ibweff
a=0.0
if(denom.gt.0.0) a=(sxy(1)-sx(1)*sy(1)/ibweff)/denom
b=(sy(1)-a*sx(1))/ibweff
var=(syy+a*a*sxx(1)+2.0*(a*(b*sx(1)-sxy(1))-b*sy(1)))/ibweff +b*b
avvar=avvar+var
stor4(4*nobs+i)=var
avvar=avvar/float(nobs)
ibw2=varbw*bandw*nobs/2.0
sx(1)=0.0
sx(1)=sx(1)+stor4(4*nobs+j)
ibweff=ibw2
if(i .gt. nobs-ibw2)goto 12841
sx(1)=sx(1)+stor4(4*nobs+i+ibw2)
ibweff=ibweff+1
if(i .le. ibw2+1)goto 12861
sx(1)=sx(1)-stor4(4*nobs+i-ibw2-1)
ibweff=ibweff-1
stor4(2*nobs+i)=sx(1)/ibweff
sxy(i)=0.0
sx(i)=sxy(i)
sy(i)=sx(i)
sxx(i)=sy(i)
ibw(i)=smofac(i)*bandw*nobs/2.0+0.5
goto 12883
x=stor4(nobs+j)
y=stor4(3*nobs+itag(j))
sxy(3)=sxy(3)+x*y
sx(3)=sx(3)+x
sy(3)=sy(3)+y
sxx(3)=sxx(3)+x*x
if(j.gt.ibw(2))goto 12881
sx(2)=sx(2)+x
sy(2)=sy(2)+y
sxx(2)=sxx(2)+x*x
sxy(2)=sxy(2)+x*y
if(j.gt.ibw(1))goto 12881
sx(1)=sx(1)+x
sy(1)=sy(1)+y
sxx(1)=sxx(1)+x*x
sxy(1)=sxy(1)+x*y
goto 12881
ibwt(i)=ibw(i)-1
if(i .le. 1)goto 12931
x=stor4(nobs+i-1)
y=stor4(3*nobs+itag(i-1))
sxy(j)=sxy(j)+x*y
sx(j)=sx(j)+x
sy(j)=sy(j)+y
sxx(j)=sxx(j)+x*x
x=stor4(nobs+i)
y=stor4(3*nobs+itag(i))
sxy(j)=sxy(j)-x*y
sx(j)=sx(j)-x
sy(j)=sy(j)-y
sxx(j)=sxx(j)-x*x
if(i .gt. nobs-ibw(j))goto 12951
x=stor4(nobs+i+ibw(j))
y=stor4(3*nobs+itag(i+ibw(j)))
ibwt(j)=ibwt(j)+1
sxy(j)=sxy(j)+x*y
sx(j)=sx(j)+x
sy(j)=sy(j)+y
sxx(j)=sxx(j)+x*x
if(i .le. ibw(j)+1)goto 12971
x=stor4(nobs+i-ibw(j)-1)
y=stor4(3*nobs+itag(i-ibw(j)-1))
ibwt(j)=ibwt(j)-1
sxy(j)=sxy(j)-x*y
sx(j)=sx(j)-x
sy(j)=sy(j)-y
sxx(j)=sxx(j)-x*x
thisfc=stor4(2*nobs+i)/avvar
if(thisfc .gt. smofac(2))goto 12991
denom=sxx(1)-sx(1)*sx(1)/ibwt(1)
a=0.0
if(denom.gt.0.0) a=(sxy(1)-sx(1)*sy(1)/ibwt(1))/denom
s1=a*stor4(nobs+i)+(sy(1)-a*sx(1))/ibwt(1)
goto 13001
denom=sxx(3)-sx(3)*sx(3)/ibwt(3)
a=0.0
if(denom.gt.0.0) a=(sxy(3)-sx(3)*sy(3)/ibwt(3))/denom
s3=a*stor4(nobs+i)+(sy(3)-a*sx(3))/ibwt(3)
if(thisfc .gt. smofac(1))goto 13021
stor4(4*nobs+i)=s1
goto 12901
if(thisfc .lt. smofac(3))goto 13041
stor4(4*nobs+i)=s3
goto 12901
denom=sxx(2)-sx(2)*sx(2)/ibwt(2)
a=0.0
if(denom.gt.0.0) a=(sxy(2)-sx(2)*sy(2)/ibwt(2))/denom
s2=a*stor4(nobs+i)+(sy(2)-a*sx(2))/ibwt(2)
if(thisfc .gt. smofac(2))goto 13061
wt=(smofac(2)-thisfc)/(smofac(2)-smofac(1))
stor4(4*nobs+i)=(1.0-wt)*s2+wt*s1
goto 13071
wt=(thisfc-smofac(2))/(smofac(3)-smofac(2))
stor4(4*nobs+i)=(1.0-wt)*s2+wt*s3
stor4(nobs+itag(i))=stor4(4*nobs+i)
avsr=0.
r=stor4(3*nobs+i)-stor4(nobs+i)
avsr=avsr+r**2
stor4(4*nobs+i)=r
avsr=avsr/float(nobs)
smuuth=avsr
return
end
SUBROUTINE PSETUP(NNXIN,XLIN,XHIN,NNYIN,YLIN,YHIN)
C INITIALIZATION FOR PLOTTING ROUTINES.
C NNXIN,NNYIN = NUMBER OF X OR Y BINS
C XLIN = LEFT LIMIT FOR X
C XHIN = RIGHT LIMIT FOR X
C YLIN = BOTTOM LIMIT FOR Y
C YHIN = TOP LIMIT FOR Y
C (THERE IS NO REQUIREMENT FOR THE SIGN OF THE BINWIDTH,
C BUT IT MUST BE NONZERO, I.E. THE LIMITS MUST NOT BE EQUAL)
C MAXIMUM NUMBER OF BINS IS GIVEN BY MAXX, MAXY, AND MAXXY, BELOW
C SET UP THE PLOTS
COMMON/PLTCM1/ STATS(3,3),XPROJ(100),YPROJ(100),ARRAY(5000)
COMMON /PLOTCM/ NNX,XL,XH,DX,NNY,YL,YH,DY
EQUIVALENCE (IVAL,XVAL)
C
C GET THE RIGHT TYPE FOR INPUT VALUES
C NNX IS INTEGER
COMMON/NTAPES/ITAP1,ITAP2,ITAP3,ITAP4
DATA LX/1HX/,LY/1HY/
C MAXX IS RESTRICTED BY FORMAT STATEMENTS AND
C TO WIDTH OF PAPER. MAXY AND MAXXY ARE LIMITED
C ONLY BY THE DIMENSION STATEMENTS IN COMMON.
DATA MAXX/100/,MAXY/100/,MAXXY/5000/
IVAL = NNXIN
* This is not portable Fortran, and apparently is needed by masa anyway.
* IF (XVAL.EQ.0.0 .OR. IABS(IVAL).GT.250 000 000) IVAL = XVAL
NNX = IVAL
C XL IS REAL
XVAL = XLIN
* IF (XVAL.NE.0.0 .AND. IABS(IVAL).LE.250 000 000) XVAL = IVAL
XL = XVAL
C XH IS REAL
XVAL = XHIN
* IF (XVAL.NE.0.0 .AND. IABS(IVAL).LE.250 000 000) XVAL = IVAL
XH = XVAL
C NNY IS INTEGER
IVAL = NNYIN
* IF (XVAL.EQ.0.0 .OR. IABS(IVAL).GT.250 000 000) IVAL = XVAL
NNY = IVAL
C YL IS REAL
XVAL = YLIN
* IF (XVAL.NE.0.0 .AND. IABS(IVAL).LE.250 000 000) XVAL = IVAL
YL = XVAL
C YH IS REAL
XVAL = YHIN
* IF (XVAL.NE.0.0 .AND. IABS(IVAL).LE.250 000 000) XVAL = IVAL
YH = XVAL
C CHECK VALUES FOR OBVIOUS ERROR
IF (XL.EQ.XH) GO TO 55
IF (YL.EQ.YH) GO TO 60
IF (NNX.LE.0 .OR. NNX.GT.MAXX) GO TO 40
IF (NNY.LE.0 .OR. NNY.GT.MAXY) GO TO 40
NNXY = NNX*NNY
IF (NNXY.GT.MAXXY) GO TO 40
C THEY LOOK OK
DX = (XH-XL)/FLOAT(NNX)
DY = (YH-YL)/FLOAT(NNY)
DO 10 I=1,NNXY
10 ARRAY(I)=0.
DO 21 J=1,3
DO 20 I=1,3
20 STATS(I,J)=0.
21 CONTINUE
DO 22 I=1,NNX
22 XPROJ(I)=0.
DO 24 I=1,NNY
24 YPROJ(I)=0.
RETURN
C ---------- ERRORS FROM HERE ON -------------
40 WRITE(ITAP3,50) NNX,NNY,MAXX,MAXY,MAXXY
50 FORMAT (' THE PLOTTING ROUTINE WAS GIVEN NNX=',I4,/
X ' AND NNY=',I4,'. IT REQUIRES NNX FROM 1 TO',I4/
X ' NNY FROM 1 TO',I4,' AND NNX*NNY.LE.',I5)
STOP
55 LY = LX
YL = XL
60 WRITE(ITAP3,70) LY,YL
70 FORMAT (' THE PLOTTING ROUTINE WAS GIVEN EQUAL UPPER ',
X ' AND LOWER LIMITS FOR ',A1,/
X ' OF',G12.5,'. ZERO BINWIDTH IS ILLEGAL.')
STOP
END
SUBROUTINE PLOT(XVAL,YVAL)
COMMON/PLTCM1/ STATS(3,3),XPROJ(100),YPROJ(100),ARRAY(5000)
COMMON /PLOTCM/ NNX,XL,XH,DX,NNY,YL,YH,DY
EQUIVALENCE (X,IXX),(Y,IYY)
DATA WEIGHT /1.0/
C
C CHECK FOR INTEGER INPUT
X=XVAL
IF (X.NE.0. .AND. IABS(IXX).LE.250 000 000)
X X=FLOAT(IXX)
Y=YVAL
IF (Y.NE.0. .AND. IABS(IYY).LE.250 000 000)
X Y=FLOAT(IYY)
IX = (X-XL)/DX
IY = (Y-YL)/DY
C FIND AREA (1-2-3 / 4-5-6 / 7-8-9)
IXI = 2
IYI = 2
IF (IY .GE. NNY) IYI = 1
IF (IY .LT. 0) IYI = 3
IF (IX .GE. NNX) IXI = 3
IF (IX .LT. 0) IXI = 1
C INCREMENT PLOT STATISTICS
STATS(IXI,IYI) = STATS(IXI,IYI) + WEIGHT
C QUIT IF NOT IN BOUNDS BOTH WAYS
IF (IXI.NE.2 .OR. IYI.NE.2) R E T U R N
C BUMP PLOT CONTENTS
NN = IX + NNX*IY
ARRAY(NN+1) = ARRAY(NN+1) + WEIGHT
XPROJ(IX+1) = XPROJ(IX+1) + WEIGHT
YPROJ(IY+1) = YPROJ(IY+1) + WEIGHT
RETURN
END
SUBROUTINE PPRINT
COMMON/PLTCM1/ STATS(3,3),XPROJ(100),YPROJ(100),ARRAY(5000)
COMMON /PLOTCM/ NNX,XL,XH,DX,NNY,YL,YH,DY
COMMON/NTAPES/ITAP1,ITAP2,ITAP3,ITAP4
C
INTEGER LINE1(116),LINE2(116),LINE3(116)
REAL WD(120)
EQUIVALENCE (WD(1),IWD(1))
C CHARACTERS, FORMAT A4
DIMENSION BINED(4),IPROJX(5)
C CHARACTERS
INTEGER CHAR(35),ICH1(15),ICH2(15),ICH3(15),ICHRAY(15)
C SIZE OF ICH ARRAYS IS ALSO IN NCHAR, BELOW
INTEGER BLANK,SPLUS,SMINUS,ASTER,EQUALS,EYE,BINED,SLASH
INTEGER CRVCHR,CRVCH2,IWD(116),IPROJY(4,3)
C
DATA IOVER /0/
DATA BLANK,SPLUS,SMINUS,ASTER /1H ,1H+,1H-,1H*/
DATA EQUALS,EYE,SLASH /1H=,1HI,1H//
C CHARACTERS FOR SHADED PLOTS
DATA NCHAR /15/
DATA ICH1 /
D 1H=,1H+,1HH,1H*,1HW,1HO,1HH,1HO,1HO,1HO,1HH,1HO,1HO,1HO,1HO/
DATA ICH2 /
D 1H ,1H ,1H ,1H ,1H ,1H-,1H=,1H=,1H*,1HI,1HI,1HX,1HX,1HX,1HX/
DATA ICH3 /
D 1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H-,1H=,1H*,1HI/
C CHARACTERS FOR NORMAL PLOTS
DATA CHAR /1H+,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,
X 1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,
X 1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ/
C
C MISC CHARACTERS
DATA IPROJY /1HP,1HR,1HO,1HJ,1HO,1HN,1H ,1HY,1HA,1HX,1HI,1HS/
C FORMAT A4
DATA IPROJX /1H ,1H ,4HPROJ,4HON X,4HAXIS/
DATA BINED /4HLEFT,4HBIN ,4HEDGE,1H /
DATA HNONE /4HNONE/
C
CRVCHR = ASTER
CRVCH2 = EQUALS
C CHECK FLAG FOR SHADED PLOTTING
IF (IOVER.EQ.0) GO TO 15
C SHADED PLOTS. FIND MAX BIN CONTENTS
CRVCHR = SMINUS
NXY = NNX*NNY
VALMAX = NCHAR
DO 5 IXY=1,NXY
VALMAX = AMAX1(VALMAX,ARRAY(IXY))
5 CONTINUE
C SET UP ARRAY
N = NCHAR
LAST = VALMAX + 1.0
8 INC = (LAST-2)/N + 1
LAST = LAST - INC
ICHRAY(N) = LAST
N = N - 1
IF (N.NE.0) GO TO 8
C --------- DO TITLES ---------
C FIRST LINE
15 WRITE(ITAP3,20)
20 FORMAT (1H )
C PLOT STATISTICS
WRITE(ITAP3,30) STATS
30 FORMAT( 18H PLOT STATISTICS =,3F5.0,3X,1H/,3F5.0,3X,1H/,3F5.0)
40 CONTINUE
IF (STATS(2,2).EQ.0.) GO TO 200
C -------------- PRINT CONTENTS OF X CHANNELS ------
DO 50 I=1,116
IWD(I) = BLANK
50 LINE1(I) = BLANK
L = NNX + 2
DO 100 J=1,5
CALL VPRINT(XPROJ,LINE1,NNX,2,6,0,J,JPRT)
IF (J.LT.3) GO TO 95
DO 90 I=1,4
90 LINE1(L+I) = IPROJY(I,J-2)
GO TO 98
95 IF (JPRT.EQ.0) GO TO 100
98 WRITE(ITAP3,310) IPROJX(J),LINE1
100 CONTINUE
C PUT ON PLOT BOUNDARY
DO 105 I=NNX,116
LINE3(I) = BLANK
LINE2(I) = BLANK
105 LINE1(I) = BLANK
DO 110 I=1,NNX
110 IWD(I) = SMINUS
DO 111 I=1,NNX,10
111 IWD(I) = SPLUS
C IWD CONTAINS CHARACTERS. IT WILL BE USED
C AGAIN FOR LOWER AXIS(SO DONT USE IWD OR WD)
WRITE(ITAP3,330) YH,BLANK,IWD
C --------- MAKE THE PLOT ----------
UFLAG = UFUN1(XL,YL)
VFLAG = UFUN2(XL,YL)
C J COUNTS LINES, IY COUNTS FROM NNY TO 1
IY = NNY + 1
YSTEPS = FLOAT(NNY)
Y1 = YH
DO 140 J=1,NNY
IY = IY - 1
LL = NNX*(NNY-J)
Y2 = Y1
YSTEPS = YSTEPS - 1.
Y1 = YL + YSTEPS*DY
LOVER2 = 0
LOVER3 = 0
C LOOP OVER X-BINS
XSTEPS = 0.
X2 = XL
C SET BIN COUNTER FOR USER CURVE
IUX = 0
IVX = 0
DO 130 I=1,NNX
XSTEPS = XSTEPS + 1.
X1 = X2
X2 = XL + XSTEPS*DX
LL = LL + 1
NMBR=ARRAY(LL)
IF (NMBR.EQ.0) GO TO 120
IF (IOVER.NE.0) GO TO 114
NMBR = MIN0(NMBR,35)
LINE1(I)=CHAR(NMBR)
GO TO 122
114 DO 116 K=1,NCHAR
IF (NMBR.LE.ICHRAY(K)) GO TO 118
116 CONTINUE
K = NCHAR
118 LINE1(I) = ICH1(K)
LINE2(I) = ICH2(K)
IF (LINE2(I).NE.BLANK) LOVER2 = 1
LINE3(I) = ICH3(K)
IF (LINE3(I).NE.BLANK) LOVER3 = 1
GO TO 122
120 LINE1(I) = BLANK
LINE2(I) = BLANK
LINE3(I) = BLANK
122 IF (UFLAG.EQ.HNONE) GO TO 127
IF (IUX.EQ.I) GO TO 125
C LAST BIN WASNT TREATED. GET U-VALUES
U1 = SIGN(1.,UFUN1(X1,Y1))
U2 = SIGN(1.,UFUN1(X1,Y2))
125 U1N = SIGN(1.,UFUN1(X2,Y1))
U2N = SIGN(1.,UFUN1(X2,Y2))
IF (U1N.EQ.U2 .AND. U2N.EQ.U1) GO TO 126
LINE1(I) = CRVCHR
LINE2(I) = BLANK
LINE3(I) = BLANK
C PLACE RIGHT-SIDE VALUES TO BE
C LEFT-SIDE ON NEXT TIME THRU LOOP
126 U1 = U1N
U2 = U2N
C AND SET BIN COUNTER FOR NEXT BIN
IUX = I+1
127 IF (VFLAG.EQ.HNONE) GO TO 130
IF (IVX.EQ.I) GO TO 128
C LAST BIN WASNT TREATED. GET U-VALUES
V1 = SIGN(1.,UFUN2(X1,Y1))
V2 = SIGN(1.,UFUN2(X1,Y2))
128 V1N = SIGN(1.,UFUN2(X2,Y1))
V2N = SIGN(1.,UFUN2(X2,Y2))
IF (V1N.EQ.V2 .AND. V2N.EQ.V1) GO TO 129
LINE1(I) = CRVCH2
LINE2(I) = BLANK
LINE3(I) = BLANK
C PLACE RIGHT-SIDE VALUES TO BE
C LEFT-SIDE ON NEXT TIME THRU LOOP
129 V1 = V1N
V2 = V2N
C AND SET BIN COUNTER FOR NEXT BIN
IVX = I+1
130 CONTINUE
C CHARACTER FOR RIGHT+LEFT BORDER
IAXIS=EYE
IF(MOD(IY,10).EQ.1) IAXIS=SPLUS
LINE1(NNX+1) = IAXIS
C ENCODE VALUE FOR Y-AXIS PROJECTION
L = NNX + 2
DO 132 I=1,5
CALL VPRINT(YPROJ(IY),LINE1(L),1,2,6,0,I,JPRT)
132 L=L+1
IF (IOVER.EQ.0) GO TO 138
C PUT ON SYMBOL DICTIONARY
IF (J - (NCHAR+1)) 136,133,138
C (NO MORE LINES. MAKE SURE SPACES ARE BLANK)
133 DO 134 I=L,116
LINE1(I) = BLANK
LINE2(I) = BLANK
134 LINE3(I) = BLANK
GO TO 138
C PUT IN VALUE AND SYMBOL FOR DICTIONARY
136 L = L + 3
DO 137 I=1,3
CALL VPRINT(ICHRAY(J),LINE1(L),1,1,3,0,I,JPRT)
137 L = L + 1
LINE1(L) = SLASH
LINE1(L+1) = ICH1(J)
LINE2(L+1) = ICH2(J)
IF (ICH2(J).NE.BLANK) LOVER2 = 1
LINE3(L+1) = ICH3(J)
IF (ICH3(J).NE.BLANK) LOVER3 = 1
C PUT OUT A LINE
138 WRITE(ITAP3,330) Y1,IAXIS,LINE1
C AND ANY OVERPRINTING
IF (LOVER2.NE.0) WRITE(ITAP3,335) LINE2
IF (LOVER3.NE.0) WRITE(ITAP3,335) LINE3
C LOOP FOR NEXT LINE
140 CONTINUE
C PUT ON LOWER AXIS
WRITE(ITAP3,310) BLANK,IWD
C PUT ON X-BIN EDGE VALUES
DO 150 I=NNX,116
150 LINE1(I) = BLANK
X2 = 0.
DO 160 I=1,NNX
WD(I)=XL + X2*DX
160 X2 = X2 + 1.
I = 0
DO 170 J=1,10
CALL VPRINT(WD,LINE1,NNX,2,10,4,J,JPRT)
IF (JPRT.NE.0) GO TO 165
C NO CHARACTERS THIS LINE. LEAVE DO-LOOP IF THRU
IF (I.NE.0) GO TO 200
GO TO 170
165 IF (I.LT.4) I=I+1
WRITE(ITAP3,310) BINED(I),LINE1
170 CONTINUE
200 R E T U R N
C
310 FORMAT (3X,A4,6X,116A1)
330 FORMAT(1X,F10.4, 1X,A1,116A1)
335 FORMAT (1H+,12X,116A1)
END
SUBROUTINE VPRINT(VIN,VOUT,NB,MODE,JW,JD,JENT,JPRINT)
C VIN = ARRAY OF NB INPUT NUMBERS
C VOUT = ARRAY OF NB OUTPUT CHARACTERS
C NB = NUMBER OF COLUMNS (BINS) TO DECODE.
C MODE 1 = INTEGER INPUT AND OUTPUT
C 2 = REAL IN, F OUT
C 3 = REAL IN, E OUT
C
C JW = PLACES
C JD = DECIMAL PLACES (F AND E FORMAT)
C JENT = THIS PLACE (VPRINT SHOULD BE CALLED
C JW TIMES, WITH JENT=1,JW
C JPRINT = 0 FOR ALL BLANK LINE, 1 OTHERWISE
DIMENSION IADD(3),TENPWR(21)
C CHARACTERS
INTEGER VOUT(NB),DIGIT(10),BLANK,SMINUS,PERIOD,SPLUS
C INPUT VALUES
REAL VIN(NB),INPXX
EQUIVALENCE(XX,IXX,INPXX)
DATA TENPWR /
X 1.E0,1.E1,1.E2,1.E3,1.E4,1.E5,1.E6,1.E7,1.E8,1.E9,1.E10
X ,1.E11,1.E12,1.E13,1.E14,1.E15,1.E16,1.E17,1.E18,1.E19,1.E20/
DATA DIGIT /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/
DATA BLANK,SPLUS,SMINUS,PERIOD / 1H ,1H+,1H-,1H./
DATA IADD/1,0,-3/
C ACCUR1=MACHINE ACCURACY. ACCUR2=0.1/ACCUR1
DATA ACCUR2/3.0E3/,ACCUR1/3.3E-5/
C SMALL IS LIMIT FOR CHECK FOR ZERO
DATA SMALL /1.E-40/
C
JPRINT=0
KD=JD
IF(MODE.EQ.1) KD=0
C FAC5 IS LOWEST SIG FIG
FAC5 = 1./TENPWR(KD+1)
C KPOW IS CURRENT PLACE (1 UNITS, 2 TENS,...)
KPOW=JW-KD-JENT+IADD(MODE)
IF (KPOW) 10,200,20
10 KPOW = KPOW + 1
C RND IS HALF OF LOWEST SIG FIG
20 RND = 0.4999*FAC5
C FAC2 IS 10. IN CURRENT DIGIT.
C FAC2 = 10.**KPOW
I=IABS(KPOW)
FAC2=1.
22 IF (I.LT.20) GO TO 25
I=I-20
FAC2=FAC2*TENPWR(21)
GO TO 22
25 FAC2=FAC2*TENPWR(I+1)
IF (KPOW.LT.0) FAC2=1./FAC2
C FAC6 IS 0.1 IN CURRENT DIGIT
C (FOR CHECKING IF NEXT IS FIRST DIGIT,
C SO THIS IS THE PLACE FOR THE MINUS SIGN)
FAC6 = 0.01*FAC2
FAC8 = FAC6 * 5.
C FAC1 BRINGS CURRENT DIGIT TO UNITS PLACE
FAC1=10.0001/FAC2
C FAC7 IS LARGEST NUMBER FOR WHICH THE
C CURRENT DIGIT IS SIGNIFICANT.
FAC7 = FAC2*ACCUR2
C SIG FIG IS IN UNITS PLACE
C NOW LOOP OVER ALL NB NUMBERS
DO 110 IB=1,NB
INPXX=VIN(IB)
IF (XX.NE.0.) GO TO 28
IF (KPOW.EQ.1) GO TO 27
VOUT(IB)=BLANK
GO TO 110
27 VOUT(IB)=DIGIT(1)
GO TO 100
28 GO TO (50,60,30),MODE
C FIGURE EXPONENTIAL FORMAT
30 IF(ABS(XX).GT.SMALL) GO TO 31
XX=0.0
ITMP=0
GO TO 32
31 TMP=ALOG10(ABS(XX))
IF(TMP.LT.0.0) TMP=TMP-1.0
ITMP = IFIX(TMP+0.000001)
32 IF(JENT-JW+2) 33,37,40
C PUT IN VALUE TIMES EXPONENTIAL BIAS
C SO I WANT XX=XX*10.0**(-ITMP)
33 XXFAC=1.
I=IABS(ITMP)
34 IF (I.LT.20) GO TO 35
I=I-20
XXFAC=XXFAC*TENPWR(21)
GO TO 34
35 XXFAC=XXFAC*TENPWR(I+1)
IF (ITMP.GT.0) XXFAC=1./XXFAC
XX=XX*XXFAC
GO TO 60
C PUT IN SIGN OF EXP
37 VOUT(IB)=SPLUS
IF(ITMP.LT.0) VOUT(IB)=SMINUS
GO TO 100
C PUT IN EXPONENTIAL
40 ITMP = IABS(ITMP)
IF (JENT.EQ.JW) GO TO 42
C TENS PLACE
JX = ITMP/10
GO TO 44
C UNITS PLACE
42 JX = MOD (ITMP,10)
44 VOUT(IB) = DIGIT(JX+1)
GO TO 100
C INTEGER FORMAT. INSTALL REAL
50 XX = FLOAT(IXX)
C VALUE NOW IN XX FOR INTERPRETATION
60 X=ABS(XX)
C CHECK SIGNIFICANCE
IF (X.GT.FAC7) GO TO 68
X = X + RND
C KNOCK OFF HIGH DIGITS
XM = AMOD(X,FAC2)
IF (XM.LT.FAC8) GO TO 63
C BRING TO UNITS PLACE AND TRUNCATE
JX = XM*FAC1
C ACCEPT IF NOT A ZERO
IF(JX.NE.0) GO TO 95
C ITS A ZERO (OR BLANK)
63 VOUT(IB) = DIGIT(1)
C JUMP IF PRINTING TENS PLACE OR GREATER
C (TO CHECK FOR - AND DELETE LEADING 0)
C ACCEPT IF PRINTING UNITS
IF (KPOW-1) 65,100,70
C (ACCEPTING KPOW=0 WOULD ALLOW TENTHS PLACE)
C CHECK IF THERE ARE FURTHER NON-ZERO
C SIGNIFICANT DIGITS.
65 IF (XM.LT.FAC5) GO TO 110
IF (XM.LT.X*ACCUR1) GO TO 110
GO TO 100
C DIGIT BELOW MACHINE ACCURACY.
68 VOUT(IB) = BLANK
GO TO 110
C HUNDREDS PLACE OR GREATER
C SKIP IF HIGHER PLACE EXISTS
70 IF(X.GE.FAC2) GO TO 100
VOUT(IB)=BLANK
C CHECK FOR - SIGN
IF (XX.GE.0.0) GO TO 110
IF (KPOW.EQ.2) GO TO 90
IF (X.LT.FAC6) GO TO 110
90 VOUT(IB)=SMINUS
GO TO 100
C GET CORRESPONDING CHARACTER
95 VOUT(IB)=DIGIT(JX+1)
100 JPRINT=1
110 CONTINUE
R E T U R N
C
C KPOW=0. THIS PLACE FOR DECIMAL POINT.
200 DO 210 I=1,NB
210 VOUT(I) = PERIOD
C SET PRINT FLAG IF EXPONENTIAL FORMAT
IF (MODE.NE.2) GO TO 230
C SET PRINT FLAG IF ANY SUCCEEDING
C CHARACTERS WILL BE NON-BLANK
DO 220 I=1,NB
INPXX=VIN(I)
X = ABS(XX)
XM = AMOD(X,1.)
IF (XM.LT.FAC5) GO TO 220
IF (XM.GE.X*ACCUR1) GO TO 230
220 CONTINUE
C NO PRINT FLAG
RETURN
C SET PRINT FLAG
230 JPRINT = 1
RETURN
END
SUBROUTINE SORT (V,A,II,JJ)
C
C PUTS INTO A THE PERMUTATION VECTOR WHICH SORTS V INTO
C INCREASING ORDER. ONLY ELEMENTS FROM II TO JJ ARE CONSIDERED.
C ARRAYS IU(K) AND IL(K) PERMIT SORTING UP TO 2**(K+1)-1 ELEMENTS
C
C THIS IS A MODIFICATION OF CACM ALGORITHM #347 BY R. C. SINGLETON,
C WHICH IS A MODIFIED HOARE QUICKSORT.
C
DIMENSION A(JJ),V(1),IU(20),IL(20)
INTEGER T,TT
INTEGER A
REAL V
M=1
I=II
J=JJ
10 IF (I.GE.J) GO TO 80
20 K=I
IJ=(J+I)/2
T=A(IJ)
VT=V(IJ)
IF (V(I).LE.VT) GO TO 30
A(IJ)=A(I)
A(I)=T
T=A(IJ)
V(IJ)=V(I)
V(I)=VT
VT=V(IJ)
30 L=J
IF (V(J).GE.VT) GO TO 50
A(IJ)=A(J)
A(J)=T
T=A(IJ)
V(IJ)=V(J)
V(J)=VT
VT=V(IJ)
IF (V(I).LE.VT) GO TO 50
A(IJ)=A(I)
A(I)=T
T=A(IJ)
V(IJ)=V(I)
V(I)=VT
VT=V(IJ)
GO TO 50
40 A(L)=A(K)
A(K)=TT
V(L)=V(K)
V(K)=VTT
50 L=L-1
IF (V(L).GT.VT) GO TO 50
TT=A(L)
VTT=V(L)
60 K=K+1
IF (V(K).LT.VT) GO TO 60
IF (K.LE.L) GO TO 40
IF (L-I.LE.J-K) GO TO 70
IL(M)=I
IU(M)=L
I=K
M=M+1
GO TO 90
70 IL(M)=K
IU(M)=J
J=L
M=M+1
GO TO 90
80 M=M-1
IF (M.EQ.0) RETURN
I=IL(M)
J=IU(M)
90 IF (J-I.GE.11) GO TO 20
IF (I.EQ.II) GO TO 10
I=I-1
100 I=I+1
IF (I.EQ.J) GO TO 80
T=A(I+1)
VT=V(I+1)
IF (V(I).LE.VT) GO TO 100
K=I
110 A(K+1)=A(K)
V(K+1)=V(K)
K=K-1
IF (VT.LT.V(K)) GO TO 110
A(K+1)=T
V(K+1)=VT
GO TO 100
END
SUBROUTINE GENPTS(NVECT,NDIM,Y)
C
C
DIMENSION Y(1),U(2)
C
N=NVECT*NDIM
DO 1000 J=1,N,2
200 CONTINUE
CALL RANUMS(U,2)
V1=2.*U(1)-1.
V2=2.*U(2)-1.
S=V1*V1+V2*V2
IF(S.GE.1.) GO TO 200
S=SQRT(-2.*ALOG(S)/S)
Y(J)=V1*S
IF((J+1).LE.N)Y(J+1)=V2*S
RETURN
END
C
C
SUBROUTINE RANUMS(X,N)
REAL X(N)
C GENERATE N UNIFORM RANDOM NUMBERS IN THE INTERVAL (0.0, 1.0).
C URAND -- A PORTABLE RANDOM NUMBER GENERATOR, BY MICHAEL A.
C MALCOLM AND CLEVE B. MOLER, STANFORD COMPUTER SCIENCE REPORT
C STAN-CS-73-334, JANUARY 1974.
C NOTE THAT THE STATISTICAL PROPERTIES OF THIS ROUTINE HAVE NOT
C BEEN EXTENSIVELY VERIFIED ON A LARGE NUMBER OF ARCHITECTURES.
C THIS IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND
C SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL 2.
INTEGER IA, IC, ITWO, IY, M2, M
DOUBLE PRECISION HALFM, DATAN, DSQRT
DATA M2 / 0 /, ITWO / 2 /, IY /123456789/
IF (M2 .NE. 0) GO TO 20
C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORK LENGTH
M = 1
10 M2 = M
M = ITWO * M2
IF (M .GT. M2) GO TO 10
HALFM = M2
C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD
IA = 8 * IDINT(HALFM*DATAN(1.0D0)/8.0D0) + 5
IC = 2 * IDINT(HALFM*(0.5D0-DSQRT(3.0D0)/6.0D0)) + 1
C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT
S = 0.5 / HALFM
C NOW COMPUTE NEXT RANDOM NUMBER
20 DO 30 I=1,N
IY = IY * IA + IC
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE WORD LENGTH FOR
C ADDITION IS GREATER THAN THAT FOR MULITPLICATION
IF (IY/2 .GT. M2) IY = (IY - M2) - M2
C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER OVERFLOW
C AFFECTS THE SIGN BIT
IF (IY .LT. 0) IY = (IY + M2) + M2
30 X(I) = FLOAT(IY) * S
RETURN
END
C
C
C
SUBROUTINE SMOOTH(X,Y,N,BANDW,SORTX,SMO,R3,ITAG)
DIMENSION X(N),Y(N),SORTX(N),SMO(N),R3(N),ITAG(N)
C
C SORT OBS IN INCREASING VALUES OF X
C
DO 10 I=1,N
ITAG(I)=I
SORTX(I)=X(I)
10 CONTINUE
CALL SORT(SORTX,ITAG,1,N)
DO 20 I=1,N
SMO(I)=Y(ITAG(I))
20 CONTINUE
C
C COMPUTE RUNNING MEDIANS OF 3 AND FILL IN R3
C
N1=N-1
DO 30 I=2,N1
Y1=SMO(I-1)
Y2=SMO(I)
Y3=SMO(I+1)
IF(Y1.LE.Y2)GOTO 40
TEMP=Y1
Y1=Y2
Y2=TEMP
40 CONTINUE
RMED=Y2
IF(Y2.LE.Y3)GOTO 50
RMED=Y3
IF(Y3.GT.Y1)GOTO 50
RMED=Y1
50 CONTINUE
R3(I)=RMED
30 CONTINUE
R3(1)=SMO(1)
R3(N)=SMO(N)
C
C FIT RUNNING STRAIGHT LINES
C
IBW=N*BANDW
K=(IBW-1)/2
IF(K.LT.1)K=1
C
C INITIALIZE UPDATE
SX=0.
SY=0.
SXX=0.
SXY=0.
C
DO 60 I=1,K
SX=SX+SORTX(I)
SY=SY+R3(I)
SXX=SXX+SORTX(I)**2
SXY=SXY+SORTX(I)*R3(I)
60 CONTINUE
IBWEFF=K
C
C NOW LOOP OVER POINTS TO SMOOTH
C
DO 100 I=1,N
IF(I.GT.N-K)GOTO 70
C
C ADD OBSERVATION TO THE RIGHT
C
SX=SX+SORTX(I+K)
SY=SY+R3(I+K)
SXX=SXX+SORTX(I+K)**2
SXY=SXY+SORTX(I+K)*R3(I+K)
IBWEFF=IBWEFF+1
C
70 CONTINUE
IF(I.LE.K+1)GOTO 80
C
C DELETE OBSERVATION ON THE LEFT
C
SX=SX-SORTX(I-K-1)
SY=SY-R3(I-K-1)
SXX=SXX-SORTX(I-K-1)**2
SXY=SXY-SORTX(I-K-1)*R3(I-K-1)
IBWEFF=IBWEFF-1
C
80 CONTINUE
C
C NOW COMPUTE SLOPE AND INTERCEPT
C
B=0.
DENOM=SXX-SX**2/FLOAT(IBWEFF)
IF(DENOM.NE.0.)B=(SXY-SX*SY/FLOAT(IBWEFF))/DENOM
A=(SY-B*SX)/FLOAT(IBWEFF)
SMO(I)=A+B*SORTX(I)
100 CONTINUE
RETURN
END
.
