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 10010 do 10011 i=1,nobs s=0. 10020 do 10021 j=1,npred s=s+direct(j)*predic((i-1)*npred+j) 10021 continue 10022 continue stor4(i)=s stor4(2*nobs+i)=smodel(predic((i-1)*npred+1)) 10011 continue 10012 continue 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 10030 do 10031 i=imin,imax s=stor4(2*nobs+itag(i)) if(s .ge. rmin)goto 10051 rmin=s 10051 continue if(s .le. rmax)goto 10071 rmax=s 10071 continue 10031 continue 10032 continue rmax=rmax+1.e-4*(rmax-rmin) call psetup(min0(nxplt,100),xmin,xmax, min0(nyplt,50),rmin,rmax) 10080 do 10081 i=1,nobs call plot(stor4(i),stor4(2*nobs+itag(i))) 10081 continue 10082 continue 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) 10100 do 10101 i=1,nobs proj=0. 10110 do 10111 j=1,npred proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j) 10111 continue 10112 continue stor4(i)=proj itag(i)=i 10101 continue 10102 continue 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) 10120 do 10121 i=imin,imax r=bvalue(rkno((ipro-1)*maxkno+1),coeff((ipro-1)*maxkno+1),nkno(ipr *o)-korder, korder,stor4(i),0) stor4(nobs+i)=r 10121 continue 10122 continue rmin=rbig rmax=-rmin 10130 do 10131 i=imin,imax if(stor4(nobs+i) .ge. rmin)goto 10151 rmin=stor4(nobs+i) 10151 continue if(stor4(nobs+i) .le. rmax)goto 10171 rmax=stor4(nobs+i) 10171 continue 10131 continue 10132 continue 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) 10180 i=1 goto 10183 10181 i=i+1 10183 if((i).gt.(nkno(ipro)))goto 10182 call plot(rkno((ipro-1)*maxkno+i),rmin) goto 10181 10182 continue call pprint iplopr=0 write(itap1,10190)ipro,(axis((ipro-1)*npred+j),j=1,npred) rmin=rbig rmax=-rmin 10200 do 10201 i=imin,imax r=stor4(3*nobs+itag(i)) if(r .ge. rmin)goto 10221 rmin=r 10221 continue if(r .le. rmax)goto 10241 rmax=r 10241 continue stor4(2*nobs+i)=r 10201 continue 10202 continue 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) 10250 do 10251 i=1,nobs call plot(stor4(i),stor4(2*nobs+i)) 10251 continue 10252 continue 10260 i=1 goto 10263 10261 i=i+1 10263 if((i).gt.(nkno(ipro)))goto 10262 call plot(rkno((ipro-1)*maxkno+i),rmin) goto 10261 10262 continue call pprint return 10090 format( 46h1spline function and knots along direction nr ,i3/ 2 *0(7(2x,f7.4)/)) 10190 format( 40h1residuals and knots along direction nr ,i3/ 20(7(2x *,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 10281 continue 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 10290 format( 39h0information for restart saved on tape ,i3) 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 10311 continue irest=1 write(itap1,10290)itap4 return 10290 format( 40h0information for restart read from tape ,i3) 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 10370 do 10371 i=1,nobs stor4(3*nobs+i)=respon(i) 10371 continue 10372 continue inmode=mode if(mode .ne. 3)goto 10391 inmode=1 10391 continue sum=0. 10400 do 10401 i=1,nobs sum=sum+respon(i) 10401 continue 10402 continue sum=sum/float(nobs) sstot=0. 10410 do 10411 i=1,nobs sstot=sstot+(respon(i)-sum)**2 10411 continue 10412 continue asrold=sstot/float(nobs) write(itap1,10420)asrold goto 10431 10361 continue npro=npro+1 ntry=1 irest=0 10431 continue 10351 continue bandw=banfac/(maxkno-korder) 10440 continue 10441 continue write(itap1,10450)npro,ntry if(inmode .ne. 1)goto 10471 call select(asr) goto 10481 10471 continue call gensta call trimin(axis((npro-1)*npred+1),npred,asr,smuuth,maxit,nprint, * stor4(1)) 10481 continue 10461 continue write(itap1,10490)asr,(axis((npro-1)*npred+i),i=1,npred) call range(npro) if(iprint .lt. 2)goto 10511 call plores 10511 continue 10520 do 10521 i=1,npro call knopla(i) 10521 continue 10522 continue if(iprint .ne. 3)goto 10541 10550 do 10551 i=1,npro nki=itnkno(i) write(itap1,10560)i,nki,(stor2(maxkno*maxpro+(i-1)*maxkno+j),j=1,n *ki) 10551 continue 10552 continue 10541 continue 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 10591 continue write(itap1,10620) 10630 do 10631 i=1,nobs stor4(3*nobs+i)=stor4(4*nobs+i) 10631 continue 10632 continue 10640 do 10641 i=1,maxkno 10650 do 10651 j=1,npro 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) 10651 continue 10652 continue 10641 continue 10642 continue 10660 do 10661 i=1,npro nkno(i)=itnkno(i) 10661 continue 10662 continue npro=npro+1 ntry=1 asrold=asr 10611 continue 10581 continue call decide(iterm) if(iterm .ne. 1)goto 10681 goto 10442 10681 continue goto 10441 10442 continue npro=npro-1 if(iprint .lt. 2)goto 10701 10710 do 10711 i=1,npro call plospl(i) 10711 continue 10712 continue 10701 continue return 10090 format( 47h1multidimensional additive spline approximation 10 *h (4/19/80)/ 24h0parameters for this run) 10190 format( 13h nobs ,i5/ 13h npred ,i5) 10320 format( 13h mode ,i5/ 13h maxtry ,i5/ 13h maxp *ro ,i5/ 13h ppconv ,g12.6/ 13h maxit ,i5) 10330 format( 13h korder ,i5/ 13h maxkno ,i5/ 13h ban *fac ,g12.6) 10420 format( 43h0average squared residual around the mean ,g12.6) 10340 format( 13h iprint ,i5/ 13h nprint ,i5/ 13h plot *rm ,g12.6) 10450 format( 16h1projection no. ,i3, 12h trial no. ,i3) 10490 format( 48h0minimum asr found by search for new direction: ,g12. *6/ 20h solution direction:/20(7(2x,f7.4)/)) 10560 format( 15h0direction no. ,i3, 9h nknots=,i3/ 16h knot po *sitions:/20(7(2x,g12.6)/)) 10570 format( 43h0average squared residual from spline fit= g12.6) 10600 format( 24h0projection not accepted) 10620 format( 20h0projection accepted) 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. 10720 do 10721 i=1,npro proj=0. 10730 do 10731 j=1,npred proj=proj+x(j)*axis((i-1)*npred+j) 10731 continue 10732 continue if(proj .ge. xlim(1,i))goto 10751 proj=xlim(1,i) 10751 continue if(proj .le. xlim(2,i))goto 10771 proj=xlim(2,i) 10771 continue smod=smod+bvalue(rkno((i-1)*maxkno+1),coeff((i-1)*maxkno+1), nkno *(i)-korder,korder,proj,0) 10721 continue 10722 continue 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 10800 do 10801 ipred=1,npred 10810 do 10811 i=1,npred axis((npro-1)*npred+i)=0.0 10811 continue 10812 continue 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 10831 continue if(npro .ne. 1)goto 10851 write(itap1,10090)ipred,asr if(iprint .lt. 2)goto 10871 call range(npro) call plores 10871 continue 10851 continue 10801 continue 10802 continue 10880 do 10881 i=1,npred axis((npro-1)*npred+i)=0.0 10881 continue 10882 continue axis((npro-1)*npred+ix)=1.0 write(itap1,10190)ix,fmin goto 10891 10791 continue call genpts(1,npred,axis((npro-1)*npred+1)) sum=0. 10900 do 10901 i=1,npred sum=sum+axis((npro-1)*npred+i)**2 10901 continue 10902 continue sum=sqrt(sum) 10910 do 10911 i=1,npred axis((npro-1)*npred+i)=axis((npro-1)*npred+i)/sum 10911 continue 10912 continue write(itap1,10330) 10891 continue 10781 continue return 10090 format( 29h0average squared residual fori3, 14hth predictor= , *g12.6) 10190 format( 14h0predictor nr ,i3, 33h selected for starting direct *ion;/ 6h asr= ,g12.6) 10330 format( 26h0random starting direction) 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 10931 continue 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 10971 continue return 10951 continue if(mode .ne. 2)goto 10991 if((npro .le. maxpro) .and. (ntry .le. maxtry))goto 11011 iterm=1 11011 continue return 10991 continue 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 11051 continue if((npro .le. maxpro) .and. (ntry .le. maxtry))goto 11071 iterm=1 11071 continue 11031 continue 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 11080 do 11081 i=1,nobs proj=0. 11090 do 11091 j=1,npred proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j) 11091 continue 11092 continue if(proj .ge. xmin)goto 11111 xmin=proj 11111 continue if(proj .le. xmax)goto 11131 xmax=proj 11131 continue 11081 continue 11082 continue 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 11140 do 11141 i=1,npred if(npro .le. 1)goto 11161 iused=0 11170 j=1 goto 11173 11171 j=j+1 11173 if((j).gt.(npro-1))goto 11172 if(axis((j-1)*npred+i) .le. 0.1)goto 11191 iused=1 goto 11172 11191 continue goto 11171 11172 continue if(iused .ne. 1)goto 11211 goto 11141 11211 continue 11161 continue 11220 do 11221 j=1,npred axis((npro-1)*npred+j)=0. 11221 continue 11222 continue axis((npro-1)*npred+i)=1. asr=smuuth(axis((npro-1)*npred+1),npred) if(asr .ge. asrmin)goto 11241 asrmin=asr imin=i 11241 continue if(npro .ne. 1)goto 11261 write(itap1,10090)i,asr if(iprint .lt. 2)goto 11281 call range(npro) call plores 11281 continue 11261 continue 11141 continue 11142 continue 11290 do 11291 i=1,npred axis((npro-1)*npred+i)=0. 11291 continue 11292 continue axis((npro-1)*npred+imin)=1. asr=smuuth(axis((npro-1)*npred+1),npred) 10090 format( 29h0average squared residual for,i3, 13hth predictor=, *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 11300 do 11301 i=1,npro nbspli=nbspli+itnkno(i) 11301 continue 11302 continue nbspli=nbspli-(korder+1)*npro+1 11310 do 11311 i=1,nbspli 11320 do 11321 j=1,nbspli xtx((j-1)*maxpro*maxkno+i)=0. 11321 continue 11322 continue stor4(4*nobs+i)=0. 11311 continue 11312 continue 11330 i=1 goto 11333 11331 i=i+1 11333 if((i).gt.(nobs))goto 11332 istart=0 11340 do 11341 j=1,nbspli stor2(j)=0. 11341 continue 11342 continue 11350 j=1 goto 11353 11351 j=j+1 11353 if((j).gt.(npro))goto 11352 proj=0. 11360 do 11361 k=1,npred proj=proj+axis((j-1)*npred+k)*predic((i-1)*npred+k) 11361 continue 11362 continue 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 11390 do 11391 k=1,korder stor2(is+k)=biatx(k) 11391 continue 11392 continue istart=istart+itnkno(1)-korder goto 11401 11381 continue is=istart+left-(korder+1) if(left .ne. korder)goto 11421 11430 do 11431 k=2,korder stor2(is+k)=biatx(k) 11431 continue 11432 continue goto 11441 11421 continue 11450 do 11451 k=1,korder stor2(is+k)=biatx(k) 11451 continue 11452 continue 11441 continue 11411 continue istart=istart+itnkno(j)-(korder+1) 11401 continue 11371 continue goto 11351 11352 continue 11460 do 11461 j=1,nbspli stor4(4*nobs+j)=stor4(4*nobs+j)+respon(i)*stor2(j) 11470 do 11471 k=1,nbspli xtx((k-1)*maxpro*maxkno+j)=xtx((k-1)*maxpro*maxkno+j)+stor2(j)*sto *r2(k) 11471 continue 11472 continue 11461 continue 11462 continue goto 11331 11332 continue 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 11491 continue if(iprint .ne. 3)goto 11511 write(itap1,10190)(stor4(nobs+i),i=1,nbspli) 10190 format( 19h0eigenvalues of xtx/20(5(1x,g12.6)/)) 11511 continue small=ratio*stor4(nobs+nbspli) ibig=1 11520 continue 11521 if(stor4(nobs+ibig) .ge. small) goto 11522 ibig=ibig+1 goto 11521 11522 continue 11530 do 11531 i=ibig,nbspli stor4(i)=0. 11540 do 11541 j=1,nbspli stor4(i)=stor4(i)+xtx((i-1)*maxpro*maxkno+j)*stor4(4*nobs+j) 11541 continue 11542 continue stor4(i)=stor4(i)/stor4(nobs+i) 11531 continue 11532 continue 11550 do 11551 i=1,nbspli stor4(4*nobs+i)=0. 11560 do 11561 j=ibig,nbspli stor4(4*nobs+i)=stor4(4*nobs+i)+xtx((j-1)*maxpro*maxkno+i)*stor4(j *) 11561 continue 11562 continue 11551 continue 11552 continue is=0 11570 do 11571 i=1,npro it4=itnkno(i)-korder if(i .ne. 1)goto 11591 11600 do 11601 j=1,it4 is=is+1 stor2(2*maxpro*maxkno+(i-1)*maxkno+j)=stor4(4*nobs+is) 11601 continue 11602 continue goto 11611 11591 continue stor2(2*maxpro*maxkno+(i-1)*maxkno+1)=0. 11620 do 11621 j=2,it4 is=is+1 stor2(2*maxpro*maxkno+(i-1)*maxkno+j)=stor4(4*nobs+is) 11621 continue 11622 continue 11611 continue 11581 continue 11571 continue 11572 continue rss=0. 11630 do 11631 i=1,nobs stor4(4*nobs+i)=0. 11640 do 11641 j=1,npro proj=0. 11650 do 11651 k=1,npred proj=proj+axis((j-1)*npred+k)*predic((i-1)*npred+k) 11651 continue 11652 continue 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) 11641 continue 11642 continue stor4(4*nobs+i)=respon(i)-stor4(4*nobs+i) rss=rss+stor4(4*nobs+i)**2 11631 continue 11632 continue asr=rss/float(nobs-nbspli) return 10090 format( 15h1projection nr ,i3, 11h trial nr ,i3/ 33h unab *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 11660 do 11661 i=1,nobs proj=0. itag(i)=i 11670 do 11671 j=1,npred proj=proj+axis((ipro-1)*npred+j)*predic((i-1)*npred+j) 11671 continue 11672 continue stor4(i)=proj 11661 continue 11662 continue call sort(stor4(1),itag,1,nobs) i=0 11680 continue 11681 if(i.ge.nobs-1)goto 11682 i=i+1 m0=i 11690 continue 11691 if(stor4(i+1).gt.stor4(i))goto 11692 i=i+1 if(i.ge.nobs)goto 11692 goto 11691 11692 continue if(i.eq.m0)goto 11681 call ranums(stor4(2*nobs+1),i-m0) 11700 j=m0 goto 11703 11701 j=j+1 11703 if((j).gt.(i-1))goto 11702 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 11702 continue goto 11681 11682 continue avvar=0. ibw2=sambw*bandw*nobs/2.0 sxy=0.0 sx=sxy sy=sx sxx=sy syy=sxx 11710 do 11711 j=1,ibw2 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 11711 continue 11712 continue ibweff=ibw2 11720 do 11721 i=1,nobs 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 11741 continue 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 11761 continue 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 11721 continue 11722 continue ibw2=varbw*bandw*nobs/2.0 sx=0.0 11770 do 11771 j=1,ibw2 sx=sx+stor4(2*nobs+j) 11771 continue 11772 continue ibweff=ibw2 11780 do 11781 i=1,nobs if(i .gt. nobs-ibw2)goto 11801 sx=sx+stor4(2*nobs+i+ibw2) ibweff=ibweff+1 11801 continue if(i .le. ibw2+1)goto 11821 sx=sx-stor4(2*nobs+i-ibw2-1) ibweff=ibweff-1 11821 continue stor4(nobs+i)=sx/ibweff 11781 continue 11782 continue avvar=0. 11830 do 11831 i=1,nobs avvar=avvar+stor4(nobs+i) 11831 continue 11832 continue avvar=avvar/float(nobs) rhigh=avvar*smofac(3) rlow=avvar*smofac(1) avvar=0. 11840 do 11841 i=1,nobs if(stor4(nobs+i) .le. rhigh)goto 11861 stor4(nobs+i)=rhigh goto 11851 11861 if(stor4(nobs+i) .ge. rlow)goto 11871 stor4(nobs+i)=rlow 11871 continue 11851 continue stor4(nobs+i)=1./stor4(nobs+i) avvar=avvar+stor4(nobs+i) 11841 continue 11842 continue 11880 do 11881 i=1,nobs stor4(nobs+i)=stor4(nobs+i)/avvar 11881 continue 11882 continue nkint=maxkno-2*korder 11890 do 11891 i=1,korder stor2(maxkno*maxpro+(ipro-1)*maxkno+i)=stor4(1) 11891 continue 11892 continue nk=korder nsame=korder sum=0. ind=0 11900 do 11901 i=1,nkint 11910 continue 11911 continue ind=ind+1 sum=sum+stor4(nobs+ind) if(sum.ge.float(i)/float(nkint+1))goto 11912 goto 11911 11912 continue if(stor4(ind) .ne. stor2(maxkno*maxpro+(ipro-1)*maxkno+nk))goto 11 *931 nsame=nsame+1 goto 11941 11931 continue nsame=1 11941 continue 11921 continue if(nsame .le. korder)goto 11961 goto 11901 11961 continue nk=nk+1 stor2(maxkno*maxpro+(ipro-1)*maxkno+nk)=stor4(ind) 11901 continue 11902 continue 11970 i=nk goto 11973 11971 i=i+(-1) 11973 if((-1)*((i)-(1)).gt.0)goto 11972 if(stor2(maxkno*maxpro+(ipro-1)*maxkno+i) .eq. stor4(nobs))goto 11 *991 goto 11972 11991 continue nk=nk-1 goto 11971 11972 continue if(nk .ge. korder)goto 12011 write(itap1,10090)(axis((ipro-1)*npred+i),i=1,npred) stop 12011 continue d=(stor4(nobs)-stor4(1))*10.e-6 12020 i=nk+1 goto 12023 12021 i=i+1 12023 if((i).gt.(nk+korder))goto 12022 stor2(maxkno*maxpro+(ipro-1)*maxkno+i)=stor4(nobs)+d goto 12021 12022 continue itnkno(ipro)=nk+korder return 10090 format( 43h0no knot placement possible, only one value / 6h a *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) 12041 continue nit=1 nfun=0 12050 continue 12051 continue xnorm=0. 12060 i=1 goto 12063 12061 i=i+1 12063 if((i).gt.(nvar))goto 12062 xnorm=xnorm+x(i)**2 goto 12061 12062 continue xnorm=sqrt(xnorm) scrat(1,1)=sqrt((1.-x(1)/xnorm)/2.) if(scrat(1,1) .ne. 0.)goto 12081 12090 i=2 goto 12093 12091 i=i+1 12093 if((i).gt.(nvar))goto 12092 scrat(i,1)=0. goto 12091 12092 continue goto 12101 12081 continue 12110 i=2 goto 12113 12111 i=i+1 12113 if((i).gt.(nvar))goto 12112 scrat(i,1)=-x(i)/(2.*xnorm*scrat(1,1)) goto 12111 12112 continue 12101 continue 12071 continue 12120 i=2 goto 12123 12121 i=i+1 12123 if((i).gt.(nvar))goto 12122 12130 j=1 goto 12133 12131 j=j+1 12133 if((j).gt.(nvar))goto 12132 scrat(j,2)=-2.*scrat(i,1)*scrat(j,1) goto 12131 12132 continue scrat(i,2)=scrat(i,2)+1 call cirmin(x,scrat(1,2),nvar,fun,nfun,scrat(1,3)) goto 12121 12122 continue 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) 12151 continue nit=nit+1 if(nit .le. maxit)goto 12171 goto 12052 12171 continue goto 12051 12052 continue 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) 12191 continue if(iprint .le. 0)goto 12211 write(itap2,10320) 12211 continue return 10090 format( 14h0starting axis/10(5(2x,g12.6)/)) 10190 format( 11h0iteration ,i5, 6h nfun=,i5, 7h fmin= ,g12.6, * 15h current axis =/10(5(2x,g12.6)/)) 10320 format( 21h iteration terminated) 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) 12220 i=1 goto 12223 12221 i=i+1 12223 if((i).gt.(nvar))goto 12222 scrat(i)=acos*x(i)+asin*rotvec(i) goto 12221 12222 continue 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 12241 continue 12250 continue 12251 continue idelta=idelta+1 w=v+delta(idelta) if(w .ge. 1.)goto 12271 sw=sign*w acos=cos(sw*pi) asin=sin(sw*pi) 12280 i=1 goto 12283 12281 i=i+1 12283 if((i).gt.(nvar))goto 12282 scrat(i)=acos*x(i)+asin*rotvec(i) goto 12281 12282 continue fw=fun(scrat,nvar) nfun=nfun+1 if(fw .le. fv)goto 12301 goto 12310 12301 continue u=v fu=fv v=w fv=fw goto 12321 12271 continue w=1. fw=f0 goto 12310 12321 continue 12261 continue goto 12251 12252 continue 12310 continue 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) 12330 i=1 goto 12333 12331 i=i+1 12333 if((i).gt.(nvar))goto 12332 scrat(i)=acos*x(i)+asin*rotvec(i) goto 12331 12332 continue fsol=fun(scrat,nvar) nfun=nfun+1 if(fsol .gt. fv)goto 12351 12360 i=1 goto 12363 12361 i=i+1 12363 if((i).gt.(nvar))goto 12362 x(i)=scrat(i) goto 12361 12362 continue return 12351 continue acos=cos(sign*pi*v) asin=sin(sign*pi*v) 12370 i=1 goto 12373 12371 i=i+1 12373 if((i).gt.(nvar))goto 12372 x(i)=acos*x(i)+asin*rotvec(i) goto 12371 12372 continue 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 12380 do 12381 i=1,nobs s=0. 12390 do 12391 j=1,npred s=s+axis((npro-1)*npred+j)*predic((i-1)*npred+j) 12391 continue 12392 continue stor4(i)=s itag(i)=i 12381 continue 12382 continue call sort(stor4(1),itag,1,nobs) xmin=stor4(imin) xmax=stor4(imax) xmax=xmax+1.e-4*(xmax-xmin) rmin=rbig rmax=-rmin 12400 do 12401 i=imin,imax r=stor4(3*nobs+itag(i)) if(r .ge. rmin)goto 12421 rmin=r 12421 continue if(r .le. rmax)goto 12441 rmax=r 12441 continue stor4(nobs+i)=r-stor4(4*nobs+itag(i)) 12401 continue 12402 continue rmax=rmax+1.e-4*(rmax-rmin) call psetup(min0(nxplt,100),xmin,xmax, min0(nyplt,50),rmin,rmax) 12450 do 12451 i=1,nobs call plot(stor4(i),stor4(3*nobs+itag(i))) 12451 continue 12452 continue call pprint write(itap3,10090) 10090 format(1h0/1h0/1h0) 10190 format( 47h1residuals from model of previous iteration and/ *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 12491 continue if(x .lt. stor4(nobs))goto 12511 ufun1=y-stor4(nobs+nobs) return 12511 continue min=1 max=nobs 12520 continue 12521 if(max-min.le.1)goto 12522 nt=(max+min)/2 if(x .ge. stor4(nt))goto 12541 max=nt goto 12551 12541 continue min=nt 12551 continue 12531 continue goto 12521 12522 continue if(stor4(max) .ne. stor4(min))goto 12571 ufun1=y-stor4(nobs+min) return 12571 continue ufun1=y-stor4(nobs+min)-(stor4(nobs+max)-stor4(nobs+min)) *(x-sto *r4(min))/(stor4(max)-stor4(min)) goto 12581 12471 continue xx=x if(xx .ge. xlim(1,iplopr))goto 12601 xx=xlim(1,iplopr) 12601 continue if(xx .le. xlim(2,iplopr))goto 12621 xx=xlim(2,iplopr) 12621 continue r=bvalue(rkno((iplopr-1)*maxkno+1),coeff((iplopr-1)*maxkno+1), nk *no(iplopr)-korder,korder,xx,0) ufun1=r-y 12581 continue 12461 continue 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 12630 do 12631 i=1,npred if(ax(i).ne.0.0) ic=ic+1 12631 continue 12632 continue if(ic .ne. 1)goto 12651 12660 id=1 goto 12663 12661 id=id+1 12663 if((id).gt.(npred))goto 12662 if(ax(id).ne.0.0)goto 12662 goto 12661 12662 continue 12651 continue 12670 i=1 goto 12673 12671 i=i+1 12673 if((i).gt.(nobs))goto 12672 if(ic .ne. 1)goto 12691 stor4(nobs+i)=predic((i-1)*npred+id) goto 12671 12691 continue sum=0. 12700 j=1 goto 12703 12701 j=j+1 12703 if((j).gt.(npred))goto 12702 sum=sum+ax(j)*predic((i-1)*npred+j) goto 12701 12702 continue stor4(nobs+i)=sum goto 12671 12672 continue 12710 i=1 goto 12713 12711 i=i+1 12713 if((i).gt.(nobs))goto 12712 itag(i)=i goto 12711 12712 continue call sort(stor4(nobs+1),itag,1,nobs) i=0 12720 continue 12721 if(i.ge.nobs-1)goto 12722 i=i+1 m0=i 12730 continue 12731 if(stor4(nobs+i+1).gt.stor4(nobs+i))goto 12732 i=i+1 if(i.ge.nobs)goto 12732 goto 12731 12732 continue if(i.eq.m0)goto 12721 call ranums(stor4(2*nobs+1),i-m0) 12740 j=m0 goto 12743 12741 j=j+1 12743 if((j).gt.(i-1))goto 12742 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 12742 continue goto 12721 12722 continue 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) 12750 do 12751 j=1,ibw2 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 12751 continue 12752 continue ibweff=ibw2 12760 do 12761 i=1,nobs 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 12781 continue 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 12801 continue 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 12761 continue 12762 continue avvar=avvar/float(nobs) ibw2=varbw*bandw*nobs/2.0 sx(1)=0.0 12810 do 12811 j=1,ibw2 sx(1)=sx(1)+stor4(4*nobs+j) 12811 continue 12812 continue ibweff=ibw2 12820 do 12821 i=1,nobs if(i .gt. nobs-ibw2)goto 12841 sx(1)=sx(1)+stor4(4*nobs+i+ibw2) ibweff=ibweff+1 12841 continue if(i .le. ibw2+1)goto 12861 sx(1)=sx(1)-stor4(4*nobs+i-ibw2-1) ibweff=ibweff-1 12861 continue stor4(2*nobs+i)=sx(1)/ibweff 12821 continue 12822 continue 12870 do 12871 i=1,3 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 12871 continue 12872 continue 12880 j=2 goto 12883 12881 j=j+1 12883 if((j).gt.(ibw(3)))goto 12882 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 12882 continue 12890 do 12891 i=1,3 ibwt(i)=ibw(i)-1 12891 continue 12892 continue 12900 do 12901 i=1,nobs 12910 do 12911 j=1,3 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 12931 continue 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 12951 continue 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 12971 continue 12911 continue 12912 continue 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 12991 continue 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) 13001 continue 12981 continue if(thisfc .gt. smofac(1))goto 13021 stor4(4*nobs+i)=s1 goto 12901 13021 continue if(thisfc .lt. smofac(3))goto 13041 stor4(4*nobs+i)=s3 goto 12901 13041 continue 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 13061 continue wt=(thisfc-smofac(2))/(smofac(3)-smofac(2)) stor4(4*nobs+i)=(1.0-wt)*s2+wt*s3 13071 continue 13051 continue 12901 continue 12902 continue 13080 do 13081 i=1,nobs stor4(nobs+itag(i))=stor4(4*nobs+i) 13081 continue 13082 continue avsr=0. 13090 do 13091 i=1,nobs r=stor4(3*nobs+i)-stor4(nobs+i) avsr=avsr+r**2 stor4(4*nobs+i)=r 13091 continue 13092 continue 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 1000 CONTINUE 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